Ethics

This study was approved by the Papua New Guinea Institute of Medical Research Institutional Review Board (1203) and PNG Medical Research Advisory Board (12.05).

Sample collections

We collected mosquitoes from the villages of Dimer, Wasab, Kokofine, Mirap and Matukar in the Madang province of Papua New Guinea (PNG) in June and August 2012. In each village, field technicians collected mosquitoes between 1800 and 0600 using barrier screens as described by Burkot et al [32]. These screens were manually searched every 20 minutes and resting mosquitoes were captured from the screens using an aspiration device. After collection, the sex and species of each mosquito were determined by morphology as previously described [33]. All male mosquitoes and non-Anopheles mosquitoes were discarded. We visually classified each female Anopheles mosquito as non-fed, partially-fed or fully-fed by examining the size and coloration of their abdomen. We individually stored each mosquito in an Eppendorf tube containing silica gel as desiccant.

DNA isolation and molecular species identification

We extracted DNA from individual mosquitoes using a 96 well Qiagen DNeasy blood and tissue kit as previously described [34]. Mosquito species identification was determined using a PCR-based assay that evaluates species-specific polymorphisms in the ribosomal RNA internal transcribed spacer unit 2 (ITS2) [35].

In silico assessment of mammalian mt 16S rRNA primers

To test the range of mammals that should be amplified using mt 16S rRNA primers [36], we conducted an in silico analysis using the primerTree package. We also conducted in silico analyses for two other previously published primers, cytochrome oxidase I (COI) and cytochrome b (Cytb) that have been previously used for mosquito blood meal identification [37]. Since the 16S rRNA locus appeared to be the most informative for our purposes (S1 Fig), we restricted our further analyses to this locus. Briefly, we performed primer-BLAST searches using the mammalian mt 16S rRNA primer sequences against the NCBI nucleotide database using default parameters but allowing for up to three mismatches in the primer sequences. In our search, we set the maximum number of blast hits retrieved to 10,000 and retrieved the taxonomical information of each sequence retrieved. As this search can be biased by recent release of many DNA sequences from a specific taxon, we performed this search separately for each mammalian order. We then calculated how many different species were obtained from each order to calculate the total number of mammalian species likely to be amplified by this primer pair. Note that, when conducting the search without any taxonomic restrictions, these mammalian primers were also predicted to amplify amphibian and fish 16S rRNA genes.

To estimate the total number of mammalian species for which the targeted locus has been sequenced and deposited in NCBI, we randomly selected one DNA sequence from each mammalian family and used BLAST searches to identify similar DNA sequences in the NCBI nucleotide database (accessed on July 2015). We filtered out any DNA sequence from the database that did not contain the primer sequences (allowing for up to three mismatches). We then merged the results from the searches performed in each family and counted how many unique species were observed. These analyses provided us with the total number of mammalian species that should be amplified if the primers were truly universal.

We also evaluated whether the 16S rRNA primers amplified sufficiently informative DNA sequences to support rigorous species identification (i.e., whether related species could be distinguished). First, we retrieved the mammalian DNA sequence alignment from the primerTree analysis and calculated the number of nucleotide differences (including deletions) between every pair of DNA sequences using the dist.dna program of the Ape package [38]. We then calculated the average proportion of nucleotide differences between species belonging to the same mammalian order and between species belonging to different orders. Second, we used the same approach to determine, for each mammalian order, how often two different species (or genera) have the exact same DNA sequence for the targeted region of the 16S rRNA gene.

Targeted high-throughput sequencing of mammalian mt 16S rRNA genes and human mt genome hypervariable region I haplotypes

To interrogate the mammalian species composition of individual mosquito blood meals we amplified a 140 bp of the mammalian mt 16S rRNA gene using universal mammalian primers [36] modified to include a 5’-end tail complementary to the Illumina sequencing primers (S1 Table). We also attempted to amplify a subset of 192 mosquitoes with universal avian primers ([39] and S1 Fig) using a pooled approach but failed to detect any bird DNA in these samples. To identify individual differences among human blood meals, we designed PCR primers to amplify 300 bp of the human mt hypervariable region I. We first aligned 795 whole mt genomes of individuals from Oceania [40] using MAFFT version 7 [41] to evaluate the extent of sequence variation across the mt genome hypervariable region I and then designed primers positioned in conserved flanking sequences with Primer3 [42]. As described above, we added a 5’ tail to each primer for sample barcoding and high-throughput sequencing (S1 Table).

For each sample and amplicon, we performed two rounds of PCR amplification to prepare products for Illumina sequencing (Fig 1). First, we performed a locus-specific amplification (i.e., targeting either the mammalian mt 16S rRNA or the human mt hypervariable region) using the Promega GoTaq PCR kit protocol (50 μL reaction) with 1μL of genomic DNA, 0.2mM of each dNTP, 1.25 units of GoTaq DNA polymerase, 4mM of magnesium and 0.2 μM primers. PCR amplification was carried out under the following conditions: 3 minutes at 94°C followed by 30 cycles at 94°C for 45 seconds, 50°C for 45 seconds, 72°C for 30 seconds and a final elongation step at 72°C for 3 minutes. We then purified these PCR products using the QIAquick 96 PCR purification kit protocol (QIAGEN). Second, we incorporated Illumina adaptors, including unique 6-nucleotide sample identification barcode sequence through 10 additional PCR cycles, using barcoding primers complementary to the 5’-end tail incorporated during the first PCR amplification (Fig 1). For these reactions, we used the Promega GoTaq protocol as described above with 1uL of PCR product being added to each reaction. The same thermocycling conditions as described above were used but for an annealing temperature of 56°C. Predicted sizes for the mammalian mt 16S rRNA amplicons ranged from 265 to 343 bp; sizes for the human mt hypervariable region I amplicons ranged from 440 to 444 bp (amplicon sizes include Illumina sequencing primers, unique barcode sequence and Illumina adaptors, Fig 1). Finally, we pooled the barcoded amplification products from 96 individual mosquitoes and simultaneously sequenced them on an Illumina MiSeq instrument (Sequences deposited in NCBI SRA: SRP062959).

Download:

• PPT
  PowerPoint slide
• PNG
  larger image
• TIFF
  original image

Fig 1. Overview of the sequencing assay used to characterize blood meal composition of individual mosquitoes.

(i) A first PCR amplification is performed on DNA extracted from each mosquito targeting ~140 bp of the mammalian mt 16S rRNA (gray) using primers modified with a 5’-end tail complementary to the Illumina sequencing primers (red). (ii) A second PCR amplification incorporates the Illumina adaptors and a 6-nucleotide barcode unique to each mosquito at the ends of the individual blood meal PCR products. (iii) After pooling amplification products from 96 samples together, the PCR products are simultaneously sequenced on an Illumina MiSeq to (iv) generate paired-end reads (in grey) and barcode sequences (grey box). (v) Paired-end reads are then merged to provide error-corrected consensus sequence reads. The dotted black line indicates that the 6-nucleotide barcode corresponding to each read is known but is sequenced independently.

https://doi.org/10.1371/journal.pntd.0004512.g001

Bioinformatic assessment of blood meal composition from individual mosquitoes

We discarded from further analyses any read that did not carry the exact barcode and primer sequences. After recording the read origin using the barcode sequence, we removed the primer and barcode sequences to only keep the amplified DNA sequences. We discarded any resulting read smaller than 50 bp as these likely represent primer dimers. Since each amplified molecule was sequenced in both directions using paired-end reads, we merged each pair of sequencing reads using PANDAseq [43] (Fig 1) keeping, at each position, the nucleotide with the highest sequencing quality. We then analyzed 16S RNA and human mtDNA sequences separately.

Using all 43,743,363 16S rRNA sequences generated from the 442 mosquitoes, we identified all unique DNA sequences using Mothur [44] and recorded the number of reads carrying each unique DNA sequence. We removed any DNA sequence that was observed less than 10 times across all samples, as these likely resulted from sequencing errors. We compared the remaining unique DNA sequences against all DNA sequences present in the NCBI nucleotide database using blastn. For each DNA sequence, we recorded the best match(es), only considering sequences with > 90% identity over the entire sequence length. We then retrieved the taxonomic information from each best-matched sequence using the ‘get_taxonomy’ function in PrimerTree. When an amplified sequence matched multiple species equally well, we recorded all species names associated with that sequence. We then summarized the blood meal of each mosquito by calculating the proportion of reads matching each species. As a small number of reads generated could reflect low level PCR contamination or an error in the sequence barcode identification, we only analyzed mosquito samples with at least 1,000 reads (S2 Fig). For the same reason, we considered a mosquito as having fed on a single mammalian host if >90% of the sequencing reads aligned to the 16S rRNA of that species. Alternatively, if >10% of the sequencing reads aligned to a second species, we considered the mosquito to have fed on multiple mammalian hosts.

For the human mt hypervariable region, we aligned the consensus reads to the human mitochondrial reference genome sequence (NC_012920.1) using bowtie2 [45] and calculated, for each sample, the number of reads supporting each haplotype. Only haplotypes supported by more than 500 reads were considered to avoid incorporating sequencing or PCR errors (i.e., rare haplotypes that differed from an abundant haplotype by one nucleotide difference) in the analyses. Finally, we reconstructed a phylogenetic tree with all identified human mt haplotypes using MEGA version 6 [46].

In silico assessment of the amplification range and specificity of the universal mammalian 16S rRNA gene primer pairs

We first conducted extensive in silico analyses to confirm that the primer pair selected [36] could amplify DNA sequences from a wide range of mammalian orders including Primates, Rodentia (rodents), Artiodactyla (even-toed ungulates), Carnivora (carnivorans), Chiroptera (bats), Cetacea (cetaceans), Insectivora (insectivores) and Marsupials (Table 1). Overall, in silico analysis predicted that these primers should amplify 1,752 of the 1,779 mammalian species (98.5%) sequenced at this locus and present in the NCBI nucleotide database (Table 1). Besides mammals, these primers were predicted to also amplify Actinopteri (bony-fishes) and Amphibia (amphibians) (S3 Fig).

Download:

• PPT
  PowerPoint slide
• PNG
  larger image
• TIFF
  original image

Table 1. Summary of the amplification range and discriminatory power predicted for the mammalian 16S rRNA primers.

The table indicates, for each mammalian order, the number of species deposited in NCBI for the 16S rRNA genes, the number of species predicted to be amplified by the universal primers as well as the percentage of genera and species that would carry a unique sequence for this locus (enabling their rigorous identification).

https://doi.org/10.1371/journal.pntd.0004512.t001

In addition to amplifying a wide range of targets, our approach requires primers to amplify DNA sequences containing enough information to identify each species specifically. We tested this parameter by comparing the DNA sequences predicted to be amplified by this primer pair (see Methods for details). Despite the short amplified DNA sequence (~140 bp), these primers enabled rigorous differentiation of most mammalian species as illustrated by the average proportion of nucleotide differences (including deletions) between sequences of species belonging to the same or different order (S2 Table). For example, 27% of the nucleotide positions at this locus differ, on average, between one Carnivora and one Primate species and 17% of the nucleotides differ between the sequences of two Carnivora species. This high discriminating ability is also shown by the long branch lengths displayed by the phylogenetic tree reconstructed using these sequences (S1C Fig). In fact, we found that one DNA sequence amplified by these primers typically matches a single genus and, in 86% of the cases, a single species (Table 1).

Application to field-caught female Anopheles mosquito blood meals

We analyzed mosquitoes collected in five villages in the Madang Province in PNG: Dimer (n = 45), Wasab (n = 81), Kokofine (n = 83), Mirap (n = 171) and Matukar (n = 62). These mosquitoes included several species of the Anopheles punctulatus group: An. punctulatus s.s., An. koliensis, An. farauti 4 and An. farauti s.s. We characterized the blood meal composition of a total of 442 female Anopheles by amplifying the mammalian mt 16S rRNA genes from DNA extracted from these mosquitoes, pooling the PCR products of 96 samples after individual barcoding, and simultaneously sequencing the samples on an Illumina MiSeq instrument (Fig 1). We generated a total of 43,743,363 paired-end reads of 150 bp (includes added primers). For 42,198,573 DNA sequences (96.5%), we were able to collapse the overlapping paired-ends (Fig 1) and thus correct many sequencing errors. After combining the reads generated from all samples together, we identified 2,436,277 unique DNA sequences. We discarded from further analyses 2,404,684 unique DNA sequences that were carried by less than 10 reads across all samples as these likely represent DNA sequences with rare sequencing errors (accounting for a total of 4,432,784 reads or 10.1% of the total number of reads generated). We then compared the remaining 31,593 unique DNA sequences, accounting for 39,310,579 reads (89.9%), to all DNA sequences deposited in the NCBI database. 28,999 of these DNA sequences (representing 38,375,616 reads) had > 90% nucleotide identity to at least one mammalian DNA sequence present in NCBI: 18,814 unique DNA sequences best matched a single mammalian species sequence while 10,185 unique DNA sequences matched equally well to multiple mammalian species sequences (S3 Table).

Overall, we generated an average of 82,528 reads per mosquito. The number of reads generated from each mosquito varied considerably (S2 Fig) as it depends on several factors including: the amount of starting template (i.e., quantity of mammalian DNA present in the mosquito), the amplification efficiency and uneven pooling or variations in sequencing output between MiSeq runs. For further analyses, we only considered mosquito samples with more than 1,000 reads. None of the 30 extraction controls (i.e., water samples that have been processed in parallel through DNA extraction, PCR and sequencing) reached this cutoff illustrating the low level of cross-contamination or read mis-assignment due to errors in the barcode sequence (if any). Overall, we analyzed mammalian DNA from 314 blood fed mosquitoes, including 258 out of the 337 mosquitoes characterized as fully-fed (76.6%) and 56 out of the 86 mosquitoes visually-classified as partially-fed (65.1%). Only 5 out of the 19 mosquitoes visually classified as non-fed yielded mammalian 16S rRNA sequences: four yielded exclusively human 16S rRNA sequences, the last one a mix of human and pig sequences. These DNA sequences could indicate possible contamination either during field collection or in the laboratory, or detection of DNA from a previous, partially digested, blood meal. There was no statistical difference between fully-fed and partially-fed mosquitoes, however the number of sequencing reads generated for mosquitoes visually classified as fully-fed or partially-fed were significantly different from those classified as non-fed (p<0.05, Wilcoxon Rank-Sum test, S4 Fig). In total we successfully amplified and sequenced mammalian DNA from 319 Anopheles mosquitoes.

We identified 201 Anopheles mosquitoes that carried human DNA, 111 carried pig DNA, 60 carried dog DNA and 5 carried mouse DNA (Table 2; further details in S3 and S4 Tables). In addition to these expected hosts, we identified one mosquito that carried DNA from two different bat species: 7.2% of the reads matched perfectly Dobsonia moluccenis, a fruit bat commonly found in PNG, while 5.1% of the reads were most similar (94.4% identity) to another megabat species, Dobsonia praedatrix, also endemic to PNG (Table 2). These bat DNA sequences were clearly distinct (8 nucleotide differences between them) and unlikely to have been derived from sequencing errors, indicating that the mosquito fed on two different bats (S5 Fig). Additionally, in another mosquito 13% of the total reads (7,599 reads) were most similar to the common spotted cuscus (Spilocuscus maculatus, 98% similarity), a marsupial found in the forests of PNG (S6 Fig). Note that, consistent with our in silico analysis, we were not always able to identify the exact species that was fed upon. For example, we could not differentiate Canis lupus from Canis aureus (S3 Table). Overall, these finding illustrate the unbiased nature of this sequencing approach to identify host species regardless of expectations for mosquito blood meal feeding (as long as a closely related species has been sequenced).

Download:

• PPT
  PowerPoint slide
• PNG
  larger image
• TIFF
  original image

Table 2. Summary of the hosts identified in the mosquito blood meals.

For each host, the number of Anopheles mosquitoes carrying a corresponding DNA sequencing is indicated as well as the highest percent identity between the read generated and the host DNA sequence in NCBI and the average number of reads per sample carrying each DNA sequence.

https://doi.org/10.1371/journal.pntd.0004512.t002

Out of 319 mosquitoes analyzed, 52 (16.3%) showed clear evidence of having fed on more than one host species (with >10% of the reads supporting the minor host): 44 mosquitoes carried DNA from two species and eight carried DNA from three species (Fig 2).

Download:

• PPT
  PowerPoint slide
• PNG
  larger image
• TIFF
  original image

Fig 2.

Composition of the blood meals for mosquitoes collected in Mirap (A), Kokofine (B), Wasab (C), Dimer (D) and Matukar (E). Each vertical bar shows the composition of the blood meal for one mosquito: each color represents a different host species and the height of each stacked bar corresponds to the proportion of reads matching this host DNA sequence. Gray corresponds to human DNA, turquoise to pig, blue to dog, white to mouse, red to bat and orange to cuscus.

https://doi.org/10.1371/journal.pntd.0004512.g002

Within each village, we identified three major mammalian hosts—humans, dogs and pigs—accounting for 37 to 100% of each mosquito blood meal. However, the proportion of mosquitoes that fed on each host varied within and between villages (Fig 2). For example, in Mirap, only 31 of the 127 Anopheles mosquitoes (24%) fed on humans while 62 (49%) fed on pigs, 11 fed on dogs (9%) and 23 fed on two or three species (18%) including one mosquito that fed on two bat species and one mosquito that fed on a common spotted cuscus (Fig 2A). By contrast, in Kokofine, 52 out of the 62 mosquitoes fed on humans (84%) while the remaining 10 mosquitoes fed on dogs (n = 3), pigs (n = 3) or on multiple species (n = 4) (Fig 2B). The data for the three other villages are presented in Fig 2C–2E. Note that as host density information is not available for these villages, we were unable to test whether the observed differences in blood meal composition were caused by differences in mosquito feeding behavior among locations or species.

Evidence of mosquito blood meals containing multiple human hosts

Since we observed that 16.3% of the mosquitoes analyzed had fed on multiple mammalian hosts, we hypothesized that mosquitoes could also be feeding on multiple human individuals. We therefore investigated the number of different human DNA sequences present in 157 human-fed mosquitoes, using the same approach to sequence ~300 bp of the human mt hypervariable region. We generated an average of 26,721 sequencing reads of 250 bp for each sample and successfully amplified 102 of the 157 mosquitoes for the human mt hypervariable regions yielding a total of 20 different human mtDNA sequences (S7 Fig). While a single DNA sequence was amplified from 78.5% (n = 80) of the human-fed mosquitoes analyzed, 21.5% (n = 22) mosquitoes carried two distinct DNA sequences (S7 Fig). One sequence, identified in 14 of these potential mixed human blood meal, was always present at low frequency (<8% of the reads) and was actually more similar to a region of human chromosome 11 (98% similarity) than to the mitochondrial genome sequence (91%). This DNA sequence likely resulted from the amplification of the nuclear insertion of the mitochondrial sequence (numt, [47]) and was excluded from further analyses. Nine mosquitoes, belonging to two species and collected in three locations, showed presence of two human mtDNA sequences (S5 Table). For four of these mosquitoes, only one substitution (out of the 300 bp amplified) differentiated the two sequences and these could possibly be caused by a PCR error occurring at an early cycle. However, for the remaining five mosquitoes, 5–14 nucleotide substitutions differentiated the two sequences amplified and indicated that the mosquito successively fed on multiple individuals (Fig 3 and S5 Table).

Download:

• PPT
  PowerPoint slide
• PNG
  larger image
• TIFF
  original image

Fig 3. Neighbor-joining tree showing the relationships among the human mtDNA haplotypes amplified from mosquitoes.

Each symbol represents one DNA sequence amplified from one mosquito. Different shapes represent different Anopheles species (squares-An. punctulatus s.s., triangles-An. farauti 4) and is colored according to the collection site (green-Dimer, blue-Wasab, purple-Kokofine). Mixed blood meals are highlighted by boxes of the same color: for example, the two red boxes show two human mtDNA haplotypes amplified from a single An. farauti 4 mosquito collected in Kokofine.

https://doi.org/10.1371/journal.pntd.0004512.g003

DNA profiling of human maternal lineages from field collected mosquitoes

Previous studies have used microsatellites to compare the attractiveness of different individuals or group of individuals [22, 24, 25], examine the blood feeding patterns of mosquitoes [20, 53, 54] or determine the effectiveness of insecticide treated bed nets [55–58]. DNA profiling with microsatellites allows for the identification of unique genetic profiles from human individuals fed on and can be a very powerful method to differentiate DNA from unrelated individuals. However, microsatellites can only detect the simultaneous presence of multiple individual DNAs (typically two) if their proportion in one sample is relatively similar. Otherwise, the signal from the less abundant DNA is typically obscured and not distinguishable from background noise. Rigorously identifying whether a disease vector feeds on a single or multiple individuals is however essential for disease control as vectors that feed on multiple individuals are more likely to rapidly spread the disease than those that only feed on a single individual.

As an alternative to microsatellites, our approach relies on identifying unique human mitochondrial haplotypes carried by a mosquito by analyzing 300 bp of the mt hypervariable region I. We showed that at least five (out of 102 mosquitoes analyzed) carried human mitochondrial DNA sequences from more than one person. It is important to emphasize here that the number of mixed human blood meals is clearly underestimated as only maternal lineages can be detected by this approach: all offspring will carry the same DNA sequence as their mother and therefore it would not be possible to distinguish between siblings (or cousins from mothers who are sisters). However, one could, at least partially, circumvent this limitation by including additional polymorphic nuclear loci in the assay and sequence them together with the mt hypervariable region locus (and the 16S rRNA). Overall, our approach allows for a rapid evaluation of the number of maternal lineages a mosquito has fed on that can be added to the characterization of the blood meal at no additional costs, and could be used to determine if mosquitoes preferentially feed on some individuals and avoid other individuals.