Multiplex vs. singleplex assay for the simultaneous identification of the three components of avian malaria vector-borne disease by DNA metabarcoding

Accurate detection and identification of vector-host-parasite systems are key to understanding their evolutionary dynamics and to design effective disease prevention strategies. Traditionally, microscopical and serological techniques were employed to analyse arthropod blood meals for host/parasite detection, but these were limited in taxonomic resolution and only to pre-selected taxa. In recent years, molecular techniques have emerged as a promising alternative, offering enhanced resolution and taxonomic range. While singleplex PCR assays were used at first to identify host, vector and parasite components in separate reactions, today multiple primer pairs can be combined in a single reaction, i.e. multiplex, offering substantial time and cost savings. Nonetheless, despite the potential benefits of multiplex PCR, studies quantifying its efficacy compared to singleplex reactions are scarce. In this study, we used partially digested mosquito blood meals within an avian malaria framework to jointly identify the host, vector and parasite using multiplex DNA metabarcoding, and to compare it with separate singleplex PCRs. We aimed to compare the detection probabilities and taxonomic assignments between both approaches. We found both to have similar performances in terms of detection for the host and the vector, but singleplex performed better than multiplex for the parasite component. We suggest adjusting the relative concentrations of the PCR primers used in the multiplex assay can increase the efficiency in detecting all the components of the studied multi-species system. Overall, the results show that DNA metabarcoding is an effective approach that could be applied to any vector-borne interaction involving blood-feeding arthropods. Our insights will not only refine laboratory procedures, but also enhance research efforts and medical diagnosis of vector-borne diseases.

Methods
Sample collection

36 wild gravid female mosquitoes were captured between June and September 2019 in two distinct locations in Switzerland (Dorigny (46°31’N, 6°34’E) and Monods (46°34’N, 6°24’E) forests, Canton Vaud). In addition, 15 Culex pipiens mosquitoes reared in the laboratory and fed on experimentally Plasmodium-infected captive birds (Serinus canaria, Fringillidae, and Passer domesticus, Passeridae) were included as positive controls. The abdomens of all 51 mosquitoes were dissected under a binocular microscope using sterilised Vannas-Tübingen Spring Scissor (n°15003-08, Fine Science Tools) and stored at -80°C until DNA extraction.

DNA extraction

Prior to DNA extraction, mosquito abdomens were grinded with homogenisation beads on a tissue homogeniser (Bertin Technologies, Montigny-le-Bretonneux, France). DNA extraction was done using a DNeasy Blood & Tissue Kit (Qiagen, Hilden, Germany) following the manufacturers' protocol but with the elution step modified as follows: (i) 75 µL of elution buffer and 15 minutes of incubation at room temperature, followed by centrifugation, and (ii) 100 µL of elution buffer and one minute of incubation followed by centrifugation.

Primer description and preliminary tests

Three primer pairs targeting different mitochondrial DNA (mtDNA) genes were used to identify the blood meal origin (i.e. host), the mosquito species (i.e. vector) and the parasite based on the DNA extracted from the 51 mosquito abdomens.

First, a metabarcode was amplified to identify the blood source using Aves02 (Taberlet et al., 2018), targeting the 12S rDNA gene in birds. As Aves02 primers were only tested using in silico PCR so far, preliminary PCRs with these primers were done on tissue DNA extracts (liver, muscle, blood) of 21 avian species from a wide range of orders (Table S1). The presence of DNA bands at the expected molecular size was verified on an 2% agarose gel. The results of this test showed that the targeted Aves02 metabarcode was not amplified for the Eurasian blue tit (Cyanistes caeruleus), which is a common species in the study sites. Therefore, based on Passeriformes sequences complementary to Aves02 primers found in GenBank, two new versions of the forward Aves02 primer targeting the blue tit Cyanistes caeruleus (Aves02_Cyan_F) and the great tit Parus major (Aves02_Parus_F; also common in the study sites) were designed (Table S2) to better cover the avian taxonomic panel. These two modified versions were mixed, in a proportion of 25% each, with the standard forward Aves02 primer for the final PCRs (the mixture of these three forward and the original reverse PCR primers is thereafter noted as Aves02_mix). In case of failed amplifications with Aves02_mix primers, samples were tested using a more generalist vertebrate primer (Vert01,(Riaz et al., 2011; Taberlet et al., 2018)) to check if the host was a non-avian vertebrate or to confirm the negative detection of the host in the blood. Second, to identify the mosquito vector, we used a primer pair targeting the 16S rDNA gene in Culicidae (Culi01;(Schneider et al., 2016)). Third, primers targeting the genus Plasmodium and amplifying a 37-44bp portion of the cytochrome b gene were newly designed for this study (Plas01) using theOBITools package (Boyer et al., 2016). In brief: we used ecoPrimers software(Ficetola et al., 2010; Riaz et al., 2011) to find the best suitable primer for this purpose based on all complete mtDNA genomes available in GenBank. It was validated with a in silico PCR using ecoPCR software(Ficetola et al., 2010) for the targeted taxonomic group, using the EMBL database release 142 (January 2020). As Plas01 primers were also only tested using in silico PCR so far, preliminary PCRs were performed on DNA extracts of 18 Plasmodium mtDNA lineages (Table S1). The presence of DNA bands at the expected molecular size was verified on an agarose gel.

All blood meal samples were amplified in triplicate with the three primer pairs in both singleplex and multiplex assays.

Singleplex DNA metabarcoding

To test for the presence of PCR inhibitors, preliminary PCRs were performed with each primer pair on 8 randomly selected samples with different dilutions, and the amplicons band intensities checked on agarose gels. Based on these results, DNA extracts from wild blood-fed mosquitoes were diluted 2-fold and those from laboratory blood-fed mosquitoes 2000-fold before being amplified with Aves02_mix. DNA extracts were diluted 10-fold before being amplified with Culi01, and no dilution was applied for the amplification with Plas01.

The PCR reagents and their final concentrations were the following: AmpliTaq Gold 360 Master Mix 1x, tagged forward and reverse primers at 0.5 µM, BSA at 0.16 mg/mL. The final volume per well was 20 µl including 2 µL of template DNA. The PCR thermal profile started with denaturation at 95°C for 10 minutes, followed by 40 cycles with Aves02_mix, 45 cycles with Plas01 and 50 cycles with Culi01. Each cycle was composed of 30 seconds at 95°C, 30 seconds at 56°C, 55°C and 60°C for Aves02_mix, for Plas01 and for Culi01 primers, respectively, and one minute at 72°C, before a final elongation at 72°C for seven minutes.

Multiplex DNA metabarcoding

The three primer pairs used for singleplex were mixed to perform a single PCR reaction using the following final concentrations: Aves02_mix 0.4 µM, Plas01 0.4 µM, Culi01 0.15 µM. Final concentrations of the other PCR reagents were the same as in the singleplex assay, as well as final and template DNA volumes. We reduced the concentration of Culi01 because the preliminary tests revealed an overamplification of mosquito sequences with equal primer concentration (data not shown). For the multiplex PCR reactions, all DNA extracts were diluted 10-fold. The PCR conditions were as follows: initial denaturation at 95°C for 10 minutes, followed by 45 cycles with 30 seconds at 95°C, 30 seconds at 55°C and 1 min at 72°C; final elongation at 72°C for seven minutes.

Purification, library preparation and sequencing

For each assay, extraction negative, PCR negative and positive controls (Table S3) as well as blanks were included in each PCR plate (for details on plate layout, see Taberlet et al., 2018). Before pooling amplicons per plate, amplification success was verified on an agarose gel for a subset of samples.

Amplicon pools were purified using the MinElute PCR Purification Kit (Qiagen, Hilden, Germany) and quantified using a Qubit 2.0 Fluorometer (Life Technology Corporation, USA). Library preparation was performed using the TruSeq DNA PCR-Free Library Prep Kit (Illumina, San Diego, CA, USA) starting at the repair ends and library size selection step, with an adjusted beads ratio of 1.8 to remove small fragments. After adapter ligation, libraries were validated on a fragment analyser (Advanced Analytical Technologies, USA). According to the results, one or two, depending on the libraries, SPRIselect bead purifications were performed again to better select the fragments of interest. Final libraries were quantified by qPCR, normalised and pooled before 150 paired-end sequencing on the Illumina Miniseq Sequencing System with a Mid-Output Kit (Illumina, San Diego, CA, USA).

Bioinformatics

Processing of raw sequences was conducted separately for each library using the OBITools package(Boyer et al., 2016). Initially, reads were assembled with a minimum quality score of 40. Sequences were then assigned to samples based on unique tags combinations. OTUs with less than 100 reads per library were discarded as well as those not fitting the range of metabarcode lengths. Afterwards, pairwise dissimilarities between OTUs were computed and lesser abundant ones with single nucleotide dissimilarity were clustered into the most abundant ones. We used the sumaclust algorithm(Mercier C, 2013) to further refine the resulting clusters based on a sequence similarity of 97 %. It uses the same clustering algorithm as UCLUST(Prasad, D.V., 2015) and identifies erroneous sequences produced during amplification and sequencing. Remaining sequences were assigned to taxa using a reference database. We built a database using the ecoPCR for Aves02_mix, Culi01 and Plas01 by running an in silico PCR based on all the potentially amplifiable sequences available in the EMBL database (European Molecular Biology Laboratory). Since not all of these bird species were present in our database, missing sequences for frequent species occurring in the study area were generated by DNA Sanger-sequencing and added manually to the database (Aves02_mix Database, Supplementary data). Using these databases, we retained OTUs with a similarity match above 90% for Aves02_mix, 95% for Culi01 and 90% for Plas01. After taxonomic assignment using the abovementioned databases, all remaining OTUs were then double-checked using BLASTn on the NCBI database.

Further data cleaning and statistical analyses were performed in R 4.02. First, we used the metabaR package(Zinger et al., 2021) to assess the coverage, tag-jump rate and contamination of remaining PCR replicates. PCR replicates with too small reads count were also discarded. Removal of tag-leaked sequences was done independently for each library. This approach allowed us to discard single OTUs instead of whole PCR replicates. Remaining PCR replicates were grouped by sample and primer pair, and the mean number of sequences and mean relative read abundance (RRA) was calculated for each. We also re-calculated RRA to obtain the amount of reads per sample and primer pair relative to the total amount of sequences retrieved combining the three primers. We then visualised the differences between singleplex and multiplex in terms of primer pair, sample and taxonomy.