Monitoring mosquito richness in an understudied area: can environmental DNA metabarcoding be a complementary approach to adult trapping?

Abstract

Mosquito surveillance programmes are essential to assess the risks of local vector-borne disease outbreaks as well as for early detection of mosquito invasion events. Surveys are usually performed with traditional sampling tools (i.e., ovitraps and dipping method for immature stages or light or decoy traps for adults). Over the past decade, numerous studies have highlighted that environmental DNA (eDNA) sampling can enhance invertebrate species detection and provide community composition metrics. However, the usefulness of eDNA for detection of mosquito species has, to date, been largely neglected. Here, we sampled water from potential larval breeding sites along a gradient of anthropogenic perturbations, from the core of an oil palm plantation to the rainforest on São Tomé Island (Gulf of Guinea, Africa). We showed that (i) species of mosquitoes could be detected via metabarcoding mostly when larvae were visible, (ii) larvae species richness was greater using eDNA than visual identification and (iii) new mosquito species were also detected by the eDNA approach. We provide a critical discussion of the pros and cons of eDNA metabarcoding for monitoring mosquito species diversity and recommendations for future research directions that could facilitate the adoption of eDNA as a tool for assessing insect vector communities.

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, 2000; Segan et al., 2016; Eriksson and Hillebrand, 2019), with important consequences for the emergence of infectious diseases that affect wildlife, livestock and human populations (Caminade et al., 2019; Smith et al., 2019). Among global emerging infectious disease events, vector-borne diseases are disproportionately over-represented (Swei et al., 2020) and constant efforts for monitoring insect vector populations should be carried out in locations at risk (Pedersen et al., 2009). Among insect vectors, mosquitoes (Culicidae), with more than 3500 described species widely distributed around the world (Harbach, 2013), 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., 2013; Iwamura et al., 2020). 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, 2009).

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, 2004). Skilled entomologists are able to identify specimens using morphological traits (Besansky et al., 2003; Hajibabaei et al., 2007), however some species are indistinguishable morphologically (e.g., cryptic species of Anopheles) (Coetzee and Koekemoer, 2013). 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., 2007). 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., 2008; Goldberg et al., 2011; Jerde et al., 2011; Minamoto et al., 2012; Thomsen et al., 2012; Spear et al., 2015; Egeter et al., 2018), and more recently for invasive invertebrates (Clusa et al., 2017; Klymus et al., 2017; Mychek-Londer et al., 2020). In natural habitats, eDNA is affected by a variety of factors, such as temperature, microbial activity, pH (Seymour et al., 2018), conductivity (Collins et al., 2018), 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., 2011; Barnes et al., 2014; Pilliod et al., 2014). Hence, most eDNA detection is expected to indicate a current or recent colonization of the habitat (Piaggio et al., 2014), 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. (2016) 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., 2019; Krol et al., 2019). These studies both used eDNA primers targeting the cytochrome oxidase I (COI) gene. Boerlijst et al. (2019) 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., 2019). 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., 2019). 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., 2019). 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., 2012; Taberlet et al., 2012; Deiner et al., 2017; Hering et al., 2018). 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., 2012).

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., 2017) 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).

Figure 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., 2016), 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.

Figure 2. 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.

Table 1. 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).

Figure 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 (1941), Hopkins (1952), Gillies and Coetzee (1987), Service (1990) and Ribeiro et al. (1998). 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., 2016; Lopes et al., 2017) 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., 2021). The supernatant was discarded (Valiere and Taberlet, 2000) 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., 2015). 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., 2019). 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, 2011) and the quality of the FASTQ files was checked using the software FastQC (Andrews, 2010). Plots summarizing the quality across bases of R1 and R2 reads were generated by using MultiQC (Ewels et al., 2016) (see Supplementary file). The merging of the R1 and R2 reads was performed with FLASH2 (Magoč and Salzberg, 2011). 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, 2011) 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., 2016). De novo chimaera detection was carried out using the UCHIME algorithm (Edgar et al., 2011) 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, 2007) (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., 2021) (last version on July 2020) based on the BOLD reference database (Ratnasingham and Hebert, 2007).

We employed the script feature-classifier classify-consensus-vsearch implemented in Qiime2 (Bokulich et al., 2018) and the VSEARCH algorithm (Rognes et al., 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., 2015; Bartram et al., 2016; Guardiola et al., 2016; Illumina, 2017). 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., 2015; Thomas et al., 2016; Matesanz et al., 2019).

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, χ ² = 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., 2019). 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., 2021), 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).

Figure 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.

Table 2. 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., 2017). 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, 2008; Turner and Foster, 2009; Fayle et al., 2010; Ghazali et al., 2016). Recently, Young et al. (2021) 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., 2017), 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., 1998). 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., 2017; Loiseau et al., 2022). This replacement pattern has been largely documented in Florida, USA (Yang et al., 2021) but is less evident in mainland Central Africa (Simard et al., 2005; Paupy et al., 2010; Kamgang et al., 2013; Tedjou et al., 2019). Nonetheless, our eDNA metabarcoding approach corroborates the actual known distribution of these two Aedes species on the island (Loiseau et al., 2022). 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., 2017), and is not known to be present on São Tomé Island (Loiseau et al., 2022). 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., 2022), 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., 2016), 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., 2019). In addition, sampling small volumes of water can lead to false-negative detection when the density of targeted organisms is low (Ulibarri et al., 2017). 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., 2017). 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., 2018). 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., 2018). 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, 2014).

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, 1999). 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., 2022) 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., 2016), 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, 1990). 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., 1998). 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., 2015).

Concluding remarks

To date, numerous studies have demonstrated that eDNA sampling generally provides greater detection probabilities than traditional techniques (Thomsen et al., 2012; McKelvey et al., 2016; Valentini et al., 2016), 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, 2020).