Introduction
Factors associated with global change, such as temperature increase, land-use change and the increasing spread of invasive species, are leading to a considerable loss and reorganization of biodiversity (Harold, Reference Harold and Richard2000; Segan et al., Reference Segan, Murray and Watson2016; Eriksson and Hillebrand, Reference Eriksson and Hillebrand2019), with important consequences for the emergence of infectious diseases that affect wildlife, livestock and human populations (Caminade et al., Reference Caminade, McIntyre and Jones2019; Smith et al., Reference Smith, Machalaba, Seifman, Feferholtz and Karesh2019). Among global emerging infectious disease events, vector-borne diseases are disproportionately over-represented (Swei et al., Reference Swei, Couper, Coffey, Kapan and Bennett2020) and constant efforts for monitoring insect vector populations should be carried out in locations at risk (Pedersen et al., Reference Pedersen, Stolk, Laney and Michael2009). Among insect vectors, mosquitoes (Culicidae), with more than 3500 described species widely distributed around the world (Harbach, Reference Harbach2013), are considered among the main insect vectors involved in the transmission of pathogens including viruses, protozoans and nematodes. Three main genera, Anopheles, Aedes and Culex are considered of medical importance for humans and transmit pathogens causing diseases to more than 700 million people annually, resulting in over 1 million deaths (World Health Organization, 2020). In the last few decades, the rapid worldwide spread of the invasive yellow fever mosquito Aedes aegypti and the Asian tiger mosquito Aedes albopictus is producing novel epidemiological scenarios (Bonizzoni et al., Reference Bonizzoni, Gasperi, Chen and James2013; Iwamura et al., Reference Iwamura, Guzman-Holst and Murray2020). Early detection of mosquito invasion events, as well as continued surveillance of such invasions, is becoming essential to assess the risks of local mosquito-borne disease outbreaks. In addition, it seems essential to understand the ecological interactions between mosquito species at breeding sites to evaluate the competitiveness of indigenous species (Juliano, Reference Juliano2009).
To date, surveys of mosquito species have been performed with traditional sampling using ovitraps and dipping method for immature stages or with light/decoy traps and human landing catches for adults (Focks, Reference Focks2004). Skilled entomologists are able to identify specimens using morphological traits (Besansky et al., Reference Besansky, Severson and Ferdig2003; Hajibabaei et al., Reference Hajibabaei, Singer, Hebert and Hickey2007), however some species are indistinguishable morphologically (e.g., cryptic species of Anopheles) (Coetzee and Koekemoer, Reference Coetzee and Koekemoer2013). In addition, the identification of different mosquito stages (i.e., eggs, larvae and adult mosquito specimens) needs solid knowledge from experts in entomology. Identification can be time consuming, especially when there are many species such as in the tropics, but also characterization at the species level becomes very difficult, if not impossible, when specimens are too damaged (Foley et al., Reference Foley, Rueda and Wilkerson2007). Developments in molecular techniques over the past decade, coupled with reduced sequencing costs, have made use of environmental DNA (eDNA) as an approach with a huge potential to survey micro-biodiversity in the field. eDNA is DNA that is shed by organisms (e.g., through faecal waste, dead skin, gastrointestinal tract cells, gametes or via post-mortem degradation), and it has formed the basis of numerous studies focussed on vertebrate detection (Ficetola et al., Reference Ficetola, Miaud, Pompanon and Taberlet2008; Goldberg et al., Reference Goldberg, Pilliod, Arkle and Waits2011; Jerde et al., Reference Jerde, Mahon, Chadderton and Lodge2011; Minamoto et al., Reference Minamoto, Yamanaka, Takahara, Honjo and Kawabata2012; Thomsen et al., Reference Thomsen, Kielgast, Iversen, Wiuf, Rasmussen, Gilbert, Orlando and Willerslev2012; Spear et al., Reference Spear, Groves, Williams and Waits2015; Egeter et al., Reference Egeter, Peixoto, Brito, Jarman, Puppo and Velo-Antón2018), and more recently for invasive invertebrates (Clusa et al., Reference Clusa, Miralles, Basanta, Escot and García-Vázquez2017; Klymus et al., Reference Klymus, Marshall and Stepien2017; Mychek-Londer et al., Reference Mychek-Londer, Balasingham and Heath2020). In natural habitats, eDNA is affected by a variety of factors, such as temperature, microbial activity, pH (Seymour et al., Reference Seymour, Durance, Cosby, Ransom-Jones, Deiner and Ormerod2018), conductivity (Collins et al., Reference Collins, Wangensteen and O'Gorman2018), water chemistry or ultraviolet radiation. It is degraded over time, but can remain at detectable levels weeks after an organism's removal (Dejean et al., Reference Dejean, Valentini, Duparc, Pellier-Cuit, Pompanon, Taberlet and Miaud2011; Barnes et al., Reference Barnes, Turner, Jerde, Renshaw, Chadderton and Lodge2014; Pilliod et al., Reference Pilliod, Goldberg, Arkle and Waits2014). Hence, most eDNA detection is expected to indicate a current or recent colonization of the habitat (Piaggio et al., Reference Piaggio, Engeman, Hopken, Humphrey, Keacher, Bruce and Avery2014), making it a potentially suitable method for contemporary surveillance of aquatic populations, such as mosquito aquatic stages.
Although studies have shown the usefulness of eDNA metabarcoding for the monitoring of numerous invertebrate species, to the best of our knowledge, only few studies have demonstrated the usefulness of this technique for detection of mosquito species in particular. Schneider et al. (Reference Schneider, Valentini, Dejean, Montarsi, Taberlet, Glaizot and Fumagalli2016) analysed the potential of eDNA for the detection of invasive Aedes mosquitoes in Europe. They collected water samples in the field and used both quantitative real-time polymerase chain reaction (qPCR) and eDNA metabarcoding of a short fragment of the 16S rRNA gene of the Culicidae family. Both molecular methods gave comparable results and performed better than the traditional survey methods, however, the detection capacity decreased by half 10 days after the removal of the larvae. Those authors recommended for the eDNA approach to be used as a complement to traditional captures. Two other studies compared the effectiveness of eDNA approaches with traditional sampling techniques to detect mosquito larvae diversity in the field (Boerlijst et al., Reference Boerlijst, Trimbos, Van der Beek, Dijkstra, Van der Hoorn and Schrama2019; Krol et al., Reference Krol, Van der Hoorn, Gorsich, Trimbos, Bodegom and Schrama2019). These studies both used eDNA primers targeting the cytochrome oxidase I (COI) gene. Boerlijst et al. (Reference Boerlijst, Trimbos, Van der Beek, Dijkstra, Van der Hoorn and Schrama2019) found that 98% of the Culicidae species were correctly identified using eDNA, suggesting that eDNA-based approaches are reliable and can be even more reliable than traditional dipping methods for certain species. However, both studies yielded only a subset of the adult community known in their field sites. Species that were detected with eDNA were generally the most abundant species in the traps indicating that the eDNA metabarcoding method was more likely to pick up more abundant species than rare mosquito species (Krol et al., Reference Krol, Van der Hoorn, Gorsich, Trimbos, Bodegom and Schrama2019). Although eDNA metabarcoding can increase the accuracy of identification, while reducing the cost and time, compared to classical barcoding, it must be integrated with classical taxonomy and molecular methods for comprehensive ecological studies (Ruppert et al., Reference Ruppert, Kline and Rahman2019). The use of eDNA is a booming technique, but also has many limitations, including the degradation of eDNA in the environment, especially in tropical regions, as well as the sample preservation methods. In addition, one of the important considerations in eDNA metabarcoding studies is the primer design (Ruppert et al., Reference Ruppert, Kline and Rahman2019). Primers for different genes vary in coverage, resolution and inter-taxon bias. COI gene is the standard gene for barcoding, but other regions such as 12S or 16S ribosomal RNA may be more appropriate for different taxa (Epp et al., Reference Epp, Boessenkool, Bellemain, Haile, Esposito and Riaz2012; Taberlet et al., Reference Taberlet, Prud'Homme, Campione, Roy, Miquel and Shehzad2012; Deiner et al., Reference Deiner, Bik, Mächler, Seymour, Lacoursière-Roussel and Altermatt2017; Hering et al., Reference Hering, Borja, Jones, Pont, Boets and Bouchez2018). Primers for eDNA metabarcoding must be short enough to amplify degraded samples, identical for the same species, but variable between species, allowing amplification of a variety of species (Epp et al., Reference Epp, Boessenkool, Bellemain, Haile, Esposito and Riaz2012).
In our study on São Tomé Island, Gulf of Guinea (Africa), we wanted to evaluate the richness of mosquito species along a gradient of anthropogenic disturbances in order to compare the assemblage of species between human habitation areas (i.e., village with domestic animals), intensive agricultural areas (i.e., oil palm plantations) and natural neighbouring forested areas. To assess the mosquito richness at these three habitat types, we collected (i) water from larval breeding sites to perform eDNA metabarcoding using COI and (ii) adult specimens using CDC light traps set up in trees. The aims of this study were (i) to compare our metabarcoding results with the visual identification of larvae and the light traps captures, taking into account the samples’ characteristics (i.e., water turbidity, containers), (ii) to identify the assemblage of mosquito species along a gradient of anthropogenic disturbance, (iii) to detect the presence of the invasive tiger mosquito Ae. albopictus which recently colonized the island (Reis et al., Reference Reis, Cornel, Melo, Pereira and Loiseau2017) and finally (iv) to perform a short review of the pros and cons of the eDNA metabarcoding as a complementary methodological approach to traditional ones.
Materials and methods
Study sites and sampling
Water sampling took place in three different types of habitats in October 2019 on São Tomé Island (Gulf of Guinea, Africa): (i) a small village located in the middle of the oil palm plantation (0°6′57.308″ N; 6°35′33.414″ E), (ii) the oil palm plantation that surrounds the village and (iii) the secondary rainforest adjacent to the plantation at 1 km from the village (fig. 1).
Left: Map of São Tomé Island (Gulf of Guinea, Africa), with the black frame representing the sampling area in the southeast of the island. Right: Satellite picture of that area, with the village (circled in red), surrounding by the oil palm plantation; the green line being the border between the oil palm plantation and secondary forest.
We collected 37 water samples (30 ml each, with 10 ml of Longmire solution added for preservation) (Williams et al., Reference Williams, Huyvaert and Piaggio2016), from a variety of containers, either natural or artificial, that presented variation in water turbidity (defined as either clean or dirty; fig. 2, table 1). Eighteen (48.65%) of the water samples were taken in larval development sites where larvae were present, while 19 samples (51.35%) came from sites with no larvae detected. When larvae were visually detected, they were identified at least at the genus level (table 1), except for three samples for which a correct de visu identification was not possible.
Photograph representing the sampling methods used in our study: (A) sampling water in an artificial container, (B) sampling in a natural rock hole and (C) a CDC light trap in the oil palm plantation.
Characteristics of water samples (n = 37) collected on São Tomé Island (village n = 17; oil palm plantation n = 8; forest n = 12; with A or N for artificial or natural containers respectively) and the species identification, either visually or by metabarcoding (COI marker).
In five sequenced samples, we did not detect Culicidae species but other arthropod families (see fig. 3).
Stacked bar plots of the various arthropods detected along the anthropogenic gradient using eDNA metabarcoding (COI marker): (a) order level, (b) family level for the Diptera order and (c) species level for the family Culicidae (V, village; P, plantation; F, forest).
A total of 47 CDC light traps were set up to collect adult mosquitoes three consecutive nights in each habitat in parallel of the water sampling (fig. 2). Eleven traps were in the village, 18 in the oil palm plantation and 18 in the forested areas. Every morning, traps were gathered and placed in a freezer for 15 min. Then all arthropods were sorted and dipterans of interest were identified morphologically using a Leica S9E stereomicroscope (Leica Microsystems GmbH, Germany). Adults and larvae mosquito were identified to species or species group using different morphological keys and detailed descriptions are provided in Edwards (Reference Edwards1941), Hopkins (Reference Hopkins1952), Gillies and Coetzee (Reference Gillies and Coetzee1987), Service (Reference Service1990) and Ribeiro et al. (Reference Ribeiro, Da Cunha Ramos, Capela and Alves Pires1998). Our sources of data on species naming were based on that recorded in the Walter Reed Biosystematics Unit Mosquito Catalogue (http://www.mosquitocatalog.org).
Molecular methods
DNA extractions were performed in a laboratory ‘clean-room’ (at CIBIO, Portugal) equipped with UV radiation where strict protocols are followed for the prevention of contaminations (disposable laboratory clothing, UV sterilization of all equipment before entering the laboratory and laboratory cleaning with a 60% dilution of bleach between extraction batches). Prior to filtration, the water samples were manually shaken for 5 min (Civade et al., Reference Civade, Dejean, Valentini, Roset, Raymond and Bonin2016; Lopes et al., Reference Lopes, Sasso, Valentini, Dejean, Martins, Zamudio and Haddad2017) to homogenize the water column within the 50 ml Falcons tube. To concentrate material to a suitable volume for subsequent extraction, we filtered each sample (40 ml; water + Longmire) by pouring it into a sterile container (100-ml filtering cup; Nalgene Polysulphone Filter Holder with funnel, Thermo Scientific, USA) through sterile 47 mm nitrocellulose disc filters, 0.45 μm pore size (Whatman, UK), using a vacuum pump. The disc filters were cut into small pieces and placed in a 50 ml Falcon tube with 1.5 ml 3 M sodium acetate and 33 ml absolute ethanol for the water samples. These samples were placed in a carousel rotating shaker for 2 h at room temperature to homogenize the samples. Subsequently, the water samples were stored for 24 h at −20°C. Filter manipulation was performed with sterilized forceps between samples. Subsequently, the samples were centrifuged at 3184g for 45 min, at 10°C to recover the precipitated DNA and/or cell debris (Peixoto et al., Reference Peixoto, Chaves, Velo-Antón, Beja and Egeter2021). The supernatant was discarded (Valiere and Taberlet, Reference Valiere and Taberlet2000) and we performed DNA extraction on the pellet using the Dneasy Blood and Tissue Kit following the manufacturer's instructions (Qiagen, Hilden, Germany) (Gutiérrez-López et al., Reference Gutiérrez-López, Martínez-de la Puente, Gangoso, Soriguer and Figuerola2015). The pellet was exposed to enzymatic lysis using proteinase K in a carousel rotating shaker for 1 h at 56°C and the supernatant was spun through the column purification of DNA. We included a negative control (distilled water) in each set of extractions to monitor potential contaminations. The DNA was eluted in 80 μl of ultrapure sterilized MilliQ water. After extraction, DNA was quantified using the Qubit high-sensitivity dsDNA assay (Thermo Fisher Scientific). DNA metabarcoding libraries were prepared by amplifying a 200 bp fragment of the COI genomic region using the following primers: eCul-F (5′ GGRKCHGGDACWGGDTGAAC 3′) and eCul-R (5′ GATCAWACAAATAAAGGTAWTCGATC 3′) (Krol et al., Reference Krol, Van der Hoorn, Gorsich, Trimbos, Bodegom and Schrama2019). Illumina sequencing primer sequences were attached to the 5′ ends of PCR primers with i7 and i5 as indexes (Index 1 (i7) Adapter: P7-P5’CAAGCAGAAGACGGCATACGAGAT[i7]GTGACTGGAGTTCAGACGTGTGCTCTTCCGATC; Index 2 (i5) Adapter: P5-P7’AATGATACGGCGACCACCGAGATCTACAC[i5]ACACTCTTTCCCTACACGACGCTCTTCCGATCT). PCRs were carried out in a final volume of 25 μl, containing 2.5 μl of template DNA, 0.5 μM of each primer, 12.5 μl of Supreme NZYTaq 2× Green Master Mix (NZYTech) and ultrapure water up to 25 μl. The thermocycler programme for DNA amplification started with an initial denaturation step at 95°C for 5 min, followed by 40 cycles of 95°C for 30 s, 58°C for 45 s, 72°C for 30 s and a final extension step at 72°C for 10 min.
The oligonucleotide indices, which are required for multiplexing different libraries in the same sequencing pool, were attached in a second PCR round with identical conditions but for only ten cycles and 60°C as the annealing temperature. We used in-house designed indexes, which are a combinatorial set of 24 i5 and 24 i7 indexes, which we have pre-mixed and randomized. They are 8-bp long and the Levenshtein distance between any two indexes is at least 3. A negative control containing no DNA (MilliQ water) was included in every PCR round to check for contamination during library preparation. The libraries were run on 2% agarose gels stained with GreenSafe (NZYTech), and imaged under UV light to verify the library size. Libraries were purified using the Mag-Bind RXNPure Plus magnetic beads (Omega Biotek). We pooled the samples and purified the resulting pool following the same method (1× of magnetic beads). The purified pool was run through a size-select eGel to precisely select the band of interest. Libraries were quantified using the Qubit high-sensitivity dsDNA assay (Thermo Fisher Scientific).
Very low library quantification was detected in 18 samples that were removed for sequencing. These samples corresponded to water samples in which only one larva (n = 2) or none were detected visually (n = 16; table 1). Therefore, 19 samples were selected for sequencing and were pooled in equimolar amounts and re-purified for double size selection in an e-gel system (Life Technologies) for primer dimer elimination. The pool was sequenced in a fraction (1/16) of a MiSeq PE300 run (Illumina). Library preparation and sequencing were carried out by AllGenetics & Biology SL (www.allgenetics.eu).
Bioinformatics analyses and taxonomic assignment
Illumina paired-end raw files consist of forward (R1) and reverse (R2) reads sorted by library and their quality scores. The indices and sequencing primers were trimmed from the samples using software CUTADAPT (Martin, Reference Martin2011) and the quality of the FASTQ files was checked using the software FastQC (Andrews, Reference Andrews2010). Plots summarizing the quality across bases of R1 and R2 reads were generated by using MultiQC (Ewels et al., Reference Ewels, Magnusson, Lundin and Käller2016) (see Supplementary file). The merging of the R1 and R2 reads was performed with FLASH2 (Magoč and Salzberg, Reference Magoč and Salzberg2011). The mismatch resolution in the overlapping region (minimum overlap of 30 base pairs) was accomplished by keeping the base with the higher quality score. We used CUTADAPT 1.3 software (Martin, Reference Martin2011) to remove sequences that did not contain the PCR primers (allowing up to two mismatches) and sequences that ended up being shorter than 145 nucleotides and larger than 210 nucleotides. The sequences were quality-filtered (minimum Phred quality score of 20), then were dereplicated (-derep fulllength) and clustered at a similarity threshold of 97% (-cluster fast, -centroids option) and sorted (-sortbysize) using VSEARCH (Rognes et al., Reference Rognes, Flouri, Nichols, Quince and Mahé2016). De novo chimaera detection was carried out using the UCHIME algorithm (Edgar et al., Reference Edgar, Haas, Clemente, Quince and Knight2011) implemented in VSEARCH.
We conducted the taxonomic assignment of each operational taxonomic unit (OTU) (I) using a customized taxonomic COI reference database. The database included (i) newly generated mosquito sequences of four species sampled during the fieldwork using light traps (Aedes nigricephalus, Culex cambournaci, Uranotaenia micromelas and Uranotaenia connali; Genbank accession numbers ON504276–ON504279), and (ii) sequences downloaded from the National Center for Biotechnology Information (NCBI) and the Barcode of Life Data System (BOLD) databases (Ratnasingham and Hebert, Reference Ratnasingham and Hebert2007) (accessed on March 2022). These mosquito sequences (from mosquito species known to be present on São Tomé; table S1) were added to the database build using RESCRIPt (Robeson et al., Reference Robeson, Rourke, Kaehler, Ziemski, Dillon, Foster and Bokulich2021) (last version on July 2020) based on the BOLD reference database (Ratnasingham and Hebert, Reference Ratnasingham and Hebert2007).
We employed the script feature-classifier classify-consensus-vsearch implemented in Qiime2 (Bokulich et al., Reference Bokulich, Kaehler, Rideout, Dillon, Bolyen, Knight, Huttley and Gregory2018) and the VSEARCH algorithm (Rognes et al., Reference Rognes, Flouri, Nichols, Quince and Mahé2016) with a sequence similarity threshold of 95%. In addition, we used the top-hits-only option in the VSEARCH command to recover only the hit with the highest percentage of identity. In spite of the multiple top hits used in the consensus taxonomic assignment carried out by VSEARCH, this option allows the assignment of the query to the closest reference sequence. The table resulting from this step lists the number of sequences from each OTU found in each sample and their corresponding taxonomic information (table S2 – before OTU filtering). Subsequently, based on the results of this table, we applied several different filters. We removed singletons (i.e., OTUs containing only one-member sequence in the whole data set). In DNA metabarcoding studies, it has been observed that a low percentage of the reads of a library can be assigned to another library. This phenomenon, referred to as mistagging, tag jumping, index hopping, index jumping, etc. is the result of the misassignment of the indices during library preparation, sequencing and/or demultiplexing steps (Esling et al., Reference Esling, Lejzerowicz and Pawlowski2015; Bartram et al., Reference Bartram, Mountjoy, Brooks, Hancock, Williamson, Wright, Moppett, Goulden and Hubank2016; Guardiola et al., Reference Guardiola, Wangensteen, Taberlet, Coissac, Uriz and Turon2016; Illumina, Reference Illumina2017). In order to correct for this phenomenon, OTUs occurring at a frequency below 0.01% in each sample were removed. Finally, only the OTUs that matched any reference sequence in the database at a minimum similarity threshold of 85% were kept in the OTU table. Therefore, the unidentified OTUs (‘Unassigned’) were removed from the OTU table for downstream analysis (table S3 – after OTUs filtering). Six samples (V16, V17, P4, F1, F7, F9) had no OTUs assigned to the family Culicidae.
The alpha rarefaction plots show the number of OTUs obtained with a rarefied number of sequences in each sample. These plots were generated using the OTU table before (table S2) and after (table S3) the OTU filtering (fig. S1). The vertical axis displays the number of OTUs observed at different subsampling depths. When the rarefaction curves tend towards saturation, the sequencing depth is considered to be sufficient to retrieve most of the taxa diversity. We have to note that curve from sample V8 did not reach the plateau in the number of OTUs observed (see fig. S1.B, rarefaction plot after the OTU filtering).
In order to easily visualize the breakdown of taxonomic classification, stacked bar plots showing the relative abundance of each OTU in each sample were generated at the order, family and species levels (fig. 3). In DNA metabarcoding studies, OTU relative abundance is defined as the number of reads assigned to that OTU divided by the total number of reads. Note that the PCR may cause biases due to differences in primer specificity. These biases can cause taxa with low representation in the original DNA sample to become more abundant in the final results. As a result, this bias prevents from correctly inferring the abundance of species in the original DNA sample. For example, if SPECIES A is represented by the 35% of the sequences in SAMPLE 1, and SPECIES B is represented by the 50% of the sequences in the same sample, we cannot reliably conclude that there was more SPECIES B DNA in the original sample. That being said, it is expected that, within the same study, the PCR bias always go in the same direction. Therefore, it is possible to compare how the abundance of a given taxon varies across different samples with a similar composition. For example, if SPECIES A is represented by the 35% of the sequences in SAMPLE 1 and by the 10% in SAMPLE 2, we can conclude that there was less SPECIES A DNA in SAMPLE 2 (Geisen et al., Reference Geisen, Laros, Vizcaíno, Bonkowski and De Groot2015; Thomas et al., Reference Thomas, Deagle, Eveson, Harsch and Trites2016; Matesanz et al., Reference Matesanz, Pescador, Pías, Sánchez, Chacón-Labella, Illuminati, de la Cruz, López-Angulo, Marí-Mena and Vizcaíno2019).
Finally, we extracted the representative sequences for each of the picked OTUs before and after the OTU filtering process. For the particular case of the taxonomic assignment of OTUs to Eretmapodites intermedius, we performed a blast in NCBI and the results are shown in fig. S2.
DNA metabarcoding analyses were carried out by AllGenetics & Biology SL (www.allgenetics.eu).
Results
Visual and genetic detection
Of the 19 water samples collected from sites where no larvae were detected visually, one was positive for Ae. albopictus (5%; table 1), two others were found with chironomids (Diptera) or coleopterans (10%) and 16 could not be sequenced because of the low library DNA quantities (84%). Of the 18 water samples in which larvae were seen, eDNA metabarcoding detected Culicidae in 13 (72%), three of which had detections of other dipterans and branchiopodans (16%), and two could not be sequenced because of the low library DNA quantities (11%; table 1). When larvae were present at the collection site, one or two Culicidae genera were identified visually in each sample, whereas eDNA metabarcoding detected up to four genera per sample (table 1).
We recovered DNA sequences in 14 water samples out of the 26 considered as clean (53.8%), and in four out of seven considered as dirty (57%). Although our sample sizes remain small, we found that the turbidity of the water did not appear to be a limitation for eDNA metabarcoding (chi-squared test, χ 2 = 0.33).
Overall, the taxonomic assignments revealed four orders of arthropods that comprised of 13 families. Within Culicidae, taxonomic assignments at the species level for the genus Anopheles returned An. coluzzii, the main human malaria vector on the island (Chen et al., Reference Chen, Lien, Tseng, Cheng, Lin, Wang and Tsai2019). For the genus Aedes, the taxonomic assignments at the species level returned the invasive tiger mosquitoes Ae. albopictus and Ae. aegypti. All OTUs that matched the genus Eretmapodites, an endemic genus of the Afrotropical region and vector of various viruses (Bamou et al., Reference Bamou, Mayi, Djiappi-Tchamen, Nana-Ndjangwo, Nchoutpouen, Cornel, Awono-Ambene, Parola, Tchuinkam and Antonio-Nkondjio2021), were assigned to Er. intermedius (Supplementary fig. S2). As for the Culex genus, OTUs were assigned to Cx. cambournaci, Cx. decens and Cx. sasai.
In summary, 12 species of Culicidae were detected, seven with eDNA metabarcoding and nine with CDC light traps. Four species were common to both approaches: Ae. albopictus, An. coluzzii/gambiae, Cx. Cambournaci and Cx. decens, all collected in the village (table 2; fig. 4).
List of mosquito species by habitat recovered using both methods: CDC traps (depicted by the adult mosquito) and eDNA metabarcoding (depicted by the water bottle). Species detected uniquely with eDNA metabarcoding are in bold (see also table 2). Icons from Freepik.
List of mosquito species detected in the water samples, visually and with eDNA metabarcoding, and with CDC traps along the gradient of anthropogenic disturbance in São Tomé Island
See also fig. 4 for illustration.
a Incorrect taxonomic assignment likely due to incomplete molecular reference database.
Habitat effects on species detection
In the village, five orders and eight families of arthropods were found. The Culicidae was the dominant family found in the village, with 78% of the total reads from the village attributed to the genera Aedes, Anopheles, Culex and Eretmapodites. The invasive mosquito Ae. albopictus and the malaria vectors An. coluzzii were present respectively in 57% (N = 8) and 50% (N = 7) of the samples collected in the village that led to amplification. Ae. albopictus was found in both artificial and natural breeding sites, while Ae. aegypti was totally absent from the village, a pattern that had already been noted in previous surveys (Reis et al., Reference Reis, Cornel, Melo, Pereira and Loiseau2017). Culex spp. were present in half of the village samples that could be sequenced (7 out of 14; fig. 3).
In the plantation, in the eight potential breeding sites that were sampled, we did not detect any larvae visually. The only sample whose amplification worked gave two OTUs affiliated to the Chironomidae family (order Diptera; see Supplementary tables S2 and S3).
In the forest, four orders and four families of arthropods were found, with the Chironomidae being the dominant family with 73% of the reads (fig. 3). In the forest, Cx. sasai and Ae. aegypti were detected in the same sample (fig. 3).
Discussion
Our study showed that eDNA metabarcoding could be a complementary method to the light or decoy traps to recover mosquito diversity, and help to evaluate the assemblage of species using the same breeding sites. In particular, eDNA metabarcoding was able to detect species that were not captured with light traps and picked up different assemblage of mosquito species associated with the degree of anthropogenic disturbance.
In the oil palm plantation, we found larvae of mosquitoes by de visu at one sampling location. eDNA metabarcoding detected only one family of DIPTERA (Chironomidae) with very few reads, but no mosquito species. This result is not surprising and is consistent with the view that the core of oil palm plantations is overall poor in terms of arthropod diversity (Koh and Wilcove, Reference Koh and Wilcove2008; Turner and Foster, Reference Turner and Foster2009; Fayle et al., Reference Fayle, Turner, Snaddon, Chey, Chung, Eggleton and Foster2010; Ghazali et al., Reference Ghazali, Asmah, Syafiq, Yahya, Aziz, Tan, Norhisham, Puan, Turner and Azhar2016). Recently, Young et al. (Reference Young, Buenemann, Vasilakis, Perera and Hanley2021) also found that mosquito abundance in oil palm plantations in Borneo was lower than that in the forest. On the contrary, in the village, the arthropod diversity was much higher than that in the surrounding plantations with eight families of Diptera recorded. Culicidae was the predominant family: Ae. albopictus accounted for 36% of the reads, followed by Culex species (33.5%), while Anopheles genus was the least abundant, with 3.3% of the reads. Although more surveys are needed, Ae. albopictus, which recently colonized the island (Reis et al., Reference Reis, Cornel, Melo, Pereira and Loiseau2017), shared breeding sites with Culex, Eretmapodites and Anopheles species. Co-occurrence with the latter was less expected since these species do not usually use the same niche. Finally, in the forest, among the four families of Diptera detected, Chironomids were the predominant one, with 73% of the reads, while mosquito species were found in lower abundance (17%). Interestingly, the yellow fever mosquito Ae. aegypti was detected in only one sample, inside a bamboo stalk. It used to be very common and widespread on the island, and found equally in both natural and artificial breeding sites (Ribeiro et al., Reference Ribeiro, Da Cunha Ramos, Capela and Alves Pires1998). However, recent on-going mosquito projects and, surveys on the island revealed that Ae. aegypti became quite rare and seems to have been replaced by Ae. albopictus in lowland and disturbed habitats (Reis et al., Reference Reis, Cornel, Melo, Pereira and Loiseau2017; Loiseau et al., Reference Loiseau, Gutiérrez-López, Mathieu, Makanga, Paupy, Rahola and Cornel2022). This replacement pattern has been largely documented in Florida, USA (Yang et al., Reference Yang, Borgert, Alto, Boohene, Brew and Deutsch2021) but is less evident in mainland Central Africa (Simard et al., Reference Simard, Nchoutpouen, Toto and Fontenille2005; Paupy et al., Reference Paupy, Ollomo, Kamgang, Moutailler, Rousset, Demanou, Hervé, Leroy and Simard2010; Kamgang et al., Reference Kamgang, Ngoagouni, Manirakiza, Nakouné, Paupy and Kazanji2013; Tedjou et al., Reference Tedjou, Kamgang, Yougang, Njiokou and Wondji2019). Nonetheless, our eDNA metabarcoding approach corroborates the actual known distribution of these two Aedes species on the island (Loiseau et al., Reference Loiseau, Gutiérrez-López, Mathieu, Makanga, Paupy, Rahola and Cornel2022). Finally, the other Culicidae species found in the forest was Cx. sasai. It is highly unlikely that this mosquito is present on the island, since to date it has been detected only in Asia (Phanitchakun et al., Reference Phanitchakun, Wilai, Saingamsook, Namgay, Drukpa, Tsuda, Walton, Harbach and Somboon2017), and is not known to be present on São Tomé Island (Loiseau et al., Reference Loiseau, Gutiérrez-López, Mathieu, Makanga, Paupy, Rahola and Cornel2022). Because Cx. sasai belongs to the Culiciomyia subgenus, we probably detected here a mosquito species belonging to this same subgenus. There are actually four species of this subgenus on São Tomé Island: Cx. cambournaci, Culex nebulosus, Culex cinerellus and Culex macfiei (Loiseau et al., Reference Loiseau, Gutiérrez-López, Mathieu, Makanga, Paupy, Rahola and Cornel2022), with only two having barcoding sequences on online databases (Cx. cambournaci and Cx. nebulosus). One could speculate that the species found in this forest sample could be either Cx. cinerellus or Cx. macfiei and not Cx. sasai. This error highlights one of the limitations of the eDNA metabarcoding approach which is discussed below, i.e., incomplete reference databases.
Challenges of eDNA metabarcoding: sample quality and taxonomic assignment issues
As with any new methods, some weaknesses and concerns need to be addressed. Some critical factors for the application of eDNA methods to detect aquatic species have already been reviewed (Goldberg et al., Reference Goldberg, Turner, Deiner, Klymus, Thomsen and Murphy2016), including contamination in the field and in the laboratory, choosing appropriate sample analysis methods, validating assays or testing for sample inhibition. Here, we highlight concerns that are specific for insect vector monitoring using eDNA approaches.
First, mosquito larvae are mostly found in small and turbid breeding sites or in stagnant water bodies. While water from some larval breeding sites (e.g., rock pools, puddles, artificial containers) is easy to sample, it can be difficult to collect from other sites (e.g., tree holes, plant axils). Traps and sampling procedures, such as aspiration of resting mosquitoes, collection on human or animal bait, allow collecting a greater diversity of species. For inventory purposes, eDNA techniques may need a great water sampling effort in order to be comparable to other techniques (Krol et al., Reference Krol, Van der Hoorn, Gorsich, Trimbos, Bodegom and Schrama2019). In addition, sampling small volumes of water can lead to false-negative detection when the density of targeted organisms is low (Ulibarri et al., Reference Ulibarri, Bonar, Rees, Amberg, Ladell and Jackson2017). Another potential sampling issue is the large amount of soil and humic substances found in breeding sites that may act as PCR inhibitors, increasing the chance to obtain false-negative results (Buxton et al., Reference Buxton, Groombridge and Griffiths2017). In our case, we managed to amplify COI even from dirty samples, although these samples contained many larvae. One study experimentally tested the success of PCR detection of eDNA samples from containers with two different water volumes (50 ml and 1 litre) (Odero et al., Reference Odero, Gomes, Fillinger and Weetman2018). They found that the volume of water required in relation to the density of larvae has an effect on the mosquito detection by eDNA analysis. The detection was better when the samples had many larvae at low densities than few larvae at higher densities (Odero et al., Reference Odero, Gomes, Fillinger and Weetman2018). In addition, the effect of different substrates in the eDNA analysis as well as the preservation methods are parameters that should not be overlooked since metabarcoding analyses require good DNA quality (Ball and Armstrong, Reference Ball and Armstrong2014).
Second, it seems appealing to evaluate and compare mosquito diversity from different types of samples (water vs. bulk samples) using the metabarcoding approach because traditional dipping methods to survey larvae in breeding sites may not always reflect the adult diversity that can be found with CDC traps (and inversely). In fact, in our survey, only four species were shared between the two techniques (eDNA vs. CDC traps). It is worth noting that some species may be very difficult to detect with traditional trapping because not all insect vector species are equally attracted to dry ice or light (Reisen and Lothrop, Reference Reisen and Lothrop1999). It is especially true for daytime biting mosquitoes. On the other hand, it might be difficult to sample water in breeding sites, such as plant axils or tree holes, which can be high up. More investigations under controlled conditions are needed to compare the efficacy of metabarcoding water samples with trapped adults to characterize insect-vector communities.
Third, in the BOLD, of about 3500 species of Culicidae known globally, barcodes are only available for 1329 species (38%; accessed on 2021-05-25) and, among the 41 known mosquito genera, three genera alone (Aedes, Anopheles and Culex) account for 78% of the occurrences. Similar patterns are found when gathering data on different genes in NCBI (COI, 18S rRNA and 28S rRNA). While Aedes, Culex and Anopheles species account for only 60% of the total mosquito species, 90% of the sequences on average correspond to these three genera (see fig. S3 for illustration of these data). Sequences belonging to unknown taxa are still a common problem in eDNA barcoding and therefore, when starting a new monitoring programme to assess the mosquito diversity in a region or locality, creating a good quality reference database is an indispensable first step. This means that a considerable amount of essential taxonomic work is required to setup eDNA-based monitoring protocols. In this study, we managed to get DNA sequences of four mosquito species that were not deposited in online databases yet. Eleven species out of the 34 known on the island (Loiseau et al., Reference Loiseau, Gutiérrez-López, Mathieu, Makanga, Paupy, Rahola and Cornel2022) still have to be captured and sequenced to have a full reference database for future research work. Taking all this into account, and considering that certain limitations can be surpassed, then eDNA metabarcoding can have significant advantages for mosquito surveys.
Advantages of eDNA metabarcoding: easy sampling and less entomological expertise required
Sampling for eDNA can be as simple as collecting freshwater samples in tubes and adding preservation buffers (Williams et al., Reference Williams, Huyvaert and Piaggio2016), which drastically reduces the cost and time allocated to fieldwork, as well as equipment and resources required for sampling. This is particularly relevant for research projects carried out in remote regions. The effort required for the traditional trapping methods is substantial. Logistically it requires the transport of traps and batteries (which are voluminous and heavy), the availability of freezers (to kill mosquitoes before identification) and of high-quality stereomicroscopes. Once this material is in the field, traps must be set up for several hours, with light that attracts mosquitoes together with a wide range of flying insects, or with traps containing odour products to attract more specifically females (BG-Sentinel or Gravid Mosquito traps). Since light traps are not selective, a great amount of time is spent on sorting all the flying insects from the mosquitoes, separating engorged individuals and labelling individual tubes. Once back in the laboratory, experts may spend a great amount of time at the microscope identifying and dissecting individuals. Identification of mosquito eggs and larvae implies mounting, which is time consuming, and require a specific training. Although an alternative solution could be rearing larvae into adults for unambiguous identification, this is logistically challenging when doing fieldwork in remote places. In addition, for the identification of many adult insect vectors, dissecting male genitalia is required, which is the case for example for most of the species of the African genus Eretmapodites (Service, Reference Service1990). Molecular identification of eDNA is able to circumvent time-consuming morphological investigation and to detect the presence of species without requiring a strong entomology expertise. The efficacy of eDNA-based surveys will increase as reference databases become more complete. Interestingly, in our study, we detected the species Er. intermedius for the first time on the island, as until now Eretmapodites chrysogater was the only known representative of this genus on the island (Ribeiro et al., Reference Ribeiro, Da Cunha Ramos, Capela and Alves Pires1998). This detection would have been almost impossible using traditional light traps since Eretmapodites species are day-biting mosquitoes and males are generally less attracted to them. Finally, the ease of water sampling procedures for eDNA protocols will allow developing large-scale citizen science monitoring programmes and integrating non-specialists in research projects (Biggs et al., Reference Biggs, Ewald, Valentini, Gaboriaud, Dejean, Griffiths, Foster, Wilkinson, Arnell, Brotherton, Williams and Dunn2015).
Concluding remarks
To date, numerous studies have demonstrated that eDNA sampling generally provides greater detection probabilities than traditional techniques (Thomsen et al., Reference Thomsen, Kielgast, Iversen, Wiuf, Rasmussen, Gilbert, Orlando and Willerslev2012; McKelvey et al., Reference McKelvey, Young, Knotek, Carim, Wilcox, Padgett-Stewart and Schwartz2016; Valentini et al., Reference Valentini, Taberlet, Miaud, Civade, Herder and Thomsen2016), but it still remains to be formally demonstrated for mosquito communities. In fact, eDNA methods could surely help in applied medical and veterinary entomology and significantly improve (i) the detection of invasive species and (ii) the evaluation of the composition of mosquito communities in understudied regions. In our study, we showed that CDC light traps and adult identification methods recovered more species than the eDNA metabarcoding per habitat. However, eDNA metabarcoding was able to detect (i) more species at a mosquito breeding site than de visu larval identification, and (ii) different species than traditional methods. Therefore, our results highlight the fact that it is best to use in conjunction traditional survey methods and eDNA metabarcoding to enhance detection rates and increase confidence in the monitoring results.
Like any ecological survey tool, eDNA metabarcoding will always suffer biases and uncertainties which have to be taken into account at each step of the study (i.e., fieldwork, labwork, bioinformatics analyses). The building up of the BOLD is required to expand the potential of eDNA metabarcoding, a task where taxonomic expertise will be essential. However, the relative simplicity of field sampling protocols can create opportunities to collect samples using volunteers and even to develop citizen science programmes such as (i) for monitoring and surveillance of invasive species such as Ae. albopictus, and (ii) for improving our understanding of ecological systems (competition and predation at breeding sites) that could definitely help in vector control management (Dambach, Reference Dambach2020).
Supplementary material
The supplementary material for this article can be found at https://doi.org/10.1017/S0007485323000147.
Acknowledgements
We are grateful to the field assistants on São Tomé: Ricardo ‘Mito’ Fonseca, Martim Veiga and Sidney ‘Dulay’ Samba, and we thank Arlindo Carvalho, former Director of the Department of the Environment of São Tomé and Príncipe for granting us the permits to conduct the research. We thank Antón Vizcaíno, Ania Pino-Querido and Neus Marí-Mena for their work in the lab and with the bioinformatics analyses. We also thank the two anonymous reviewers for their constructive comments on the manuscript.
Financial support
This work is funded by through FCT – Foundation for Science and Technology (Portugal) under the PTDC/BIA-EVL/29390/2017 DEEP Research Project (C.L.) and via structural funding for CIBIO-InBIO (UIDB/50027/2021). R.G.L. was supported by the project PTDC/BIA-EVL/29390/2017 DEEP Research Project (FCT – Foundation for Science and Technology (Portugal)) and by the grant Juan de la Cierva 2019 Formación from the Ministry of Science and Innovation (REF: FJC2019-041291-I). B.E. was supported via the European Union's Horizon 2020 research and innovation programme under grant agreement No. 668981. C.P., N.R. and D.J. were supported by the French National Research Agency (ANR PRC TIGERBRIDGE, grant No. 16-CE35-0010-01). M.M. was supported via the European Union's Horizon 2020 research and innovation programme under grant agreement No. 854248.
