RNA Sequencing data from field caught Haemagogus mosquitoes. All mosquitoes from the twelve forested locations included in the RNA sequencing analysis were identified as Hag. janthinomys based on standard taxonomic keys 23–25 and wing geometric morphology26. However, a subsequent barcoding study conducted on Hag. janthinomys from Trinidad showed that the population is a species complex with three distinct genotypes (Ali et al., unpublished data).
Total RNA sequencing of the Haemagagus mosquito pools from the 12 sites yielded a total of 857,590,224 paired end reads after quality control processing, with an overall mean of 71,465,852 and a range of 57,884,090 to 87,768,438 reads as shown in Table.1. Direct Bowtie2 mapping of the total RNA sequences onto representative nucleic acid databases showed average mapping rates of 2415.93 reads per million (RPM) per location against reference viral sequences, 4.60 RPM against potential mosquito feeding host genomes and 1211.60 RPM against the Silva reference 16S bacterial database. The sites with the highest viral species richness were the Ecclesville and Mamon Village locations which included alphaviruses, insect-specific, mosquito-associated, and other viral types (Fig. 1); with remaining collection sites having lower viral species richness (Fig. S2).
Table 1
RNA read mapping against a reference viral RVDB database, ARB/SILVA 16S rRNA gene sequences and potential animal feeding host reference partial or complete mitochondrial genomes.
|
|
|
Mapped reads per million
(RPM)
|
|
|
|
|
Location
|
Total Number of Reads
|
Viral
|
Bacterial 16S
|
Feeding hosts
|
1
|
Chaguaramas
|
87 768 438
|
1.40
|
7966.13
|
17.78
|
2
|
Hollis
|
64 472 200
|
5.38
|
190.55
|
0.33
|
3
|
Cumana
|
65 560 414
|
4.84
|
326.40
|
0.46
|
4
|
Caroni Swamp
|
78 559 144
|
2.02
|
341.48
|
10.48
|
5
|
Mamoral
|
72 054 356
|
28 638.50
|
243.19
|
0.00
|
6
|
Mamon Village
|
65 353 074
|
22.83
|
811.55
|
1.85
|
7
|
Ecclesville
|
66 759 470
|
101.53
|
575.91
|
0.67
|
8
|
Claxton Bay
|
79 587 104
|
203.26
|
556.64
|
0.87
|
9
|
Rousillac
|
76 745 030
|
3.06
|
404.69
|
0.42
|
10
|
Catshill
|
65 843 962
|
0.00
|
2301.52
|
3.72
|
11
|
Morne Diablo
|
57 884 090
|
7.50
|
570.90
|
18.64
|
12
|
Quinam
|
77 002 942
|
0.90
|
250.12
|
0.00
|
Viral sequences. Mamoral (28,638.50 RPM) overwhelmingly had the highest proportion of mapped viral reads followed by Claxton Bay (203.26 RPM) and Ecclesville (101.53 RPM) (Table 1). The other sites had relatively low proportions of viral reads (0.90 to 22.83 RMP) except for Catshill, where no viral read was detected. Generally, the majority of the viral RNA reads mapped to known insect specific viruses which included the Phasi Charoen-like and Humaita-Tubiacanga viruses. Additionally, RNA sequencing reads mapped to several viruses associated with a range of insects and feeding hosts and included some of human health importance. Trinity v2.9.1 assembly of mapped viral reads resulted in a total of 169 viral contigs and 5 unassembled reads (Table 2).
Table 2
Identification of viral contigs from Haemagogus mosquitoes collected in Trinidad.
Classification
|
Virus of nearest match
(% nucleotide identity)
|
Main reported host
|
Contig size range (nt)
|
Contig identity
|
# of contigs
|
Mosquito borne viruses
|
|
Family; Togaviridae Alphavirus ssRNA(+)
|
Mayaro Virus (99.29–100%)
|
Haemagogus spp
|
149–277
|
polyprotein gene, nsp1-nsp4
|
5a
|
Known insect specific viruses
|
|
Genus; Phlebovirus ssRNA(−)
|
Phasi Charoen-like virus
(94–100%)
|
Aedes aegypti
|
243–6675
|
glycoprotein, nucelocapsid gene
|
84
|
Unclassified ssRNA(+)
|
Humaita-Tubiacanga virus
(97–99%)
|
Aedes aegypti
|
243–1998
|
capsid, RNA dependant RNA polymerase, segment 1
|
14
|
Mosquito associated viruses
|
|
Family; Bunyaviridae Phasmavirus like. ssRNA(−)
|
Wuhan mosquito virus 6
(95–97%)
|
Culex pools
|
2162
|
complete cds
|
1
|
Unclassified
|
Kaiowa virus (99–100%)
|
Culex & Stegomiya spp
|
826–1172
|
Glycoprotein
|
2
|
Unclassified
|
Hubei virga-like virus 2
(94–95%)
|
Culicidae spp
|
547–5596
|
hypothetical protein, RdRp & putative coat protein
|
3
|
Family Rhabdoviridae ssRNA(−)
|
Merida virus
(91–95%)
|
Culex & Ochlerotatus spp
|
234–671
|
nucleoprotein, glycoprotein,phoshphoprotein,RNA dependent RNA polymerase gene
|
10
|
Unclassified
|
Cumbaru virus (93%)
|
Culicidae spp
|
203
|
Glycoprotein
|
1
|
Unclassified
|
Guato virus (98–99%)
|
Culex spp
|
256–818
|
hypothetical glycoprotein
|
5
|
Sobemovirus related ssRNA (+)
|
Guadeloupe mosquito virus strain (99.09 %)
|
Aedes aeqypti
|
2978
|
segment 1
|
1
|
Unclassified Picornavirales ssRNA (+)
|
Ecclesville picorna-like virus
|
No previous report
|
320–4519
|
polyprotein gene
|
5
|
Table 2
Classification
|
Virus or nearest match
|
Main reported host
|
Contig size range (nt)
|
Contig identity
|
# of contigs
|
Other insect associated viruses
|
|
Genus: Cripavirus ssRNA(+)
|
Black Queen cell virus
(96.02–97.14)
|
Apis spp
|
208–670
|
capsid protein, nonstructural & polyprotein gene
|
5
|
Other viruses
|
|
Genus: Alpharetrovirus ssRNA(RT)
|
Avian leukosis virus
(99–100%)
|
Birds
|
217–1020
|
envelope polyprotein, receptor binding protein, transmembrane envelope protein
|
33
|
Unclassified
|
Bat sobemovirus
(93.79–95.36%)
|
Bats
|
289–321
|
capsid protein
|
2
|
Unclassified
|
Bat luteovirus
(95.46– 96.37%)
|
Bats
|
441–442
|
capsid protein
|
2
|
Plant viruses
|
|
Genus:Tobamovirus ssRNA(+)
|
Pepper mild mottle virus
(98.68%)
|
Capsicum spp.
|
269
|
coat protein, replicase
|
1
|
Assembly was performed on mapped sequencing reads from the twelve locations and identified using the blastn suite from the NCBI nucleotide collection nt/nr database. The complete data set, which includes range of percentage nucleic acid identities and contig information, is shown in Table S1. aIncludes 3 paired reads from Ecclesville that assembled into 3 contigs and 2 unassembled paired reads from Mamon village. |
Known Insect Specific Viruses (ISVs). The metaviromes (Fig. S1) were dominated by one ISV, the Phasi Charoen-like virus (PCLV), which have been previously reported in Aedes aegypti (Ae. aegypti) mosquitoes 7,27,28. Overall, PCLV was identified in nine locations (1.65 to 17 064.28 RPM, mean = 1442.02 RPM) while the Humaiata-Tubiacanga virus (HTV), another well-known ISV, was identified in eight (0.09 to 710.30 RPM, mean = 59.33 RPM) of the twelve locations sampled (Fig. S2). De novo assembly of the viral reads using Trinity v2.9.1 resulted in 3 to 27 contigs (243 to 6675 bp) of PCLV per site with similarity levels to reference sequences in GenBank ranging from 95.44–99.87%. Phylogenetic analysis of segments of the S (551 bp), L (451bp) and M (439 bp) regions of representative sequences from the different sites showed most of the local sequences clustering into one main branch for all the segments with high bootstrap support (Fig. 2). However, the sequences for all three segments from the Caroni Bird Sanctuary clustered separately from other Trinidad sequences into well-supported clades (81–99% bootstrap support) including reference sequences (Fig. 2a – c), with segments S and M having closest links to sequences from Grenada. Additionally, the local sequences of segment L from Morne Diablo also separated and formed a separate clade (Fig. 2c; 100% bootstrap support).
Assembly of the HTV viral reads showed 1–3 contigs per site (243 to 1998 bp) with similarities to reference sequences in GenBank ranging from 97.80-99.13%. Phylogenetic analysis of a 188 bp region of segment 1 from the different sites showed all local sequences clustering separately from four references from Guadeloupe, with five (5) of the Trinidad sequences in one distinct sub-clade with 91% bootstrap support, the Caroni Bird Sanctuary and Mamoral sequences in a second strongly supported sub-clade (97% bootstrap support), and the Morne Diablo sequence on a weakly-supported branch (Fig. 2d).
Mosquito-borne viruses. A relatively small proportion of viral reads (4.80 RPM) were classified as alphaviruses (Togaviridae family), which accounted for approximately 0.0005% of all viral reads from the twelve locations (Fig. 1). Among these, 10 reads (4 from Mamon Village and 6 from Ecclesville at 0.06 and 0.10 RPM, respectively) mapped onto the genome of the MAYV. Trinity assembly of the mapped MAYV viral reads showed no overlapping among the read pairs. Blastn analysis revealed, all the reads had 99.29 %-100% similarity to MAYV non-structural (nsp1-3 and nsp4) and structural polyprotein coding regions (Table S1). Phylogenetic analysis revealed that all the MAYV reads from both locations diverged from the local reference 1954 Trinidad MAYV sequence and clustered with reference South American MAYV sequences from Venezuela, Bolivia, and Peru, with good bootstrap support (88%) for clades including MAYV from Peru and Venezuela in the case of Ecclesville read 1 and Mamon Village read 2 (Fig. 3).
Mosquito-associated viruses. A relatively small number of reads mapped to reference genomes of mosquito-associated viruses (Fig. 1) including the Wuhan mosquito virus 6 (Bunyaviridae) (35.13 RPM), Hubei virga-like virus 2 (Rhabdoviridae) (8.84 RPM) and the Merida virus (Unclassified) (1.86 RPM). These mosquitoes-associated viruses were identified in the Ecclesville, Chaguaramas, Claxton Bay and Mamoral locations. Overall, the Ecclesville metavirome (Fig. 1) was affiliated with a higher number of mosquito-associated viral reads (90.31 RPM) and showed the greatest diversity when compared to other sampling locations including Claxton Bay (1.92 RPM) and Chaguaramas (1.40 RPM).
All mosquito-associated viral reads identified in the Chaguaramas (Fig. 1a) and the Claxton Bay metavirome (Fig. S2) belonged to the unclassified viruses. The Guadeloupe mosquito virus sequence (Sobemovirus) was the only mosquito-associated sequence identified in the Mamoral location. Assembly of the mapped viral reads using Trinity v2.9.1 resulted in 1 contig (2162 bp) for the Wuhan mosquito 6 virus, 3 contigs (547–5596 bp) for Hubei virga-like virus, and 10 contigs (234–671 bp) for the Merida virus. Blastn analysis of the contigs of the three viruses showed highest identity of 96.85%, 99.14–99.63% and 97.17–99.15%, respectively, to GenBank reference sequences. Based on phylogenetic analysis for the segment 3 region for the Wuhan mosquito 6 virus (2133 bp), the local Ecclesville sequence diverged from the reference sequences from Australia (Fig. S3a). Furthermore, 1 contig (2978 bp) for the Guadeloupe mosquito virus was identified only in the Mamoral metavirome. Phylogenetic analysis of an 1859 bp region of segment 1 (Fig. S3b) showed the sequence separating from all reference sequences from Guadeloupe. However, it should be noted that bootstrap values were not significant for separation of the Trinidad viral sequences from reference sequences for Wuhan mosquito 6 and the Guadeloupe mosquito virus.
Other viruses. RNA sequencing reads from Eccleville also mapped to the picorna-like virus (44.40 RPM) and reads from the Chaguaramas site mapped to unclassified viruses including Kaiowa virus (1.16 RPM), Cumabru-like virus (0.08 RPM) and Guato-like virus (0.16 RPM). The Kaiowa virus (1.30 RPM) and Guato-like virus (0.62 RPM) were also detected in mosquitoes from Claxton Bay. Additionally, reads mapped to other viruses including the avian leukosis virus in Ecclesville (1.05 RPM), Mamon Village (2.57 RPM), Claxton Bay (0.58 RPM), Cumana (1.40 RP), Hollis (1.43 RPM), Morne Diablo (1.57 RPM) and Rousillac (0.69 RPM); the black queen cell virus in Quinam (0.90 RPM); and the pepper mild mottle virus in Claxton Bay (0.06 RPM).
Trinity v2.9.1 assembly of the reads resulted in 2 contigs (826–1172 bp) of the unclassified Kaiowa virus (99 to 100% similarity levels) in the Claxton Bay and Chaguaramas sites, as well as 1 contig (203 bp) of the Cumbaru-like virus (similarity level = 93%) and 4 contigs (256–300 bp) of the Guato-like viruses (similarity level = 98–99%) from the Chaguaramas location (Table S1). Kaiowa and Guato viruses were previously found in mosquitoes from Brazil 29.
Furthermore, 5 contigs (320–4519 nt) from Ecclesville were unique, with highest similarity levels (94.94 % to 96.94 %) to the Atrato picorna-like virus. Phylogenetic analysis of a 350 bp region of the polyprotein segment (Fig. S3c) resulted in all Ecclesville contigs forming a single clade that diverged from the GenBank reference Atrato picorna-like virus sequence. A relatively large number (33) of contigs (217 to 1020 nt) from 7 locations (Hollis, Cumana, Mamon Village, Ecclesville, Claxton Bay, Rousillac and Morne Diablo) aligned with high nucleotide identity (99–100%) to the avian leukosis virus strains from China. Maximum likelihood phylogeny of the avian leukosis viral contigs resulted in local sequences separating into two main clades with sequences from Cumana and Hollis clustering together (Fig. S3d) and sequences from Rousillac and Mamon village forming a strongly supported clade. This suggests the local avian leukosis viruses are highly diverse and may include unique lineages as compared to the reference sequences available in the GenBank database. Four contigs (289-442bp) from the Ecclesville location aligned with 94–96% nucleotide identity to unclassified bat sobemoviruses and luteoviruses (Table S1).
Host sequences
Mapping of RNA reads to potential mosquito bloodmeal host genomes showed that, although levels were relatively low, the most hits in all twelve locations were against birds. Several reads mapped to the Gallus gallus representative genome (0.61 RPM for Claxton Bay, 1.85 RPM for Mamon Village, 0.49 RPM for Rousillac and 0.75 RPM for Hollis). Assembled contigs aligned with 100 % nucleotide identity to Gallus gallus (Table S3 and S4). Additionally, some reads also mapped to rodent genomes (17.73RPM for Chaguaramas and 3.72 RPM for Catshill).
Bacterial sequences
In all twelve locations, a relatively moderate level of reads (717.28 RPM) mapped to the SILVA bacterial 16S rRNA gene sequence database (Version 132). Most reads were classified within the Proteobacteria phylum and specifically within the Gammaproteobacteria class (Fig.4). The bacterial profile generally showed some similarity throughout all sampling locations (Fig.4), with a core microbiome comprising of several genera including Pseudomonas, Acinetobacter, Enterobacter and Serratia genera. Furthermore, the microbiome from the Claxton Bay (61.13 RPM), Rousillac (269.48 RPM), Chaguaramas (3777.92 RPM) and Ecclesville (188.45 RPM) locations all showed reads that mapped against the Wolbachia genus. Blastn analysis of assembled 16S contigs (234 -1569 nt) from these four sites showed 98.79 – 100 % similarity to Wolbachia species.