Abundance and distribution of sylvatic dengue virus vectors in three different land cover types in Sarawak, Malaysian Borneo

Study area

The study was conducted in and around Kampung Puruh Karu (1°15′19.09″N, 110°16′58.61″E), approximately 43 km southeast of Kuching city in Sarawak, Malaysian Borneo (Fig. 1). The climate of Borneo is typified by warm temperatures, averaging between 25 and 30 °C, high relative humidity (typically over 70%) and rainfall throughout the year (Fig. 2). The study area is predominately composed of rural communities surrounded by agricultural fields and forest patches. Communities comprise homesteads where extended families live together in separate households and maintain small gardens including fruiting trees, small rice paddies, and other agricultural crops. Most homesteads also own separate agricultural plots away from the community where they conduct large scale farming of either sustenance or economic crops. These latter plots have mainly been created by slash and burn methods within forests close to the community. According to local people, most of the forested area surrounding Kampung Puruh Karu had been cleared within the last twenty to thirty years, and these forests are now dominated by several different bamboo species, dipterocarp trees, fruiting trees, and tall grasses.

Climate data and land cover delineation

Data from the Tropical Rainfall Measuring Mission (TRMM) 3B42 V7 product and the level-3 Moderate Resolution Imaging Spectroradiometer (MODIS) Land Surface Temperature (LST) and Emissivity A2 product were used to compare annual precipitation and temperature patterns for 2013 and 2014. Daily-accumulated rainfall measures from TRMM for a central location within the study area (1°17′03.00″N, 110°16′52.00″E) were used to generate estimates of monthly averages of daily rainfall for the study area. The MODIS LST data were retrieved as averaged values of daytime clear-sky LST for 8-day periods at a 1-km resolution covering the study area. Data from 8-day periods with >50% missing pixels due to cloud cover were excluded. The 8-day average temperatures were then used to generate approximate monthly average temperatures for the study area. To compare the patterns observed for these two years to the normal trends in the area, we included data on average temperature and precipitation for the study area from WorldClim interpolated climate layers, which show averages from 1950 to 2000 at 1 km spatial resolution.

The land cover map was created through on-screen digitizing in ArcMap 10.2, using Google Earth imagery from April 22nd, 2014. The three land cover types were assessed visually and delineated. Areas that appeared to be either bare earth or non-forested vegetation were classified as “Cleared.”

Mosquito collection

Mosquito sampling was conducted using both gravid traps (Bioquip, Rancho Dominguez, California, USA) and a backpack aspirator (Bioquip). The order of sampling each site was haphazard as sampling could only commence after permission was granted from landowners. Collections were taken for 3 consecutive days at a given sampling site. Gravid traps were run continuously over these 3 days and samples were collected daily from 9 gravid traps per site at 08:00 and 16:00 h. Temperature at the time of collection was recorded for each trap. Backpack collections were performed twice daily at 08:00 and 16:00 h at each of the gravid trap locations. The backpack was run for 5 min while standing next to the trap location and moving the collection receptacle through the air, across vegetation if near the trap location and around the body of the collector.

Nine traps of each type were placed at each site; one trap at the centroid of the site, and two traps along four 90 degree transects: one half way between the centroid and the perimeter and one at the intersection of each transect and the perimeter. Mosquitoes were removed from each trap and anaesthetized using chloroform. Individual mosquitoes were identified morphologically to genus and, when possible, species, using taxonomic keys [17–19] and photographed. Samples were then pooled (< 40 mosquitoes/pool) by location, trap type, date collected, genus, known or putative species group (i.e. Aedes sp. A for a putative Aedes species), and sex. The pooled samples were held in liquid nitrogen until transported to the laboratory at the Institute of Health and Community Medicine at the Universiti Malaysia Sarawak (UNIMAS) where they were stored at -80 °C. All sampling locations were briefly resampled in April and May of 2016 using the backpack aspirator methods described above; however, each location was only sampled for 1 day rather than 3 consecutive days. Mosquito identification during resampling was based solely on morphological examination using available taxonomic keys; as well as the WRBU online interactive dichotomous key and unpublished reference guides made by the Sarawak Vector Control Department [7, 17–19].

Virus screening

At UNIMAS, mosquito pools were homogenized by bead beating for 3 min at 26 rounds per min (rpm) in 1 ml of prepared stock homogenization media (500 ml DMEM hi-glucose (Gibco, Pittsburgh, Pennsylvania, USA) supplemented with 2.5 μg/ml amphotericin B (Gibco) 10% fetal bovine serum (FBS) (Gibco) and 1× penicillin-streptomycin mixture (Gibco). After homogenization, samples were centrifuged for 5 min at 8000 rpm, and 400 μl of supernatant was passed through a 0.25 μm syringe filter and stored at -80 °C until screened for arboviruses. Samples were then shipped on dry ice to University of Texas Medical Branch (UTMB) for arbovirus assay, after which the remaining homogenate was returned to -80 °C and shipped to New Mexico State University (NMSU) for DNA barcoding.

At UTMB, 100 μl of each filtered mosquito homogenate was inoculated into 12.5 cm² culture flasks of African green monkey kidney (Vero) or Ae. albopictus (C6/36) cells. Vero cell cultures were held for 14 days at 37 °C and examined regularly for evidence of viral cytopathic effect (CPE). The C6/36 cell cultures were incubated at 28 °C for 7 days and observed daily for CPE. If CPE was observed, the cell culture supernatant was collected from the flask, clarified by centrifugation at 3200 rpm for 5 min, and transferred to a clean 1.5 ml Eppendorf tube. Total RNA was then extracted from an aliquot of the clarified cell culture supernatant through the Trizol method and resuspended in 50 μl RNase/DNase and protease-free water (Ambion, Austin, Texas). Additional aliquots were stored at -80 °C to attempt isolation of viruses from samples revealed to contain viral genomes by next-generation sequencing. Sequencing for viral screening was conducted on 14 CPE positive Aedes pools total.

Next-generation sequencing to detect viruses

RNA (0.1–0.2 μg) quantified by NanoDrop 1000 spectrophotometer (Thermo Fisher Scientific, Pittsburgh, Pennsylvania, USA), was fragmented by incubation at 94 °C for 8 min in 19.5 μl of fragmentation buffer (Illumina, San Diego, California, USA). Following fragmentation, first and second strand synthesis, adapter ligation and amplification of the library were performed using the TruSeq RNA Library Prep Kit v2 under conditions prescribed by the manufacturer (Illumina). The samples were tracked using the index tags incorporated into the adapters as defined by the manufacturer.

Cluster formation of the library DNA templates was performed with the Illumina TruSeq PE Cluster Kit v3 and the Illumina cBot workstation using conditions recommended by the manufacturer. Paired end 50 base sequencing by synthesis was performed using Illumina TruSeq SBS kit v3 on an Illumina HiSeq 1000 using protocols defined by the manufacturer. Cluster density per lane was 645–980 k/mm² and post-filter reads ranged from 148 to 178 million per lane. Base call conversion to sequence reads was performed using CASAVA-1.8.2. The de novo assembly program ABySS [20] was used to assemble the reads into contigs, using several different sets of reads, and k values from 20 to 40. In certain cases, pre-filtering of host-derived reads by mapping to Ae. albopictus and Ae. aegypti reference sequences enhanced the assembly process. Reads were mapped back to the contigs using bowtie2 [21] and visualized with the Integrated Genomics Viewer [22] to verify that the assembled contigs were correct.

Mosquito DNA barcoding

Mosquito homogenates were subject to mitochondrial cytochrome c oxidase subunit 1 (cox1) barcoding at NMSU. The steel bead was removed from individual pools by magnet and each sample was centrifuged at 13,000× g for 5 min at 4 °C. The supernatant was removed and discarded; total DNA and RNA were extracted from the remaining homogenates using Trizol reagent (Invitrogen, Pittsburgh, Pennsylvania) according to the manufacturer’s protocol. The resulting DNA and RNA pellets were air-dried for up to 1 h and reconstituted in 10 μl of DEPC treated nuclease free water (Invitrogen).

The cox1 gene region has been used extensively to identify Aedes and Culex mosquito species and to perform evolutionary analyses. The universal forward primer LCO1490 (5′-GGT CAA CAA ATC ATA AAG ATA TTG G-3′) and reverse primer HC02198 (5′-TAA ACT TCA GGG TGA CCA AAA AAT CA-3′) described by Folmer et al. [23] were used to amplify a 710 base-pair amplicon of cox1 for all Aedes species collected. PCR was performed on samples with greater than 20 ng/μl of DNA using the New England biolabs Quick-Load Taq 2× Master Mix (New England Biolabs, Ipswitch, Massachusetts, USA) following the manufacturer’s instructions. The following PCR conditions were used for amplification: an initial denaturing cycle for 2 min at 95 °C, 30 cycles of melting 95 °C for 25 s, annealing at 57 °C for 37 s (with annealing temperature decreasing by 0.4 °C per cycle), extension at 72 °C for 1 min 30 s, followed by an additional 4 cycles with denaturing at 95 °C for 25 s, annealing at 45 °C for 30 s and decreasing by 0.4 °C per cycle, extension at 72 °C for 1 min and 30 s and finished with a final extension step at 72 °C for 10 min. For samples with DNA concentrations beneath the cut-off, RT-PCR was performed on isolated RNA. All RT-PCR reactions were performed by following the recommended protocol of the Roche Titan one-tube RT-PCR system (Roche, Indianapolis, Indiana, USA) with the primers described above. The following conditions were used for all RT-PCR reactions: 1 RT cycle at 45 °C for 45 min followed by 1 denaturing cycle at 94 °C for 4 min, annealing at 45 °C for 1 min, extension at 68 °C for 1 min, followed by 3 cycles of denaturing at 94 °C for 20 s, annealing at 45 °C for 1 min, extension at 68 °C for 1 min, followed by 10 cycles of denaturing at 94 °C for 20 s, annealing at 50 °C for 30 s, extension at 68 °C for 1 min, followed by 16 cycles of denaturing at 94 °C for 20 s, annealing at 50 °C for 30 s, extension at 68 °C for 1 min and increasing by 20 s per cycle and finishing with one extension at 68 °C for 5 min.

All PCR and RT-PCR products were visualized on a 1.5% agarose gel and amplified products were purified using the High Pure PCR product purification kit (Roche) following the manufacturer’s recommended protocol. Following purification, samples were sent for sequencing using both forward and reverse primers (Eurofins genomics, Huntsville, Alabama).

Sequence alignment and phylogeny reconstruction

Forward and reverse cox1 sequences were reviewed for quality and if both reads were successful then these were combined in a single contig; if only one of the two was available then this read was used. The sequences were trimmed and compared using BLAST against the NCBI database. BLAST results with the highest max score were compared to pool sequences and results with a minimum sequence identity ≥ 98% were considered conspecific and sequences with sequence identity > 80% were considered congeners, following the cutoffs laid out by Wang et al. [24].

To review the relatedness of BLAST results to pool sequences, a phylogeny was reconstructed from an alignment of select Aedes pool sequences and reference sequences from a wide geographic range retrieved from GenBank. Representative sequences were used for sequences that had 100% sequence identity to each other in order to simplify the phylogeny. All representative sequences and references were aligned using ClustalW and a maximum likelihood phylogeny was created using the Tamura-Nei parameter model and 500 bootstrap replicates. All sequence and phylogenetic analyses were conducted using Geneious version 9.1.4 [25].

Statistical analyses

Mosquito abundance in different land covers (n = 5 sites/land cover) was tested for normality and then compared using a one-way ANOVA followed by a Tukey-Kramer post-hoc test to reveal specific differences among land cover types. The proportion of a particular genus or species in two different time periods or two different land covers or groups of land covers was compared using a Fisher’s exact test. Finally, the number of sites that were positive or negative for viral detection in mosquito pools was compared with a Fisher’s exact test, using exceptions described by Freeman & Halton [26].

Climate during the 2013–2014 sampling period

Although temperature remained relatively consistent between the two sampling periods (Fig. 2), the monthly average precipitation was higher during the 2013 sampling effort compared to 2014, as expected due to seasonal changes in rainfall.

Mosquito diversity in the 2013 and 2014 collections

A total of 2164 mosquitoes were collected over the two field seasons, 1201 in 2013 and 963 in 2014, from 5 sites each in homestead, agriculture and forest land covers. Fifty-seven percent were female. All 2164 mosquitoes were pooled into 328 total pools. Based on morphological and molecular identification, the mosquitoes collected belong to 10 genera (Fig. 3). The relative frequencies of genera collected remained consistent over the two years with Aedes comprising 79.5% of mosquitoes collected in 2013 and 78.6% in 2014. Culex was the second most common taxonomic group, constituting 7.6% of specimens in 2013 and 16.1% in 2014 (Fig. 3). The remaining genera, Anopheles, Armigeres, Coquillettidia, Lutzia, Mansonia, Toxorhynchites, Uranotaenia and Zeugnomyia were collected at low frequencies, ≤ 4%, over the two sampling years. Mosquitoes that were unidentifiable, or unknown, accounted for 4.1 and 3.1% of the collection, respectively, in the two years. The combined frequencies of the genera collected for 2013 and 2014 are represented in Fig. 3.

Land cover change and mosquito diversity in the 2016 collection

Of the 15 sites sampled in 2013 and 2014, 2 forests and 1 farm had been altered in 2016. The forests were cleared using the slash and burn method and the farm was converted from solely black pepper to mixed agriculture, demonstrating the quickly shifting landscape in this area. A total of 283 mosquitoes were collected from all sites in 2016, which now included 3 forest sites, 5 agricultural sites, 5 homestead sites, and 1 site each in two new land cover categories: highly degraded forest and barren. Of this collection, 93.2% were Aedes and the remaining 6.7% was comprised of Armigeres, Coquillettidia, Culex and Uranotaenia. To test whether alteration to sampling sites effected the proportion of Aedes sampled in 2016, we compared the proportion of Aedes in the total collection of mosquitoes, sampled in 2013 and 2014 to the proportion sampled in 2016. The proportion of Aedes collected in the 2016 resampling was significantly greater than in previous sampled years, 2013 and 2014 combined, before the alteration of the sampling sites (N = 2447, P = < 0.0001). This difference persisted even after excluding the altered sites from the 2016 resampling (N = 2437, P < 0.0001). We then compared the proportion of Aedes sampled in 2013 and 2014 to the proportion sampled in 2016 at each of the 3 altered sites. A total of 179 Aedes mosquitoes were sampled at forest site, FSTAT in 2013 and 2014, prior to its conversion to barren land; however, in 2016, no mosquitoes were collected at this site. The proportion of Aedes sampled at forest site FSTST (39% of 74), was significantly lower than from the same site after conversion to highly degraded forest in 2016 (86% of 7) (N = 81, P = 0.04. Finally, the proportion of Aedes sampled at agriculture site, FRMCO was higher before (100% of 80) than after (67% of 3) conversion from black pepper monoculture to mixed agriculture, however the extremely small sample size of the post-conversion collection precluded statistical analysis.

Mosquito DNA barcoding

While studies of the ecology of sylvatic mosquitoes in southeast Asia have generally relied exclusively on morphological identification, e.g. [27–30], we felt that in a biodiversity hotspot like Borneo it would also be prudent to barcode specimens. DNA or RNA was extracted from 138 of 193 pools morphologically identified as Aedes and subjected to DNA barcoding. A total of 122 of these pools had been identified as Ae. albopictus based on morphology; given the high morphological similarity in these individuals, we limited barcoding to 44 pools, which together included 498 mosquitoes. All of these samples, which covered all three land cover types, homestead, forest and agriculture, had between 98 and 99% sequence identity to Ae. albopictus sequences when BLAST against the NCBI database (Additional file 1: Table S1). When sequences from all 44 Ae. albopictus pools were aligned, the average percent sequence identity for all pairwise comparisons was 99.9%, indicating a high degree of genetic similarity between Ae. albopictus populations in the study area. When representative Ae. albopictus pool sequences were phylogenetically compared to reference sequences from GenBank, all of the morphologically identified Ae. albopictus pools migrated together and showed a close relationship to reference Ae. albopictus sequences (Fig. 4).

Of 8 pools of Ae. niveus collected, DNA or RNA was successfully amplified from 1 pool, which contained a single mosquito. This sequence did not return a match when BLAST against the NCBI database; an extensive search revealed that no cox1 sequences of Ae. niveus have been uploaded to GenBank. Aedes niveus is a highly distinctive group and we have high confidence in our morphological identification; moreover, it has previously been reported in Sarawak by Macdonald et al. in 1965 [28]. One cox1 sequence of an unconfirmed Ae. niveoides isolate from China is available in the NCBI database. Aedes niveoides is a species in the niveus complex that was also identified in Sarawak by MacDonald et al. [28]. A pairwise comparison of this sequence with the Ae. niveus sequence indicates 90.2% identity, indicating the two species to be congeners. Our phylogenetic comparison of representative Aedes sequences and NCBI references also shows this close relationship between the Ae. niveus collected in Sarawak and the Ae. niveoides isolate from China (Fig. 4). Thus we here report the cox1 sequence of Ae. niveus spp. for the first time.

A total of 275 morphologically identical Aedes that we were unable to identify to species levels were designated Aedes sp. A. The average percent sequence identity among cox1 sequences from 18 pools of Aedes sp. A, which together included 227 mosquitoes, was 100%, confirming that they are a single species. A BLAST search with these sequences returned a result of Ae. cogilli with between 89 and 91% sequence identity (Additional file 1: Table S1). Currently, Ae. cogilli has only been described on the Indian subcontinent; however, it is related to the Aedes niveus species complex [31]. A phylogeny of representative Aedes sp. A and reference Ae. cogilli sequences further supports this description, as the collected mosquito pools are more closely related to other Aedes sp. A sampled in Sarawak than reference sequences from Pakistan and India (Fig. 4). However, in light of new mosquito collections under way in Sarawak from 2016 and the advancement of available keys from Sarawak, Aedes sp. A has now been morphologically identified as Ae. desmotes. The average percent sequence identity between an unconfirmed Ae. desmotes reference sequence retrieved from GenBank and the Aedes sp. A sequences was 94.9% which ranks these specimens as only congeners. However, the reconstructed phylogeny supports this morphological identification over Ae. cogilli since the representative Aedes sp. A sequences migrate with the reference Ae. desmotes sequence from NCBI. Sequencing of other regions of the genome is needed to adequately identify Aedes sp. A. A second group of 6 morphologically identical mosquitoes that we were unable to identify were designated Aedes sp. B. The cox1 sequence from one pool of this species, which included 2 mosquitoes in total, returned a BLAST result of Aedes japonicus with 88% sequence identity (Additional file 1: Table S1). Thus this species is closely related to Ae. japonicus, a medically important vector of WNV, albeit not a species that has been documented in Borneo [7, 32]. Phylogenetically, this species is more closely related to the Ae. niveus specimen collected in Sarawak than to reference Ae. japonicus sequences (Fig. 4).

Two groups of morphologically similar Aedes species remain unidentified and have been designated as Aedes sp. C, which includes 14 mosquitoes in total, and Aedes sp. D, which includes 3 mosquitoes in total. Unfortunately, a PCR product was not obtained for either species group using the primers described and therefore a molecular identification for these species could not be obtained. Although these species groups were not identified via barcoding, they were nonetheless included in the rarefaction analysis described below.

Mosquito abundance and distribution in the 2013 and 2014 collection

There was no significant difference in the total number of mosquitoes collected among the three land cover types (Fig. 5) (F _(2,12) = 0.75, P = 0.49). Since it was dry prior to sampling in 2014 (Fig. 2), we also analyzed the data from 2013 only (3 sites per land cover); in this comparison also there was no significant difference in the total number of mosquitoes collected among land cover types (Fig. 5) (F _(2,6) = 0.23, P = 0.80). Although there was more limited sampling in 2016, in line with our previous findings, there was no significant difference in the total number of mosquitoes collected among land cover types (F _(2,12) = 0.23, P = 0.13).

To test whether our sampling efforts were intensive enough to explain the richness of Aedes species in all three land cover types, we performed a rarefaction analysis on the total number of Aedes collected in all homestead, agriculture and forest sites. This rarefaction analysis indicates that much of the Aedes species diversity in the forests remain to be discovered; however, our sampling of agricultural and homestead sites detected most of the Aedes species diversity in these land cover types for this area (Fig. 6).

Sampling in 2016 was too limited for statistical analysis.

Abundance and distribution of target Aedes vectors

The number of Ae. albopictus collected among the three land cover types in 2013 and 2014 (Fig. 7) did differ significantly (F _(2,12) = 4.08, P = 0.04); more Ae. albopictus were collected in agricultural sites than forests. When the analysis was limited to just 2013, this difference remained marginally significant (F _(2,6) = 4.96, P = 0.05) and Ae. albopictus showed a tendency to be more abundant in agricultural sites. We did not analyze the 2014 data by itself as only 2 sites per land cover were sampled in 2014. Counter to our previous findings in 2013 and 2014, there was no significant difference in the mean number of Ae. albopictus collected among the three land cover types in 2016 when the altered sites were excluded (mean abundance, HOM = 8.0 ± 8.2, FRM = 41.5 ± 9.1, FST = 15.0 ± 10.5; F _(2,9) = 3.95, P = 0.06). Additionally, there was no significant difference in the percent of Ae. albopictus collected among Aedes mosquitoes at forest site FSTST, before (45%) and after (33%) conversion to highly degraded forest (N = 35, P = 0.68).

Virus detection

Cell culture for arbovirus isolation was performed on 120 Aedes pools collected in 2013, including 70 Ae. albopictus, 8 Ae. niveus, 34 Aedes sp. A, 3 Aedes sp. B, 3 pools of Aedes sp. C and 2 pools of Aedes sp. D, and 73 Aedes pools collected in 2014, including 52 Ae. albopictus, 11 Aedes sp. A and 10 pools of Aedes sp. C. No virus isolations were made from Aedes mosquito pool homogenates placed into Vero cell culture; therefore, all data is based on genome detection by next generation sequencing (NGS) of C6/36 cell cultures showing CPE. NGS was performed on 14 Aedes pools that were CPE-positive in C6/36 cell cultures. Genomes of 14 different viruses were detected in 12 of these pools, which included 11 Aedes albopictus pools and 1 unidentified Aedes sp. C pool. Two pools of Ae. albopictus contained genomes from 2 viruses (Table 3). The majority of viruses detected were insect-only viruses including Euprosterna elaeasa virus, Aedes pseudoscutellaris reovirus, and the newly described Kampung Karu virus [33] (Table 3). Kampung Karu virus is a flavivirus that was isolated from a pool of Anopheles tessellatus during the 2013 sampling of this study and lies within a clade of insect only flaviviruses closely related to other flaviviruses that infect humans including DENV [33]. Viruses were detected from mosquito pools collected in all land cover types; however, the majority of identified viruses were either collected in homestead or agriculture land cover types. There was no difference among homestead, agriculture or forest land cover types in the proportion of sites in which a virus was detected (N = 15, P = 0.80). The small number of samples in which viruses were detected precluded statistical analysis; however, Fig. 8 shows the mean number of viruses detected per land cover type.

Mosquito species	Land class of collection	Sex (M/F) of mosquitoes in pool	No. of mosquitoes in pool	Viruses detected
Aedes albopictus	HOM	M	31	Euprosterna elaeasa virus; Aedes pseudoscutellaris reovirus
Aedes albopictus	HOM	F	11	Euprosterna elaeasa virus
Aedes albopictus	HOM	F	13	Euprosterna elaeasa virus
Aedes albopictus	HOM	M	27	Kampung Karu virus
Aedes albopictus	HOM	M	1	Euprosterna elaeasa virus
Aedes albopictus	FRM	M	32	Euprosterna elaeasa virus; Aedes pseudoscutellaris reovirus
Aedes albopictus	FRM	M	30	Euprosterna elaeasa virus; Aedes pseudoscutellaris reovirus
Aedes albopictus	FRM	M	22	Euprosterna elaeasa virus
Aedes albopictus	HOM	M	8	Euprosterna elaeasa virus
Aedes albopictus	FRM	F	22	Euprosterna elaeasa virus
Aedes unknown	FST	F	1	Euprosterna elaeasa virus
Aedes albopictus	FRM	F	7	Euprosterna elaeasa virus

Abbreviations: F, female; M, male HOM, homestead FRM, agriculture FST, forest