Diversity of insect-specific viruses in Culex mosquitoes from Papua New Guinea

Mosquito samples

Mosquitoes were collected between March 2019 and August 2021 from five provinces in PNG (Figure 1, Supplementary Table 1a). Collection dates varied between provinces. Thirteen collection sites were selected near piggeries and poultry farms, taking into account typical mosquito habitats in the local area. Mosquitoes were collected using CDC Light traps (normal and UV), BG Sentinel traps (with and without lure), VioLEDs-UV traps, and the human landing catch method [39, 40]. Collections were conducted every 24 h, and mosquitoes were transported to the laboratory for speciation using morphological features. Following identification, mosquitoes were surface-sterilized by rinsing in 70% ethanol and pooled (up to n=10) by location, species, and feeding status (blood-fed or unfed, determined by visual inspection of abdominal engorgement). Each mosquito pool was placed into 100 μL RNAlater (Thermo Fisher Scientific) and stored at until further testing.

Homogenization and RNA extraction

Mosquito samples preserved in RNAlater solution were thawed at room temperature and transferred to a sterile Petri dish within a Class II biosafety cabinet. Individual mosquito samples were prepared by placing each mosquito into a 2 mL microcentrifuge tube containing 250 μL of silicon carbide shards and 600 μL of Qiagen RNeasy RLT buffer. Secondary pools of up to 500 mosquitoes were prepared by first making sub-pools of no more than 100 mosquitoes each belonging to the same collection location, dates, mosquito species and feeding status. Sub-pools were created by adding up to 100 mosquitoes into a 2 mL microcentrifuge tube containing 250 μL of silicon carbide shards and 1 mL of Qiagen RNeasy RLT buffer.

Samples were homogenized using a MP FastPrep-24 5G bead beater at 4.0 m s for two cycles of 20 s each. Samples were kept on ice for 30 s between cycles. Tubes were centrifuged at for 5 min, and supernatants from each sub-pool were collected into separate tubes for nucleic acid extraction. Total RNA was extracted using the Qiagen RNeasy Mini Kit according to the manufacturer’s instructions. Eluates from the sub-pools were combined to form each secondary pool.

A total of 86 study samples, consisting of 22 individual mosquito samples and 64 secondary pooled samples, were obtained for sequencing. Total RNA was quality-checked for concentration (Qubit RNA Assay Kit, Thermo Fisher Scientific) and purity (NanoDrop, Thermo Fisher Scientific).

Shotgun metagenomic sequencing

Sequencing libraries were prepared using the Tecan Universal Plus Total RNA-Seq with NuQuant library preparation kit, incorporating mosquito ribosomal RNA depletion probes, according to the manufacturer’s instructions. Library quality was assessed using the Agilent 2100 Bioanalyzer, and DNA concentration was measured using the Invitrogen dsDNA High Sensitivity Assay Kit on a Qubit fluorometer. Libraries were pooled and sequenced on the Illumina NextSeq 2000 platform using a P2 cartridge with 200/300 cycles.

Four NextSeq runs were performed, each comprising 24–28 samples and controls. Run performance was evaluated using Illumina internal quality metrics, including cluster density, percentage of clusters passing filter, and Q30 scores, with all runs meeting recommended specifications.

Targeted PCR amplification of vertebrate mitochondrial cytochrome oxidase I gene from blood-fed mosquitoes

When the vertebrate host identity of blood meals could not be resolved from unbiased sequencing data, targeted PCR amplification was performed using a previously described protocol employing vertebrate-specific primers targeting the mitochondrial cytochrome oxidase I (COI) gene [41, 42]. PCRs were conducted following the protocol and cycling conditions described by Reeves and Burkett-Cadena (2024), with the sole modification that reactions in this study used 2× Platinum Hot Start PCR Master Mix (Thermo Fisher Scientific). Amplicons were subsequently purified using the QIAQuick gel purification system (Qiagen) and sequenced with the BigDye terminator v3.1 kit on a 3500 Genetic Analyser (Applied Biosystems).

Bioinformatic analysis

Raw reads from the Illumina sequencing runs were processed and assembled using the MetaDIVE metagenomics pipeline (https://github.com/James-ODwyer/MetaDIVE/tree/master). In summary, raw reads were trimmed, quality-filtered, and adapters were removed using fastp (v0.23.2), with the following non-default settings: minimum read length of 85 bp, minimum base Phred quality score of 16, and a four-base sliding window with a mean quality threshold of 20 [43]. Quality control of filtered reads was performed using FastQC (v0.12.1) [44].

Reads were assembled de novo into contigs using Trinity (v2.8.5), with a minimum contig length of 301 bp [45]. Contigs were translated, and all reading frames were searched against the NCBI non-redundant protein database (June 2024 release) using DIAMOND (v2.0.15) BLASTx under sensitive settings [46], with the species of the top hit assigned as the most likely identity. Viral contigs were further confirmed by BLASTn against the NCBI nucleotide database [47]. Taxonomic assignments were determined using the R package taxonomizr (v0.7.1) based on the top-matching species from DIAMOND BLASTx. Filtered reads were then mapped back to contigs to generate read counts assigned to each species using Bowtie2 v2.4.5 [48] with the preset –very-fast setting, to require only high-identity matches. To recover additional reads associated with detected viral species, reads that did not assemble into contigs were searched directly against viral protein sequences using DIAMOND BLASTx. Searches were restricted to viral species for which at least one contig had already been identified, with the preset set to –fast and –max-target-seqs set to 1. Matching reads were clustered at 97% sequence similarity using CD-HIT v4.8.1 [49], and representative sequences were confirmed by BLASTn. Taxonomic assignments were determined using the R package taxonomizr, and cluster-level assignments were transferred to individual reads within each cluster. Read counts from unassembled viral reads were combined with counts from reads mapped to viral contigs to obtain the total number of reads assigned to each viral species.

Viral sequences identified in the study were classified at the species level using International Committee on Taxonomy of Viruses (ICTV) demarcation criteria where available (https://ictv.global/report). In the absence of formal thresholds, sequences were assigned based on amino acid identity or nucleotide identity, consistent with previous studies [50, 51]. For community composition analyses, only viral reads showing ≥85% nucleotide identity to reference sequences were retained for species-level assignment and downstream abundance and diversity analyses. However, low-identity viral sequences associated with substantial assembled contigs (>500 bp and supported by >100 reads) were retained for further investigation, including manual BLAST analyses, genome characterization and, where appropriate, segment association analyses. Remaining low-identity viral sequences were not investigated further in this study. Genome or segment completeness was assessed following NCBI Viral Genome Submission guidelines. Sequences were classified based on the presence of intact coding regions with appropriate start and stop codons, absence of assembly gaps or ambiguous bases, and sequence lengths consistent with the corresponding reference genome. Sequences that met these criteria were considered complete, whereas incomplete coding regions or substantial genomic gaps were classified as partial genomes or segments (https://www.ncbi.nlm.nih.gov/labs/virus/vssi/docs/help/#glossary-nuc-completeness). Viruses were classified as insect-specific based on phylogenetic placement within well-supported ISV clades and highest similarity to reference ISVs derived from invertebrates.

Generated viral genomes were imported into Geneious Prime (v2020.1.1) for manual inspection and refinement [52]. Diversity analyses were conducted in RStudio (v2024.04.2) [53] using the tidyverse (v2.0.0), vegan (v2.6–10), ggplot2 (v3.5.0), and phytools (v2.4-4) packages, alongside Microsoft Excel [54]. Alpha diversity was quantified using the Shannon, Simpson, and Chao1 indices [55–57]. Differences in alpha diversity among locations, mosquito species, and feeding status were assessed using the non-parametric Kruskal–Wallis test [58].

Reference sequences for phylogenetic analyses were retrieved from GenBank (June 2024 release) based on BLAST similarity to PNG-derived viral sequences. The closest matching sequences and representative viruses from major phylogenetic groups were selected to provide context for the placement of PNG-derived sequences. ICTV reference sequences were included where available. Multiple sequence alignments were generated using MAFFT (v7) with automatic model selection [59]. Alignment trimming was performed using TrimAl (v1.4.1) under default settings [60]. Phylogenetic analyses were conducted using IQ-TREE (v2.0) and MEGA (v11), with maximum likelihood inference, automatic model selection, and 1,000 bootstrap replicates [61, 62]. Analyses were based on amino acid sequences of RNA-dependent RNA polymerase (RdRP), capsid, or polyprotein regions. Tanglegrams were constructed in RStudio using the phytools package (v2.5-2) to compare relationships between corresponding phylogenetic trees and visualize associations among related viral sequences. Figures were refined using Inkscape [63], with additional visualization performed in Microreact and iTOL [64, 65].

Vertebrate host species blood meal sources were identified using mitochondrial cytochrome oxidase I (COI), large subunit (LSU), and small subunit (SSU) rRNA gene sequences. Filtered contigs >500 bp with identity to Chordata were retained. The longest contigs were selected and confirmed by BLASTn against the NCBI nucleotide database. Sanger sequencing data for the CO1 gene were quality-checked, assembled, and analyzed using DNASTAR Lasergene SeqMan Pro v18 and Geneious Prime v.2020.1.1 for the vertebrate host identification, when not resolved by unbiased sequencing data.

Associations between vertebrate blood-meal source and virus detection were assessed using Fisher’s exact tests [66], with false discovery rate correction applied using the Benjamini–Hochberg method [67]. Differences in overall virus community composition among host categories were evaluated using permutational multivariate analysis of variance (PERMANOVA) on Jaccard distances derived from presence–absence matrices, with 999 permutations [68, 69].

Data deposition

All viral genome sequences generated in this study have been submitted to GenBank under accession numbers PZ253856–PZ254055 and will be publicly available upon completion of GenBank processing. Accession numbers are provided in Supplementary Table 10. Raw sequencing reads from the Illumina NextSeq runs have been deposited in the Sequence Read Archive (SRA) under BioProject accession PRJNA1444521.

Mosquito sample description

A total of 13,871 female Culex mosquitoes were analyzed, with the highest specimen abundance recorded in the Madang and Morobe provinces during the sampling period. Cx. quinquefasciatus was the most prevalent species, accounting for 50% of the total specimens. The sample set comprised approximately equal proportions of unfed and blood-fed mosquitoes, ensuring a representative distribution across different feeding statuses (Table 1).

Shotgun metagenomic sequencing

The virome of 13,871 mosquitoes across 86 samples was characterized using unbiased high-throughput sequencing. A detailed overview of sample characteristics is provided in Supplementary Table 1b. Sequencing on the Illumina NextSeq platform generated approximately 1.6 billion paired-end reads, with read lengths of either 100 or 150 base pairs. These reads were evenly distributed across four sequencing libraries. Following stringent quality filtering, 20,158,310 paired-end reads were classified as viral, accounting for ∼1.5% of the total sequencing output. The proportion of viral reads varied across samples, ranging from 0.002% to 18.118% of the total filtered reads per sample (Supplementary Table 1b). Of these viral read pairs, 10,180,845 (50.5%) met the species-assignment threshold ( nucleotide identity to a reference sequence) and were retained for downstream community analyses. The remaining 9,977,465 viral read pairs (49.5%) showed <85% nucleotide identity to reference sequences.

ISV diversity of classified viruses

A total of 231 viral species representing 38 virus families were identified. Of these, 120 were ISV species spanning 25 virus families. ISVs accounted for 93.8% of species-assigned viral reads, demonstrating that they comprised the dominant component of the mosquito virome (Table 2). The majority of identified ISV species were RNA viruses (87%), while a smaller proportion were insect-specific DNA viruses (13%). The most abundant families were Partitiviridae, Birnaviridae, Rhabdoviridae, and Iflaviridae, collectively accounting for 86% of all classified viral reads. In addition, approximately 40% of viral reads corresponded to viruses that are not currently assigned to recognized ICTV virus families (Figure 2a). No single ISV family was detected across all samples, underscoring the high compositional diversity of the ISV community. At the species level, Hubei virga-like virus 2 (HVLV2) was the most abundant ISV, followed by species belonging to the families Birnaviridae, Partitiviridae, and Rhabdoviridae (Figure 2b).

Alpha diversity patterns

Alpha diversity was assessed using the Shannon, Simpson, and Chao1 indices (Supplementary Table 1c). No significant differences were detected among provinces or mosquito species for any of the three metrics, and Shannon and Simpson diversity did not differ significantly between blood-fed and unfed mosquitoes (Kruskal–Wallis tests, p > 0.05). However, Chao1 richness was significantly higher in unfed mosquitoes than in blood-fed mosquitoes (p = 0.0003). Overall, these results indicate broadly similar ISV diversity across locations and mosquito species, although estimated richness differed between feeding categories.

Phylogenetic analysis of selected ISVs

Specific ISV groups that demonstrated higher abundance and had sufficient genomic sequence data were selected for phylogenetic analysis.

Blood-meal host identity and virus group distribution

Molecular identification of vertebrate blood-meal sources from blood-fed mosquitoes (n = 63) indicated that pigs were the most common host (31.7%), followed by chickens (27%), other non-human vertebrates (12.7%), humans (9.5%), and mixed blood meals (14.3%), while host identity could not be resolved in 4.8% of specimens (Supplementary Figure 6).

The distribution of vertebrate hosts varied among virus groups. Cambodia Culex birnavirus, rhabdoviruses, iflaviruses, Hubei virga-like virus 2, and Culex-associated luteo-like virus were most frequently detected in chicken-fed mosquitoes, accounting for 46–67% of detections of these virus groups. In contrast, Broome luteo-like virus 1, putative Broome partiti-like virus 1, and Culex tritaeniorhynchus rhabdovirus showed a higher relative contribution from pig-derived blood meals (33.3–42.9%). Human blood meals contributed a smaller but consistent proportion of detections across all studied virus groups (∼5–29%).

Among host-identified blood-fed mosquitoes, HVLV2 was the only virus that retained a significant association with vertebrate host category after false discovery rate (FDR) correction (q = 0.0167). HVLV2 was most frequently detected in chicken-fed mosquitoes (66.7%) and was absent from pig-fed mosquitoes, with intermediate detection in other host groups. No significant differences in overall virus community composition were detected among vertebrate host categories (PERMANOVA, , p = 0.136).

Alpha diversity patterns

Alpha diversity analyses indicated that ISV community diversity was broadly similar across sampling locations and mosquito species, with no significant differences detected using the Shannon, Simpson, or Chao1 indices. This pattern is consistent with previous studies reporting broadly similar ISV diversity across mosquito populations and geographic locations [79–81]. In contrast, unfed mosquitoes exhibited significantly higher richness than blood-fed mosquitoes based on the Chao1 index, indicating a greater estimated richness of ISV taxa. However, no corresponding differences were detected using the Shannon or Simpson indices, suggesting that overall community diversity and dominance patterns remained broadly similar. As Chao1 specifically estimates richness, the biological significance of this pattern warrants further investigation. Nevertheless, the absence of significant differences in Shannon and Simpson diversity is consistent with previous studies reporting little or no influence of recent blood feeding on insect-specific viral communities [83, 84]. Additionally, unequal sampling between locations and mosquito species, including very small sample sizes for some categories, may have influenced diversity estimates and should be considered when interpreting these results. Differences in sampling times among locations may also have influenced mosquito abundance and composition through seasonal variation, potentially contributing to some of the observed variation in diversity estimates.

Novel partitivirus genome segments

A comparison of Partitiviridae family RdRP and capsid sequences from the same pooled samples generally exhibited congruent clustering. Notably, all PNG-derived sequences that clustered with Broome partiti-like virus 1 in the RdRP phylogeny had corresponding capsid sequences clustering with Aedes partiti-like virus 1 in the capsid phylogeny. Apart from three exceptions, PNG sequences clustering with Broome partiti-like virus 1 exhibited greater than 90% amino acid identity with the reference sequences, suggesting that they represent the same species [70]. The three exceptions showed less than 90% identity to the reference Broome partiti-like virus 1 and may represent divergent strains of the same species or a novel yet closely related virus, given that they clustered within the same clade (Supplementary Figure 1a). Furthermore, their reciprocal capsid sequences shared only amino acid identity with their closest relative, Aedes partiti-like virus 1, yet were more than 99% identical to each other, indicating that they originate from a distinct virus (Supplementary Figure 1b).

To date, Broome partiti-like virus 1 has only been identified and characterized using the RdRP gene [74, 85]. Based on the strict co-occurrence of these segments across samples and their low sequence similarity to all currently described viruses, we propose that the newly identified capsid segments represent the missing genomic components of Broome partiti-like virus 1, thereby revealing the complete set of genome segments of this ISV for the first time.

A similar pattern was observed for Sonnbo virus, for which only the RdRP segment has been reported to date [77, 86, 87]. In this study, a Sonnbo virus RdRP segment was detected in one pooled PNG mosquito sample (23–00945-0037), showing 99.7% amino acid identity to the reference sequence, confirming the presence of this virus. However, in the same sample, an additional genome segment encoding a putative capsid protein was recovered, whose closest available match was Zeya Brook partiti-like virus 1, sharing only amino acid similarity despite high sequence coverage. The strict co-occurrence of this capsid segment with the Sonnbo virus RdRP, together with its absence from all other samples, strongly suggests that it represents the previously undescribed capsid segment of Sonnbo virus. We therefore propose that this sequence constitutes the missing capsid genomic component of Sonnbo virus, extending its genomic characterization beyond the RdRP segment for the first time. Together with the findings for Broome partiti-like virus 1, this observation highlights a recurring pattern in partitivirus discovery, whereby RdRP segments are detected and classified in the absence of their corresponding capsid segments, which may remain unrecognized due to their extreme sequence divergence.

In addition to these paired segment discoveries, several highly divergent RdRP and capsid sequences were detected without identifiable corresponding genomic partners (uP;23–00945-0021, Ps;23–00945-0035, uDPV;23–00945-0041, uP;23–00945-0160, uP;23–00945-0164; (Figure 3)). The absence of paired segments may reflect low viral abundance, uneven sequencing depth, or extreme sequence divergence beyond the limits of similarity-based approaches [88, 89].

However, the interpretation of segment associations is limited by the use of pooled mosquito samples in this study [90]. Although consistent co-occurrence supports inferred segment pairing, it cannot be excluded that divergent RdRP and capsid sequences originated from different infected individuals within the same pool.

Cambodia Culex birnavirus is widespread in PNG Culex mosquitoes

In this study, the only member of the family Birnaviridae detected was Cambodia Culex birnavirus. This virus was originally described from Cx. quinquefasciatus mosquitoes collected in Cambodia by Mohamed Ali et al. (2023), and prior to our study, these sequences were the only Cambodia Culex birnavirus entries available in GenBank. Phylogenetic analyses showed that RdRP protein sequences from PNG viruses cluster with Cambodian reference sequences in a distinct monophyletic clade closely related to the genus Entomobirnavirus, yet share only 38.2–43.0% amino acid identity with representative Entomobirnavirus sequences. This supports the previous proposal that Cambodia Culex birnavirus represents a novel genus within the Birnaviridae family [91]. Moreover, the distinct phylogenetic relationships of the PNG strains relative to the Cambodian reference sequences suggest some degree of geographic divergence (Supplementary Figure 2). While no strongly supported clades were identified, the presence of basal clustering by sampling location at the provincial level indicates possible ongoing viral diversification and evolutionary dynamics within the region, as has been reported for ISVs in mosquitoes [92–94].

Cambodia Culex birnavirus was detected across all sampled locations and mosquito species, reinforcing its prior identification as a Culex-associated virus [91]. Notably, our study expands the range of mosquito species in which this virus has been detected to include Cx. annulirostris and Cx. sitiens, in addition to Cx. quinquefasciatus. Although present across all provinces, it was most frequently detected in Madang and Morobe, suggesting potential regional hotspots for viral maintenance or transmission. The higher detection rates in these areas may reflect localized ecological or host-related factors that favour viral persistence [95, 96]. These findings raise the possibility that Cambodia Culex birnavirus may constitute a stable component of the mosquito virome in certain regions [38, 86]. However, lower detection frequencies in the East Sepik, Central, and Western provinces may partly reflect comparatively fewer mosquito collections from these regions than in Madang and Morobe.

Distribution and host associations of PNG mosquito rhabdoviruses

Phylogenetic analysis of rhabdovirus sequences revealed that the PNG isolates clustered tightly with MERDV and GCRV, exhibiting minimal genetic divergence and little evidence of clade diversification relative to these reference lineages. This observation, combined with the limited distribution of the viruses (observed only in samples from the Western and Central Provinces), could indicate either a recent introduction of these viruses into PNG or strong purifying selection maintaining sequence conservation in functionally critical genes [97–99]. The Western and Central Provinces lie on the southern side of the PNG Highlands, which bisect the island and serve as a substantial natural barrier to mosquito dispersal and viral gene flow, separating them from study sites in the other three provinces [100, 101]. This topographic isolation would likely have restricted the geographic spread of MERDV and GCRV viruses, supported by phylogenetic evidence for their recent arrival or constrained evolution within PNG. However, similar geographic structuring was not detected for other viruses identified in this study.

The PNG sequences that clustered with GCRV exhibited comparable levels of nucleotide and amino acid identity to both the GCRV and GMRV1 reference sequences. Although these two virus species were originally described as separate species, they fulfil the ICTV criteria to be considered the same species since the pairwise alignments of their nucleoprotein (N), glycoprotein (G), and RdRP (L) proteins demonstrated greater than 99% similarity [72, 90, 102].

To date, these viruses have been detected exclusively in mosquitoes [50, 87, 90, 91, 102–104]. Most have been identified in the genus Culex, including the PNG sequences described in this study. A limited number of MERDV detections have occurred in the mosquito genera Ochlerotatus and Haemagogus, while GMRV1 has been reported in the genus Deinocerites, all of which belong to the family Culicidae [90, 102, 103]. These observations suggest a strong host preference for Culex mosquitoes while also indicating the capacity of these viruses to infect a broader range of Culicidae genera, underscoring the need to broaden vector surveillance beyond the dominant Culex species.

CTRV was first isolated from Cx. tritaeniorhynchus mosquitoes collected in Chiba, Japan [105]. All GenBank records to date either specify the host as Cx. tritaeniorhynchus or provide the host annotation as ‘Culex species’. In the present study, CTRV-related sequences were detected exclusively in Cx. annulirostris pools. Cx. annulirostris is a dominant culicine mosquito species across Australia and PNG, where it drives the transmission of multiple medically important arboviruses and can occur at very high densities [106, 107]. The restriction of CTRV to this abundant host, combined with its absence from other sympatric Culex species in our dataset, suggests a possible host-switch followed by local adaptation to Cx. annulirostris.

Furthermore, of the five recognized genes in Rhabdoviridae viruses, only the nucleoprotein (N) gene of CTRV was detected in PNG samples. This partial genome recovery may indicate degradation of viral RNA [108, 109], low viral load [110, 111], and/or incomplete genome coverage due to limitations in sequencing depth [112, 113]. Alternatively, it may reflect a divergent or reassorted viral lineage in which only the N gene shares detectable homology with reference sequences, while the remaining segments have diverged beyond current database recognition. The exclusive detection of the N gene in PNG samples may also suggest localized circulation with persistent, low-level infection, where the full virus genome may be present but only some genes are expressed at detectable levels, as observed for RNA viruses under persistent or low-abundance infection conditions [22, 111, 114]. It is also possible that some CTRV-derived sequences represent endogenous viral elements (EVEs) integrated within the mosquito genome [115]. Although the available contigs did not provide sufficient evidence to evaluate this hypothesis, it cannot be excluded based on the current data. Further investigation through targeted enrichment, full-genome sequencing, and virus isolation is required to determine whether this represents a novel rhabdovirus variant, EVE or a partial remnant of a more complete genome.

Expanded geographic range of mosquito-associated iflaviruses

The family Iflaviridae consists of single-stranded, positive-sense RNA viruses found in arthropods, with insects being the most common host [73]. The iflavirus species identified in this study, while previously documented in Culex mosquitoes, have also been detected in Aedes mosquitoes [50, 75, 86, 116–119]. This suggests a potentially wider host range and emphasizes the need for continued research into the distribution and genetic diversity of iflaviruses across different mosquito species.

To our knowledge, this represents the second published report of ICIFV since its original discovery in Japan [120]. Although an additional ICIFV sequence from China has subsequently been deposited in GenBank (PQ442367), the present study reports the first published detection of this virus outside Asia, extending its known geographic range into the Australasian region. In addition, CILV4, previously reported from America, Europe, and Africa [50, 75, 86, 116, 117, 119], and YIV1, previously detected across East Asia, Europe, and Mediterranean regions [118–120], were identified for the first time in PNG.

Phylogenetically, CILV4 and YIV1 each clustered within well-supported monophyletic clades composed of viruses from Dipteran hosts, suggesting a potential host-associated signal within these lineages. However, previous studies have shown that iflaviruses form multiple distinct clades that do not consistently reflect host taxonomy, as individual clades often include viruses from multiple insect orders [73]. Taken together, our findings are consistent with a possible Dipteran-associated signal for these PNG viruses, although broader sampling and the recovery of additional complete genomes will be required to determine whether this pattern reflects genuine host specificity or instead resulted from limited sequence availability.

Regional structuring of HVLV2

HVLV2 is an RNA virus associated with Culex mosquitoes and classified within the realm Riboviria, with its exact taxonomic status remaining unresolved. Initially identified through metagenomic analyses in Hubei Province, China [22], subsequent research has revealed its presence across various geographic locations, indicating a wider distribution than initially suspected [50, 100, 103, 121]. Consistent detections in Culex populations suggest a specific association with this mosquito genus. However, the pathogenic potential of HVLV2 in both its mosquito hosts and any potential vertebrate hosts remains uncharacterized. Notably, our analyses revealed a distinct clustering of PNG sequences by geographic origin, indicating regional genetic divergence of HVLV2 strains, a pattern also observed for Cambodia Culex birnavirus in this study. These observations underscore the need for further investigation into the evolutionary dynamics of HVLV2 and its patterns of maintenance, transmission, and interaction within local mosquito populations.

HVLV2 was detected across all mosquito samples tested, with the exception of Cx. pullus mosquitoes and samples collected from East Sepik Province. However, the absence of detections in East Sepik and Cx. pullus should be interpreted cautiously, as both were represented by comparatively limited sampling, and additional collections will be required to determine whether these patterns reflect true absence or insufficient sampling effort. HVLV2-positive samples were characterized by a high abundance of viral reads, suggesting efficient viral replication under favourable host or environmental conditions. This elevated viral load is consistent with active replication and may indicate a capacity for efficient transmission within local mosquito populations [122, 123]. The pronounced viral read abundance also raises the possibility that HVLV2 plays a biologically significant role in shaping the mosquito virome, potentially influencing host fitness or interactions with co-infecting viruses [124–127]. Further investigations are needed to determine whether this virus represents a persistent infection in specific mosquito lineages or an episodic infection with high replication efficiency.

Frequent detection of BLLV1

BLLV1 was initially detected in Cx. annulirostris mosquitoes from the Kimberley region of Western Australia and subsequently reported in the same mosquito species from Wyndham, Australia [74, 85]. In PNG, the virus was identified across all examined Culex species and sampling locations, indicating a broader host range and more extensive geographical distribution than previously recognized. This widespread detection raises the possibility that BLLV1 may constitute a stable component of the virome of Culex mosquitoes in the Australasian region.

Phylogenetic analysis revealed that PNG BLLV1 sequences formed a well-supported clade when compared only with BLLV1 reference sequences (Figure 8b). However, when additional related luteo-like viruses were included, this clustering was disrupted, with Australian BLLV1 sequences separating the PNG viruses into two distinct groups (Figure 8a). This breakdown of monophyly may reflect the influence of highly divergent sequences on tree topology but could also indicate underlying evolutionary processes such as shared ancestral diversity, incomplete lineage sorting, or multiple introductions across regions [128, 129]. Collectively, these findings underscore the importance of broad yet balanced taxon sampling in phylogenetic analyses, as the inclusion of divergent lineages can substantially alter clade resolution and interpretation [130, 131].

Genomic organization and phylogenetic placement of Culex-associated Luteo-like virus

The PNG Culex-associated Luteo-like virus identified in this study has a bipartite RNA genome, with the RNA-1 segment (∼ 2800 bp) encoding a serine protease and RdRP, while RNA-2 (∼ 1600 bp) encodes a capsid protein and a hypothetical protein (Supplementary Figure 7). Although this virus has been referred to as "Luteo-like" and grouped with members of the former Luteoviridae, (now classified within the family Tombusviridae) [50, 75, 87, 117, 132], its genome organization appears to differ from that of canonical luteoviruses. Members of the genus Luteovirus are monopartite plant viruses with a single non-segmented RNA genome. In contrast, the only recognized bipartite viruses within Tombusviridae belong to the genus Dianthovirus. However, even in dianthoviruses, genome organization differs from that observed here, RNA-1 is larger (∼ 3.9 kb) and encodes both replication-associated proteins and the coat protein, while RNA-2 encodes a movement-associated protein [133]. In the Culex-associated Luteo-like virus, the RdRP is encoded on a shorter RNA-1 segment that does not encode the capsid protein, which is instead encoded on a separate RNA-2 segment. This segmentation pattern suggests a distinct genomic organization compared to both luteoviruses and dianthoviruses.

At the time of writing, only three sequences of Culex-associated Luteo-like virus in the NCBI database were assigned within the order Tolivirales (NCBI nucleotide database search for Culex-associated Luteo-like virus). This classification is reinforced by the current RdRP protein phylogenetic analysis, indicating that Culex-associated Luteo-like virus forms a distinct clade within the Tolivirales lineage. Taken together, its genomic features and phylogenetic placement suggest that Culex-associated Luteo-like virus may represent a divergent lineage within this order, although additional related sequences will be required to better resolve its broader taxonomic context.

Effects of blood-meal host on viral abundance

The predominance of pig and chicken blood meals likely reflects local host availability rather than mosquito feeding preference, as traps were placed near pig and poultry farms, a factor known to bias blood-meal profiles toward nearby hosts [134, 135]. Accordingly, vertebrate blood-meal associations should be interpreted as indicators of mosquito feeding and transient exposure, not as evidence of vertebrate infection or host competence. Several viruses detected in this study, including Cambodia Culex birnavirus, iflaviruses, and partiti-like viruses, are recognized as insect-specific and lack evidence of vertebrate replication [73, 136, 137]. Nevertheless, documenting vertebrate blood-meal sources may be relevant for RNA viruses with high evolutionary plasticity, in which repeated exposure at vector–host interfaces may facilitate host-range expansion under appropriate ecological conditions [138, 139].

Within blood-fed mosquitoes, host-related structuring of the virome was minimal and largely virus-specific. Only HVLV2 showed significant enrichment in chicken-fed mosquitoes, while overall virus community composition did not differ among host categories. While formal statistical evaluations of viral diversity in relation to vertebrate blood-meal host remain relatively uncommon, metagenomic studies combining blood-meal identification with virome analysis have reported descriptive host–virus associations in blood-fed mosquitoes [51, 140, 141]. However, large-scale virome surveys consistently indicate that broad physiological and ecological factors tend to dominate patterns of viral abundance, rather than host-driven restructuring of the virus community [142–144]. Together, these findings suggest that vertebrate blood-meal source exerts limited influence on overall virome structure, with host effects confined to specific viruses rather than reflecting broad restructuring of the mosquito virome. The detection of human blood meals should also be interpreted cautiously, as human DNA may arise from field or laboratory contamination and can be preferentially amplified in molecular blood-meal assays [42, 145].

All viruses detected in this study were present in both blood-fed and unfed mosquitoes, indicating that their presence is not dependent on recent blood feeding and that they are unlikely to be transiently acquired through vertebrate blood meals. Even in cases where unfed mosquitoes may have previously fed and digested the blood meal, the consistent detection of viral sequences suggests active replication within the mosquito host [146]. This further indicates that these viruses are compatible with insect hosts and may represent stable constituents of the microbiome of PNG Culex mosquitoes. Their detection, irrespective of feeding status, also suggests the possibility of vertical or horizontal transmission pathways rather than acquisition through vertebrate hosts [147].

Potential co-infection versus pooling effects

Another interesting observation from this study is that several samples contained reads from more than one species within the same viral family. For example, several pooled samples contained sequences corresponding to two distinct rhabdoviruses or multiple partiti-like viruses detected within the same sample. This pattern may reflect concurrent infections with multiple ISVs within individual mosquitoes. Alternatively, given the large pool sizes (up to 500 mosquitoes per pool) used in this study, it is plausible that these detections reflect distinct viral species present in different individuals within the same pool rather than true co-infections. As most specimens were analyzed as pooled samples, definitive discrimination between these scenarios is not possible without analyses at the individual mosquito level [50, 148].

Ethics approval and consent to participate

Samples used in this study that were collected using the Human Landing Catch (HLC) methodology were obtained under the project Asia-Pacific ICEMR: Epidemiological and Entomological Malaria Surveillance, approved by the Papua New Guinea Medical Research Advisory Committee (MRAC; ethics approval number 18.16).

Consent for publication

Not applicable.

Competing interests

The authors declare no competing interests.