Culicidae-centric metabarcoding through targeted use of D2 ribosomal DNA primers

Testing primer coverage for mosquito metabarcoding

We mined GenBank (15/11/2018) for potential primer sequences that would maximize recovery of mosquito amplicons without bias and maximally exclude non-culicid insects. Initially, we selected only GenBank entries with at least a 30 bp on either side of the putative oligo to avoid database submissions that might also include original amplification-primer sequences.

From the available sequences, we identified a primer set that produced an amplicon from the D2 segment of the 28S operon. The forward primer (Mozzie.D2.Uni.F; 5′-AAGCACTCTGAATAGAGAGTC-3′) was designed at nucleotide position 67 and the reverse (Mozzie.D2.Uni.R; 5′-TGGTCCGTGTTTCAAGAC-3′) at nucleotide position 489 of the Aedes aegypti ribosomal DNA locus (GenBank accession no. MG242540).

We next tested binding site variation and the taxonomic resolution of Mozzie.D2.Uni.F and Mozzie.D2.Uni.R with only those GenBank sequences that possessed both primers (sequences initially used for primer design did not necessarily contain both priming sites) using the ecopcr module of Obitools (v1.2.13) (Boyer et al., 2016). For this analysis, we did not use the criterion for a 30 bp flank above because we sought primarily data on taxonomic resolution of the resultant amplicons.

We then compared these D2 results with primer sets previously used in mosquito metabarcoding studies. These included two sets from mitochondrial 16S and two from CO1 (Schneider et al., 2016; Talaga et al., 2017; Batovska et al., 2018; Krol et al., 2019).

We compared the primer taxonomic coverage and estimated taxonomic resolution of the resultant amplicons in culicids and non-culicids. The latter taxa were divided into separate groups based on our own experience of the bycatch content in commonly used mosquito traps and on relevant literature (Pucci, Lynch & Keiper, 2003; Medlock et al., 2018). Thus, in addition to Culicidae, primers were evaluated for non-culicid Nematocera, Hymenoptera, Lepidoptera, Coleoptera and all remaining arthropods.

We ran ecopcr with a maximum of five nucleotide differences allowed in either primer and considered the threshold for a primer set’s successful amplification in a given taxonomic group as the presence or absence of the matching 3′ nucleotide in the oligos. The resulting ecopcr output was analyzed and graphed in R (R Development Core Team 3.0.1, 2013) using the ROBITools package (http://metabarcoding.org/obitools).

Mosquito DNA and mock assemblages

Mosquitoes used in pools were collected at various times either in Sao Paulo or Rondônia states, Brazil, and chosen to represent as fully as possible the general Culicidae taxonomic topology. Species used in the mock pools included one individual each of the following taxa: Aedeomyia squamipennis, Aedes sp., Anopheles (Kerteszia) sp., Anopheles (Nyssorhynchus) sp., Coquillettidia venezuelensis, Culex (Culex) sp., Culex sp., Haemagogus sp., Limatus sp., Mansonia flaveola, Mansonia humeralis, Mansonia indubitans, Mansonia titillans, Psorophora sp., Wyeomyia (Phoniomyia) sp., and Wyeomyia sp. A Chagasia sp. DNA sample from the culicid genomic collection at the Faculty of Public Health (University of Sao Paulo) was also used. In all cases, DNA was extracted using a simple salting out protocol preceded by maceration and proteinase K digestion (Sambrook, Fritsch & Maniatis, 2001). Most samples possessed varying degrees of pre- and post-extraction DNA sheering, as verified by agarose electrophoreses (due to sample age or capture method). We therefore did not explicitly measure sample DNA concentration for individual DNA stocks as this would be meaningless in light of degradation. Rather, we utilized the same stocks across the mock pools, thus obviating the need to standardize concentrations.

The DNA from one mosquito was diluted into approximately 300-µL of water and this working stock was used for the creation of four mock genomic DNA assemblies containing only culicid DNA (hereafter Fresh Pools A–D). These mixtures were comprised of varying volumes from the 17 working stocks (ranging in contribution from 0.5 to 40 µL of each animal; Data S1). The absolute volume of component taxa varied from 0.28% to 27% of the mock solutions, i.e., two orders of magnitude. The final DNA concentration of the four Fresh Pools composite mixtures averaged 0.84-ng/µL. A replicate sequencing run for the Fresh Pool A mixture, using different PCR and Illumina sequencing runs, was done to verify the reproducibility of the results (although at a six-fold lower NGS sequencing coverage).

We tested the extent of non-target amplification from the D2 primer pair (described below) by mixing the four Fresh Pools with DNA from the combined bycatch of four CDC traps after the maximum care had been taken to remove remnant Culicidae and their tissue, such as loose scales and antennae. The ad hoc identification of CDC content estimated more than 5,000 individuals of varying size (with the bulk of the catch being 5-mm or less) and from the insect orders Blattodea, Coleoptera, Diptera, Hemiptera, Hymenoptera, Lepidoptera, Orthoptera, and several spiders.

Here, the bycatch DNA was extracted as above and used at a concentration of 135 ng/µL in all downstream assays. The CDC bycatch was mixed with of the four Fresh Pools in one of two concentrations: 1:1 and 1:10 (1-µL of each Fresh Pool to 10-µl of CDC bycatch DNA). Of the resulting eight Bycatch Pools, Bycatch Pool A1:1 and Bycatch Pool A1:10 were submitted to a replicated PCR and sequencing pipeline.

A third set of mock assemblies (named Degraded Pools I–IV) was created to test the impact of DNA degradation on the amplification of the D2. These used the same stock genomic DNA solutions for each of the 17 animals as above, after having been left at room temperature for two weeks (an average duration for the field deployment of commonly used attraction traps). Evaporated water was replaced in order to maintain the original concentration. There was also a two orders of magnitude difference between the highest and lowest individual volumes used in each pool (0.5 to 50 µL; Data S1). The degraded DNA mixture averaged 0.74-ng/µL and omitted sample Wyeomyia sp.

PCR parameters

Identical PCRs profiles were used to sequence both the 17 specimens alone and the pooled assemblages. These used Illumina overhang sequences (underlined below) on each primer in a first round PCR:

 Ill+Mozzie.D2.Uni.F (5′-TCGTCGGCAGCGTCAGATGTGTATAAGAGACAGAAGCACTCTGAATAGAGAGTC-3′) andIll+Mozzie.D2.Uni.R (5′-GTCTCGTGGGCTCGGAGATGTGTATAAGAGACAGTGGTCCGTGTTTCAAGAC-3′). The second-round PCR utilized primers designed to include the complements to the Illumina overhangs used in the first PCR, the Illumina sequencing adaptors and Illumina’s 8-bp multiplexing indices.

First-round PCRs were made to 10-µL and included final concentrations of: 1x Phusion HF buffer, 0.2 mM dNTPs, 0.2 µM of each Illumina-overhang primer, 0.2U of Phusion High-Fidelity DNA Polymerase (Thermo Fisher Scientific, USA) and one microliter of template DNA. Hotstart PCRs were initiated at 98 °C using a Eppendorf Mastercycler ep Gradient S followed by the temperature profile: 98 °C (30 sec); 30 cycles of 98 °C (10 sec), 65 °C (30 sec) and 72 °C (30 sec); and a final extension 72 °C for 5 min. Amplicons were visualized on 1% agarose gels and, after confirmation that no primer dimers were present, were diluted 10X prior to the second, indexing PCR. Here, reagent concentrations were as in PCR1, except that the concentration of the indexing primers was 0.075 µM. One microliter of the 10X diluted PCR1 product was used as template. Cycling conditions were: initial denaturation at 98 °C (10 sec); 12 cycles of 98 °C (5 sec), 55 °C (10 sec) and 72 °C (20 sec); and a final extension 72 °C for 1 min.

The products from the second PCR were multiplexed and purified using the QIAquick PCR purification kit (Qiagen, Valencia, CA) and sequenced using MiSeq Reagent v2 Nano chemistry (BPI—Biotecnologia Pesquisa e Inovação, Botucatu, Brazil).

Sequence processing and statistical analyses

We used MOTHUR v.1.36.1 (Schloss et al., 2009) to create contigs from Illumina Miseq Nano 500 paired-end reads using deltaq = 6. Subsequently, the trim.seqs command was used to filter sequences with a minimum average quality score of 25 and maxambig = 6. Primer sequences were demultiplexed and trimmed using pcr.seqs with pdiffs = 3, rdiffs = 3.

Clustering of reads into zero-radius OTUs (ZOTUs) was done using USEARCH v.11.0.667 with the module unoise3 (Edgar, 2010), which also removed chimeras based on de novo detection. For reads from the sequencing of the 17 un-pooled individuals, ZOTUs represented by fewer than 1% of the total were removed from all subsequent analyses because of the potential of contamination from parts of other species (most specimens were collected en masse in CDC traps).

We assigned ZOTUs to taxonomy using the RDP classifier within MOTHUR using classify.seqs (Wang et al., 2007) with a D2 database created from all matching GenBank entries and our own sequences derived from the Illumina sequencing of individual animals. The distance matrix used for D2 cladograms was derived using the calc_distmx command in USEARCH using default parameters.

We used linear regression to derive R² estimates that could be used to test the predictability of DNA input to read output from Illumina MiSeq sequencing. The regression equation was also used to estimate the slope of the associated regression line and y-intercept. The significance of linear comparisons was assessed using Pearson correlation. The phytools R library (Revell, 2012) was used for cladogram estimates and ggplot2 (Ginestet, 2011) for graph design.

Primer efficiency analysis

The initial evaluations for primer design and coverage used all GenBank insect sequences containing at least one of the D2 primers and a 30-bp flank (a total of 5,294 sequences for the forward priming site and 2,582 for the reverse). These included 387 culicid sequences from 233 species and 24 genera for Mozzie.D2.Uni.F. The 146 Mozzie.D2.Uni.R priming sites were form 88 species and 10 genera. Priming sites for the two oligos were found in both Anophelinae and Culicinae.

This initial analysis confirmed that Mozzie.D2.Uni.R was generally conserved across the Culicidae (Fig. 1). Although there were some mismatched base positions in GenBank culicid sequences for this priming site, we believe these were either PCR or sequence-calling errors because of the otherwise high conservation across all Arthropoda (Fig. 1).

In Mozzie.D2.Uni.F, the 3′-cytosine was nearly universal within the Culicidae, but absent in almost all other taxa, thus precluding the polymerase extension of most non-culicids during PCR. Based on the initial search parameters, Chagasia bathana was the only culicid entry that did not match the forward primer perfectly; it contained a 3′-thymine rather than cytosine in Mozzie.D2.Uni.F. This potential divergence was addressed herein by the inclusion of a Chagasia sp. individual in the pools (see below).

Based on the available GenBank sequences meeting our criteria, the only non-culicid insects with nucleotides identical in both forward (including the 3′-C) and reverse primers belonged to species from closely related Nematocera families.

In the second primer analysis, we used ecopcr to evaluate only database entries containing both priming sites. Here, we used a conservative measure of primer amplification efficacy: given that mismatches in the 3′ nucleotide of a primer will theoretically preclude any polymerase extension, and thus PCR, we characterized an oligo’s ability to amplify the markers tested by whether this nucleotide was conserved among the target taxa. The 227 Culicidae sequences meeting the ecopcr parameters were from 160 species and 21 genera, including Anophelinae and Culicinae.

As above, there were several cases in Culicidae where either the forward (five sequences) or reverse (60 sequences) primer did not match our D2 templates exactly, we believe these mismatches were likely errors for at least two reasons. First, in all of these cases, there was always another representative of the same species containing the matching template sequence. Moreover, as discussed previously, most mismatched base positions were conserved across several non-culicid taxa (Fig. 1); this was particularly evident for the reverse primer Mozzie.D2.Uni.R, which when submitted to a GenBank BLAST, is fully conserved across nearly all Arthropoda entries, and indeed across many Metazoa (including humans). Second, most of the instances of nucleotide mismatches in the primers occurred near the start of the GenBank entry. This suggests that the discordance was either a Sanger trace miscall or that the primer sequence was not trimmed prior to submission (i.e., the original oligo contained the mismatch). Indeed, our first analysis, which omitted sequences with less than 30 bp flanks, resulted in comparatively fewer mismatched oligo nucleotides.

All GenBank Culicidae D2 sequences used in the ecopcr analysis possessed the 3′-C in Mozzie.D2.Uni.F. Diptera taxa containing this nucleotide were almost exclusively Nematocera: Culicidae (17% of the 1,095 sequences with the 3′-C), Simuliidae (7%), Chironomidae (2%), and contributions of less than 5 sequences each from Blephariceridae, Ceratopogonidae, Limoniinae , Lonchopteridae, Micropezidae, Mycetophilidae , Chaoboridae, Dixidae, and Thaumaleidae. There was a sole non-Nematoceran Dipteran sequence from the Drosophilidae (Drosophila willistoni, accession number XR_002724221), although another D2 entry from this species did not possess the 3′-C.

Those non-Diptera arthropods possessing the Mozzie.D2.Uni.F 3′-C nucleotide were almost exclusively Amphipoda (Fig. 1). Although many of these taxa occur in saltwater, 63% (689 of the 1,095 sequences containing the 3′-C) included species that can also be found in freshwater: Crangonyctidae, Gammaridae, Haustoriidae, and Niphargidae.

The priming site analysis of the four other metabarcoding primer sets indicated relatively conserved annealing sites in Culicidae for 16S (Figs. S1A and S1B). In CO1, there were comparatively more mismatches due to synonymous third-nucleotide substitutions (Figs. S2A and S2B), including near the 3′-end of the oligo, a potential source of primer bias (Arnheim & Erlich, 1992). In both mitochondrial markers, however, priming sites were more conserved between the various arthropod taxa tested. The 3′ nucleotide position of all primers were generally shared extensively among all arthropods (Figs. S1 and S2).

In terms of taxonomic resolution, the D2 primer set fared well against all other markers previously used for Culicidae metabarcoding (Table 1). Among the “all arthropod” sequence dataset evaluated, only the CO1 locus used by Batovska et al. (2017) had a higher taxonomic resolution (0.910 versus 0.926, respectively). Note that this CO1 dataset contained the commonly used Folmer forward primer LCO1490 (Folmer et al., 1994), so there were relatively few Culicidae GenBank entries because of primer trimming (n = 20). Within the Culicidae-only datasets (Table 1), D2 had appreciably better species resolution than the four mitochondrial markers evaluated, including the barcode region of the CO1: 0.969 (D2) versus 0.765 (CO1; Krol et al., 2019) and 0.929 (CO1; Batovska et al., 2017).

We undertook the same evaluation of taxonomic resolution for only those species that were shared between the D2 and the other datasets in order to address biases due to different taxonomic coverage (Table S1). The primer conservation in this trimmed analysis was similar to that of the complete dataset (see ecopcr output files in Data S2). Here too, D2 outperformed the mitochondrial markers in species resolution (Table S1).

A potential explanation for D2’s apparent advantage over previously used metabarcoding primers is that the latter sets were used in ad hoc comparisons to differentiate among geographical subsets of mosquitos, rather than the amongst the entire Culicidae (for which the D2 primers were designed). Moreover, the mitochondrial amplicons were less than one half the length of the D2 sequences (Table 1), thus omitting considerable nucleotide variation.

D2 Illumina sequencing of 17 individual mosquito specimens

An average of 1,244 (SD = 479) Illumina reads were produced across the seventeen individual samples that were subsequently used to create the mock assemblages. We choose this sequencing depth to account for genomic variants in the D2 amplicon, as is common in insects (Tautz et al., 1988).

After primer trimming, the D2 amplicon length for the 17 specimens sequenced was 383 bp (SD = 9) for Culicinae and 432 bp (SD = 24) for Anophelinae. The GC percentage was identical for culicines and anophelines, 60% (SD = 2.5). The only Culicidae GenBank entry with an ambiguous binding site to the primers used herein was Chagasia bathana (accession no.: AF417831), which had thymine instead of cytosine in the a 3′-end of Mozzie.D2.Uni.F. The Chagasia sp. sample herein yielded successful PCR and sequencing results and there was no indication that it amplified with lower fidelity than other sequences (see below).

The D2 provided unique sequences for each species and also recovered the general relationships of the Culicidae topology (Fig. S3). Most species possessed multiple genomic variants within the same individual, but no sequence was shared between species and the D2 sequences of all animals were reciprocally monophyletic. The variation within individuals was due to both SNPs and indels, the latter being most common in repeating microsatellite regions. There was a total of 59 D2 sequences among the 17 specimens used. These were used to create the D2 library for the taxonomic assignments in the pooled data below. All taxonomic assignments from the Illumina sequencing of pools were based on sequences having 100% identity to those from the sequencing of these individual samples.

Results for pools

On average, 20% of sequencing output was discarded based on quality filtering parameters (Table S2). We normalized the number of sequence reads (against volume used) in all pools so that comparisons could be made between species and between pools. Following normalization, it is expected that the slope from individual pools would be 1 if all species’ DNA is equally amplifiable. Moreover, the intercept of the regression line is expected to occur near zero if there is a linear relationship between contributed DNA and reads output.

We defined a Best Estimate of the true contribution of each stock DNA by averaging the normalized read outputs for each of the 17 species from the four Fresh Pools (A–D) (Table S3). This is in effect an estimate of what percentage of reads would results if pools were created from equal initial volumes of the 17 stock solutions. The least amplicon-dense genomic DNAs were generally from older specimens or animals that had been collected in CDC traps exposed to rain.

Pearson correlation was used to test significant association between both axes in the regression calculations. There was a highly significant correlation for all comparisons (<0.01), even after adopting Bonferroni corrections for multiple tests, except in comparisons involving the 1:10 dilutions with bycatch, for which no significance was found in any of the pools (but see below).

Culicid-only pools

Table S3 and Fig. S4 present the normalized proportion of sequence reads resulting from the stock DNA averaged across the four Fresh Pools and across the four Degraded Pools. In both mixture types, the individual possessing the highest proportion of amplifiable targets was Mansonia flaveola (this specimen had been captured and its DNA preserved immediately in the field, unlike most other specimens). Although the taxon with the lowest proportion of amplified D2 in the Fresh Pools was Psorophora sp., Chagasia sp. was by far the least amplicon-dense taxon in the Degraded Pools. This, and other similar disparities vis-à-vis sample degradation may be due to differences in DNA elution solutions leading to higher degradation at room temperature.

A strong linear relationship exists in the standardized values between fresh and degraded pools (R² = 0.94, slope=1.05, y-intercept = −0.004; Fig. S4). Note however, that because of space constraints on the flow cell, the coverage for Degraded Pools was, on average, five times less than those from the Fresh Pools (Table S2 and Data S1). Thus, there is likely some variation due to sample size effects, in addition to DNA degradation.

In the five Fresh Pools, all 17 taxa were consistently recovered, including for the low-coverage replicate of Fresh Pool A (A-REP in Table S2). In these mixtures, the RDP taxonomic assignments recovered Diptera sequences in 99% of cases, the remaining balance consisted of PCR artifacts, generally of very short sequences.

The R² values of the Fresh Pools were all above 0.97. Graphical representation of the Best Estimate values against those for the four individual pools can be seen in Fig. 2.

Within the Degraded Pools, several species sometimes failed to be recovered: Chagasia sp., Culex (Culex) sp., Psorophora sp., and Mansonia humeralis (Fig. S5 and Data S1). This was not unexpected, as in all these cases of species dropout, the proportion of degraded DNA added was 2% or less of the pool volume (Data S1). Moreover, the sequencing depth for the Degraded Pools was on average five times lower than that for the Fresh Pools.

Bycatch pools

All reads from Bycatch Pools were assigned to the order Diptera (RDP), highlighting the low cross-reactivity of the designed primers. Here, comparatively fewer Culicidae reads resulted, as the primers amplified an additional five Diptera families (Table S2). By far, the main bycatch taxon represented was Chaoboridae, which has been proposed as the most phylogenetically related family to Culicidae (Kumar & Rai, 1990). The replicate samples for both the 1:1 and 1:10 mixtures of Bycatch Pool A tended to yield similar proportions of non-target reads.

Although the utmost care was taken to remove all remnant culicid tissue from the CDC bycatch prior to creating the Bycatch Pools, we believe that carryover nonetheless impacted the 1:1 and 1:10 mixtures, as it comprised 160 and 1,600 times more DNA (respectively) than the mosquitos in these pools. This is supported by the fact the more bycatch-biased pools (1:10) do not have 10 times fewer culicid reads, as would be expected (Table S2). We thus excluded from the regression analyses all specimens that had been captured sympatrically with the CDC bycatch (see Table S3).

The only 1:1 Bycatch Pool where a species was not recovered in the sequence reads was Bycatch Pool D, which could not recover the Anopheles (Kerteszia) sp. contribution. This may be a result of: low sequencing coverage (this pool yielded only 6,358 reads), low amplification density per stock DNA volume (Table S3) and because the Anopheles (Kerteszia sp.) contributed to this pool only 0.3 % of the final volume (Data S1).

The replicates for Bycatch Pools A 1:1 and A 1:10 (suffixed with “-REP” in Table S2) underwent separate PCRs and were sequenced on a different Illumina run from the original. Both replicates failed to recover any sequences for 4 of the 17 DNA samples PCRed. However, these had, on average, six times lower coverage than the original iterations (average of 2,728 reads per replicate; Table S2). The coefficient of determination of the 1:1 dilution was over 0.9. However, the R² values were very low in the 1:10 dilution because of species drop-out and the artificial enrichment for the contaminant sequences discussed above (Fig. 3, Table 2 and Data S1).

Important considerations

The adoption of modern traps in monitoring efforts represents an unprecedented ability to sample en masse and efficiently evaluate an area’s culicid biodiversity. The molecular protocol designed herein can benefit monitoring programs that utilize traps aimed specifically at hematophagous insects (e.g., Mosquito Magnet®) and can collect tens of thousands of culicids. Moreover, the primers’ amplification fidelity across primarily Culicidae also allows processing of field collections with substantial amounts of bycatch (e.g., CDC or Malaise traps, Govella et al., 2011). Notably, some of the non-culicid bycatch detected in the Illumina sequences was nonetheless also epidemiologically relevant, such as Ceratopogonidae, which comprises biting midges responsible for transmitting several human and animal diseases (Table S2). Additional work may elucidate if the D2 primer set can also be applicable in the monitoring of these dipterans.

Estimating the required sequencing depth for metabarcoding is an important consideration in terms of both time and resource efficiencies. In the Fresh Pools (which contained only Culicid DNA), each of the 17 species were consistently recovered in both pools A–D (averaging 17,000 reads per Illumina run) and the replicate pool A-REP (3,000 reads). Given the original contribution of stock DNA volume to these pools was as low as 0.2% (Data S1), within a field-collected sample (everything else being equal), the likelihood of identifying a single mosquito from a culicid-only pool of 500 is reasonable even at sequence coverage of less than 5,000 reads/pool. The MiSeq Reagent Kit v2 Nano reaction used herein yields up to one million paired-end sequences for less than US$1000, potentially allowing for multiplex sequencing of hundreds of bulk pools containing hundreds of thousands of mosquitos. In sampling scenarios with substantial bycatch, there is a need to (at least crudely) estimate the relative concentration of bycatch in order to establish the sequencing depth needed to identify all mosquito species in a collection.

Many sampling apparatuses used in culicid monitoring can be deployed during several weeks and lead to DNA and tissue loss. Here, we show that pools simulating this degradation during trap deployment will not significantly bias abundance estimates, likely due to the size of the target amplicon. Moreover, we show that even if these traps capture a substantial proportion of non-culicids, the original mosquito biodiversity will still be recoverable (with enough sequencing coverage).

Given the culicid life cycle, the primers described herein also offer substantial advantages for detection of environmental DNA (eDNA) from immatures in aquatic systems (Schneider et al., 2016). Although eDNA produces unprecedented sensitivity to mosquito presence/absence and substantially lowers sampling effort, the protocol also captures an exponentially higher ratio of bycatch DNA, given the relatively scarce aquatic culicid eDNA (Boerlijst et al., 2019). In such a context, culicid-specific primers would at worse decrease the sampling depth needed for monitoring and at best lead to the diagnosis of only the mosquito eDNA. Here, a potential drawback of the D2 primers is their co-amplification of amphipods, particularly freshwater Gammarus species (“All other arthropods” in Fig. 1). Of the GenBank sequences that possessed the 3′-C of the forward primer in the ecopcr analysis (n = 1095), 65% belonged to the amphipods (these represented the only instances where this nucleotide was present in non-Dipteran arthropods). Although this is not an issue in many scenarios, such as monitoring of mosquitos that breed in containers or temporary pools, ad hoc consideration should be given to data contamination from these animals, particularly when working with environmental DNA. Of course, this is not an issue with adult captures.

Another important caveat in ascertaining the reliability of a rDNA marker for abundance estimates is the potential for sexually dimorphic genomic D2 copies in Anophelinae. In this subfamily, the location and size of the nucleolus organizer regions (NORs) are phylogenetically inconsistent and D2 amplification may vary between homo- and heteromorphic chromosomes (Rafael, Tadei & Recco-Pimentel, 2003; Wilkins, Howell & Benedict, 2007). In Culicinae, an equal contribution of rDNA amplicons is expected, as here ribosomal DNA is concentrated in homomorphic chromosomes (Kumar & Rai, 1990).

The primer-design strategy described herein is likely adaptable to non-Culicid Metazoa taxa given the abundance of rDNA sequences on public databases. Primer design can be facilitated through available bioinformatics tools (e.g., Riaz et al., 2011; Boyer et al., 2016).