Sampling of An. stephensi
Anopheles stephensi surveys were conducted from August through November 2018 in Semera, Godey, Kebridehar, Erer-Gota, and Dire Dawa. Details about the collection sites and approach have been described previously [5]. Briefly, adult mosquitoes were collected via pyrethrum spray catch (PSC) and Centers for Disease Control and Prevention (CDC) light traps. The collected mosquitoes were sorted between culicines and Anopheles, and the latter were further sorted to distinguish An. stephensi from other Anopheles species using a standard morphological key [25] and a key prepared by Coetzee (subsequently published in 2020 [26]). Analysis of Plasmodium detection was conducted previously, and no wild-caught adult mosquitoes contained Plasmodium falciparum or P. vivax DNA [5].
Larvae and pupae of Anopheles were also collected from larval habitats including artificial water storage containers, freshwater pools, discarded tires, and plastic containers. Larvae were reared in field insectaries using water taken from larval habitats, feeding them with baker’s yeast and exposing them to sunlight during the day. Pupae were transferred into adult emergence cages and adults were morphologically identified using the same keys [25, 26]. The numbers of wild-caught adults and reared larvae are documented in Additional file 1: Table S1.
DNA extraction and molecular identification of mosquitoes
DNA was extracted from either the abdomens or whole bodies of An. stephensi mosquitoes, selected at random, using the DNeasy Blood & Tissue Kit (Qiagen, Valencia, CA, USA). PCR was performed for each individual mosquito, targeting the nuclear internal transcribed spacer 2 (ITS2) region and the mitochondrial cytochrome c oxidase subunit 1 (COI). The reagent components and final concentrations for the PCR assays were 1× Promega Hot Start Master Mix (Promega, Madison, WI, USA) and 0.5 mM for both primers, plus 1 µl of isolated DNA template. A region including the ITS2 gene was amplified via PCR using universal primers as described previously [27]. Updated annotations show that identification via these primers is based on a portion of the sequence that includes 28S in the ITS2 assay [28]. The PCR temperature protocol consisted of 95 °C for 1 min, 30 cycles of 95 °C for 30 s, 48 °C for 30 s, and 72 °C for 1 min; followed by 72 °C for 10 min. All samples successfully amplified ITS2. A region including the COI gene was amplified via PCR using universal primers as described previously [3]. The PCR temperature protocol consisted of 95 °C for 1 min, 30 cycles of 95 °C for 30 s, 48 °C for 30 s, and 72 °C for 1 min; followed by 72 °C for 10 min.
Molecular Wolbachia detection
The Wolbachia 16S rRNA-encoding gene (16S), and the Wolbachia surface protein encoding gene (wsp) were amplified to detect Wolbachia in the mosquitoes. A positive control of extracted genomic DNA from Wolbachia-infected Drosophila melanogaster (provided by The Wolbachia Project at Vanderbilt University) was used to troubleshoot PCR protocols and ensure the PCR successfully amplified Wolbachia 16S, wsp, and the genes within MLST. A negative control of no genomic DNA (gDNA) template was used to ensure no contamination of the PCR reagents. The reagent components and final concentrations for the PCR assays were 1X Promega Hot Start Master Mix (Promega, Madison, WI, USA) and 0.4 mM for each primer 16S and wsp, plus 4 µl and 2 µl of isolated DNA template for Wolbachia 16S, and wsp, respectively. For the 16S assay, a nested protocol was used.
The first set of Wolbachia 16S primers, as described in Werren and Windsor, amplified an un-nested 438-base-pair (bp) fragment [29]. The PCR cycling protocol was as follows: 95 °C for 2 min, 2 cycles of 95 °C for 1 min, 60 °C for 1 min, and 72 °C for 1 min, followed by 35 cycles of 95 °C for 30 s, 60 °C for 1 min, and 72 °C for 45 s, and lastly 72 °C for 5 min. For the nested reaction, specific internal primers were used as described by Shaw et al., which amplified a 412-bp fragment [17]. Two microliters of PCR product from the un-nested reaction was used in this reaction as template DNA. The PCR temperature protocol was as follows: 95 °C for 15 min, 35 cycles of 95 °C for 15 s, 66 °C for 25 s, and 72 °C for 30 s, followed by 72 °C for 5 min.
For the wsp assay, a 546 bp region including the gene was amplified using primers as described in Zhou et al. [30]. The PCR temperature protocol was as follows: 94 °C for 2 min, then 35 cycles of 94 °C for 1 min, 55 °C for 1 min, and 72 °C for 1 min; followed by 72 °C for 5 min. All primers used are listed in Additional file 1: Table S2. All PCR products were run on a 2% agarose gel, and the gel was visualized using a Gel Doc EZ imager (Bio-Rad Laboratories, Inc., Hercules, CA, USA) prior to Sanger sequencing by a commercial laboratory (Eurofins Genomics LLC).
Multilocus strain typing
In addition to amplifying the 16S and wsp genes, MLST was performed on each 16S-positive sample to determine whether the Wolbachia infection was present. Five loci, including ftsZ, fbpA, hcpA, coxA, and gatB, were tested, using region-specific primers. The ftsZ universal assay amplified a 775 bp region [31], the fbpA assay amplified 429 bp, the hcpA assay amplified a 444 bp region, the coxA assay amplified a 402 bp region, and the gatB assay amplified a 369 bp region [10]. The final concentrations for the PCR assays were 1X Promega Hot Start Master Mix (Promega, Madison, WI, USA), 0.4 mM for each primer, plus 2 µl of isolated DNA template for each reaction.
The PCR cycling protocol for coxA, hcpA, and gatB consisted of 94 °C for 2 min, 36 cycles of 94 °C for 30 s, 54 °C for 45 s, and 72 °C for 1 min and 30 s; followed by 72 °C for 10 min. The PCR cycling protocol for fbpA consisted of 94 °C for 2 min, 36 cycles of 94 °C for 30 s, 59 °C for 45 s, and 72 °C for 1 min and 30 s; followed by 72 °C for 10 min. Lastly, the PCR cycling protocol for ftsZ universal consisted of 95 °C for 2 min, 35 cycles of 94 °C for 30 s, 58.2 °C for 1 min, and 72 °C for 1 min; followed by 72 °C for 5 min.
A nested protocol was used for all five MLST genes. The hcpA, fbpA, and gatB nested assays were used as previously described on the PubMLST website [32]. Two microliters of PCR product from the un-nested reaction was used in this reaction as template DNA. Nested PCR cycling conditions were as follows: 95 °C for 15 min, 50 cycles of 95 °C for 15 s, 55 °C for 30 s, and 72 °C for 1 min; followed by 72 °C for 5 min. The ftsZ and coxA nested assays were used as described previously and amplified products 424 bp and 357 bp in length, respectively [33]. The nested PCR protocol for ftsZ and coxA consisted of 94 °C for 5 min, 36 cycles of 94 °C for 15 s, 55 °C for 15 s, and 72 °C for 30 s; followed by 72 °C for 10 min. All primers used are listed in Additional file 1: Table S2. All PCR products were run on a 2% agarose gel, and the gel was visualized using a Gel Doc EZ imager (Bio-Rad Laboratories, Inc., Hercules, CA, USA).
Sequence and phylogenetic analysis
Sequences were cleaned and analyzed using CodonCode (CodonCode Corporation, Centerville, MA, USA). Wolbachia 16S, Wolbachia coxA, and An. stephensi COI sequences were submitted as queries to the National Center for Biotechnology Information’s (NCBI) Basic Local Alignment Search Tool (BLAST) against the nucleotide collection under standard parameters (100 max target sequences, expect threshold 0.05, word size 28, optimized for highly similar sequences, not specific to any organism). This information was used to confirm that the amplicons produced were the Wolbachia 16S gene, Wolbachia coxA gene, or An. stephensi COI. An alignment of observed haplotypes and the nucleotide conservation in the An. stephensi Wolbachia 16S sequences were visualized on CLC Sequence Viewer 7.6 (CLC Bio Qiagen, Aarhus, Denmark).
Phylogenetic analysis was performed with the representative 16S sequences collected from this study, as well as peer-reviewed, published Wolbachia 16S sequences from NCBI GenBank in anopheline mosquitoes with an outgroup of Rickettsia japonica (Additional file 1: Table S3) [18, 19, 33,34,35,36,37,38]. Phylogenetic analysis was also performed with the coxA sequences collected from this study, as well as peer reviewed, published Wolbachia coxA sequences from NCBI GenBank in anopheline mosquitoes with an outgroup of the Wolbachia endosymbiont of Dirofilaria immitis (Additional file 1: Table S3) [18, 33, 36]. An additional phylogenetic analysis was performed with the Wolbachia 16S sequences from sub-Saharan Africa used in the first analysis, the three 16S haplotypes collected in this study, and four 16S sequences from a study in pre-print from Tamil Nadu, India (Additional file 1: Table S3) [18, 24, 33, 36]. Lastly, phylogenetic analysis was performed on the COI genes from the An. stephensi that were Wolbachia 16S-positive, with an outgroup of Anopheles maculatus (Additional file 1: Table S3). Representative sequences from each location were determined by aligning sequences and eliminating non-distinct haplotypes. Alignments were created with MAFFT version 7 [39] and trimmed to 283 bp (16S; Fig. 2), 323 bp (16S; Additional file 1: Fig. S1), 357 bp (coxA; Fig. 3), and 321 bp (COI; Fig. 4) using Mesquite version 3.61 [40]. Phylogenetic associations between all sequences were estimated using a maximum likelihood approach with RAxML [41]. The GTRGAMMA option that utilizes the general time-reversible (GTR) model of nucleotide substitution with the gamma model of the rate of heterogeneity was used. A total of 1000 runs were completed using the maximum likelihood criteria with rapid bootstrap analysis. The RAxML output was viewed in FigTree [42] with a root at the outgroup and transformed branches, and a phylogenetic tree image was made.