Single nucleotide polymorphism analysis of the ITS2 region of two sympatric malaria mosquito species in Sweden: Anopheles daciae and Anopheles messeae

Four species of the Anopheles maculipennis complex have previously been recorded in Sweden. A recent addition to the complex is Anopheles daciae, which is considered to be closely related to, but distinct from Anopheles messeae. The original designation of An. daciae was based on five genetic differences (161, 165, 167, 362 and 382) in the second internal transcribed spacer (ITS) 2 of the ribosomal RNA. Further studies have shown that only two nucleotide differences (362 and 382) robustly separate the species. Thirty‐three An. maculipennis complex mosquitoes were collected in the province of Uppland, Sweden. All were An. daciae but showed double peaks for three variable positions (161, 165 and 167). When cloned, the intra‐individual nucleotide variation was almost exclusively fixed with either TTC or AAT, originally diagnostic for An. messae and An. daciae, respectively. To further investigate the intra‐individual variation, nine An. daciae and 11 An. messeae were collected in southern Sweden and their ITS2 fragments were amplified and sequenced using Illumina MiSeq sequencing (Illumina, Inc., San Diego, CA, USA). For the diagnostic nucleotide 382 no intra‐individual variation could be detected. However, although each An. daciae specimen carried several ITS2 sequence variants for the four other nucleotides, there was no intra‐individual variation in the An. messeae specimens.


Introduction
Human malaria as a result of Plasmodium vivax was prevalent in Sweden during previous centuries until about 1860 (Jaenson, 1983). The main malaria vectors in Sweden are considered to have been Anopheles messeae and Anopheles atroparvus (Jaenson, 1983;Jaenson et al., 1986). Based on isoenzyme profiles and classical egg morphology, the geographical distributions of the member species of the Anopheles maculipennis complex were investigated in the early 1980s. It was shown that the complex in Sweden consists of four species: An. atroparvus, Anopheles beklemishevi, An. maculipennis s.s. and An. messeae (Jaenson et al., 1986). More recently a new species Anopheles Correspondence: Olle Terenius, Department of Cellular and Molecular Biology, Uppsala University, PO. Box 596, SE-751 23 Uppsala, Sweden. Tel.: +46 72 561 49 59; E-mail: olle.terenius@icm.uu.se daciae, which is closely related to An. messeae, was described from Romania (Nicolescu et al., 2004) and has subsequently been shown to be distributed in several European countries but, until now, it has not been reported from Scandinavia. The species status of An. daciae was based on egg morphology and five nucleotide differences compared with An. messeae in the DNA sequence of the rRNA internal transcribed spacer 2 (ITS2). These differences, so called single-nucleotide polymorphisms (SNPs), were confirmed in a population from England and Wales where all specimens identified as An. daciae possessed the five-nucleotide SNPs in agreement with An. daciae from Romania (Linton et al., 2005). By contrast with these findings, investigations in Italy (Di Luca et al., 2004) and in Russia (Bezzhonova & Goryacheva, 2008) showed that there is a variation in the ITS2, which contradicted the designation of An. daciae as a new species based on these five SNPs. A study by Danabalan et al. (2014) concluded that it is only one of these originally five SNPs that is diagnostic. More recent studies have shown that An. daciae also occurs in Germany (Kronefeld et al., 2012;Kampen et al., 2016), Poland (Rydzanicz et al., 2017) and Serbia (Kavran et al., 2018). To investigate the potential presence of An. daciae in Sweden and analyse the robustness in species designation based on the nucleotide highlighted by Danabalan et al., 2014, adult females of the An. maculipennis complex were collected from two geographically separate regions in Sweden several years apart. The mosquitoes were then sequenced for SNP analysis.

Mosquito collection
In June 2011 and in August 2012, 28 non-bloodfed An. maculipennis s.l. mosquito females were collected from two horse stables (and the attics above the stables) at Linné's Sävja and Kungshamn's Farm situated approximately 5 and 8 km, respectively, south of the Uppsala city centre, province of Uppland, south-central Sweden (see Supporting information, Fig. S1). Captured mosquitoes were kept at −20 ∘ C until examined. In June 2012, another five mosquitoes were collected from a cellar in Lindsbro, a location north of Uppsala, province of Uppland and kept at −20 ∘ C until examined. In December 2017, 20 mosquitoes, morphologically determined as An. maculipennis s.l., were collected in stables at two locations, Högerödja and Sjöliden, province of Småland (see Supporting information, Table S1), south Sweden and kept at −20 ∘ C until examined.

DNA amplification
For species diagnosis of the 28 specimens collected south of Uppsala city, one leg of each mosquito was used directly for a polymerase chain reaction (PCR) in 25-μL reactions with Illustra PuReTaq Ready-To-Go PCR Beads (GE Healthcare, Uppsala, Sweden) using primers for the second internal transcribed spacer of the ribosomal RNA gene (ITS2) for four sibling species of the An. maculipennis complex: An. maculipennis s.s., An. atroparvus, An. messeae and An. beklemishevi (see Supporting information, Table S2). The PCR program was: 95 ∘ C for 3 min, followed by 35 cycles of [95 ∘ C for 30 s, 60 ∘ C for 30 s, 72 ∘ C for 1 min], and then a final extension step at 72 ∘ C for 7 min. Amplification products were sequenced at Macrogen (Seoul, South Korea). For analysis of individual sequences, amplified fragments were cloned into the TOPO TA Cloning ® Kit (Invitrogen, Carlsbad, CA, U.S.A.) and sequenced using Sanger sequencing technology at Macrogen. Sequences are deposited in GenBank with accession numbers MN840043-MN840073.
For the 20 specimens collected in Småland, as well as the five specimens from Lindsbro, province of Uppland, amplification of phylogenetically determining sequences was performed as described previously (Lilja et al., 2018) using DNA extracted from legs from mosquitoes homogenized in 30 μL of Prep-Man Ultra (Applied Biosystems, Foster City, CA, U.S.A>). In brief, the homogenized solution was incubated at 100 ∘ C for 10 min, centrifuged at 15 000 g for 3 min and then 20 μL of the supernatant was collected into a new tube and used as template for PCR reactions. PCR was run with Amplitaq DNA polymerase (Invitrogen) using primers for ITS2, for the NADH dehydrogenase subunit 4 gene, for the NADH dehydrogenase subunit 5 gene, and for a fragment from the Hunchback gene (see Supporting information, Table S2). Amplification products were sequenced at Macrogen. Sequences were assembled using bioedit (http://www.mbio.ncsu.edu/BioEdit/bioedit.html). Sequences are deposited in GenBank with accession numbers MN780511-MN780580.

MiSeq
To study intra-individual variation in ITS2 sequences, tagged versions of ITS2 primer 5.8S (TCGTCGGCAGCGTCAGAT-GTGTATAAGAGACAGTGTGAACTGCAGGACACATG) and ITS2 28S (GTCTCGTGGGCTCGGAGATGTGTATAA-GAGACAGATGCTTAAATTTAGGGGGTA) were used to amplify the ITS2 fragments. The products were then indexed using Illumina Nextera XT indexes, pooled and run on Illumina MiSeq using MiSeq Reagent Kit v2 2 × 250 bp (Illumina, Inc., San Diego, CA, USA). Using blast (https://blast.ncbi.nlm.nih.gov), the resulting reads were compared with a local database of ITS2 sequences with variation only on the five nucleotides that vary between An. messeae and Anopheles daciae and sorted in accordance with their sequence variant.

Results
Mosquitoes from five different locations and three different years were analysed (for an analysis scheme, see the Supporting information, Fig. S2). From the 2011 collection, 19 An. maculipennis s.l. female mosquitoes were examined by PCR for species designation. Of these, three were An. maculipennis s.s., whereas the other 16 were designated An. messeae/daciae based on the size of the ITS2 PCR fragment. The PCR products of the 16 An. messeae/daciae mosquitoes were sent for sequencing and all specimens were found to have SNPs according to An. daciae in positions 362 and 382 (AC). However, in positions 161, 165 and 167, the sequenced PCR reactions had double peaks for each position previously used for discriminating An. messeae/An. daciae (A/T, A/T, T/C). The PCR products from eight of the An. daciae mosquitoes from the 2011 collection were cloned to determine the reason for the double peaks in the chromatograms. From each mosquito, two clones were sequenced. The sequences retrieved from each mosquito were found to be identical to either of the positions previously used for discriminating between An. messeae and An. daciae (TTC or AAT, respectively) but, commonly, both genotypes were found in the same individual (Table 1).  Table S1).
In 2012, not only another three An. messeae/daciae adult female mosquitoes were collected from the same area as in 2011 (Kungshamn's Farm), but also nine An. messeae/daciae adult female mosquitoes were obtained from another location (Linné's Sävja) 4.5 km away from the location in 2011. Although there are no data on the potential flight distance in An. messeae/daciae, the distance is anticipated to be possible to cover for An. messeae/daciae mosquitoes (at least over a couple of generations) and thus An. messeae/daciae mosquitoes from the two locations can be regarded as belonging to the same population. The PCR products of ITS2 again showed double peaks for positions 161, 165 and 167. The data were thus similar for both 2011 and 2012 indicating that the presence of both the 'An. messeae genotype' and the 'An. daciae genotype' had been stable in the population. To determine whether this existence of dual genotypes in the same individuals was restricted to the same local population, five An. messeae/daciae mosquitoes were also collected from Lindsbro, approximately 40 km from the locations just outside Uppsala. The results were identical to those from the locations south of Uppsala with double peaks for positions 161, 165 and 167 and An. daciae genotype for positions 362 and 382.
Anopheles maculipennis s.l. mosquitoes were also collected at two locations in Småland, Sweden and identified by ITS2 sequencing. Using the two diagnostic nucleotides 362 and 382, nine An. daciae and 11 An. messeae specimens were identified. For positions 161, 165 and 167 double peaks were seen in the nine An. daciae specimens. To further investigate the intra-individual variation in ITS2, next generation sequencing was used to sequence the variation in ITS2 in each individual. Reads were sorted depending on the sequence variant on the five polymorphic nucleotides 161, 165, 167, 362 and 382 using BLAST against a database of ITS2 sequences varying only on these residues. For each specimen, 4800-117,000 reads were analysed with a mean of 54,000 reads. Our analysis shows that there is intra-individual variation in An. daciae specimens but not in An. messeae on the three bases 161, 165 and 167 (Fig. 1). The analysis further shows that this variation is contained to three variants, AAT, ATC and TTC. Of these, TTC is fixed in An. messeae and it is also the most common one in An. daciae.
To further analyse potential differences in gene sequences between An. messeae and An. daciae, three additional genes were sequenced in the 20 specimens from Småland: ND4 and ND5 in the mitochondrial DNA and Hb. In all of these genes, there are polymorphic sites often with several specimens of each allelic type (ND4, 15 polymorphic sites; ND5, 19 polymorphic sites; and Hb 4, polymorphic sites). Phylogenetic analysis showed that none of these polymorphic sites could be used for separation of the two species (see Supporting information, Fig. S3). Each gene groups the specimens differently and so there is no alternative grouping of the specimens suggested by these gene sequences.

Discussion
During the recent era of molecular analysis of mosquitoes, three new species in the An. maculipennis complex have been described partly based on nucleotide differences: An. persiensis Linton, Sedaghat & Harbach (Sedaghat et al., 2003); An. daciae (Nicolescu et al., 2004) and An. artemievi (Gordeev et al., 2005). Of these, An. daciae was described as having five nucleotide differences in the ITS2 region compared with An. messeae. For three of these nucleotides, the present study shows that both haplotypes may exist in the same individual and thus cannot be used to distinguish An. messeae from An. daciae. However, the presence of double peaks in these positions is only seen in An. daciae and not in An. messeae. Danabalan et al. (2014) also showed that intra-species variation occurs on four of the nucleotides and concluded that only the fifth position of the originally five discriminatory SNPs were useful as a discriminatory SNP for An. messeae/An. daciae. Bezzhonova & Goryacheva (2008) presented evidence for intra-individual variation in the ITS2 sequence of An. messeae mosquitoes covering the same ITS2 region. They found three specimens from which they cloned and sequenced several variants on all five nucleotides differing between An. daciae and An. messeae, which might suggest that hybrids could occur. In the Sanger sequencing in the present study, no specimen revealed intermediate levels of the two diagnostic nucleotides, which would have been expected for a hybrid An. messeae/daciae specimen, and the MiSeq sequencing showed no intra-individual variation on the fifth position. Thus, even in areas where both An. daciae and An. messeae are hibernating in the very same sites, no hybrid mosquito was found, bearing in mind that the sample size was limited. A larger collection of An. messeae / An. daciae would be needed to reveal whether such hybridization may occur in Sweden. No evidence for variation was found in other nucleotide positions except for four of the five original positions that discriminated between An. messeae and An. daciae. Because the variation in other nucleotide positions reported by Bezzhonova & Goryacheva (2008) was restricted to single specimens, they need to be confirmed before they can regarded as true SNP sites.  , and in Lindsbro, Uppland (21-25). A) Variation in nucleotides 161, 165 and 167 is seen in An. daciae (specimen 1, 3-6, 8, 14, 16, 20-25) but not in An. messeae (specimen 2, 7, 9-13, 15, 17-19 Variation in the ITS2 region of An. messeae was also reported by Di Luca et al. (2004) with five different haplotypes corresponding to different regions of Europe. In their study, haplotype 1 (H1) from Northern Italy and haplotype 5 (H5) from Somerset Levels, U.K., are those represented in the Uppsala populations, apart from two specimens that differ from all the five haplotypes presented by Di Luca et al. (2004). These two mosquitoes from 2012 had clones with clear peaks with the sequence ATC. They were both collected from the same locality (Kungshamn's Farm). At this locality, no mosquito was found with the ATC variant in 2011.
The markers ND4, ND5 and Hb were also sequenced to search for other sequence variations that could differentiate the species. However, despite having polymorphic sites in the three markers, they do not separate the species phylogenetically but, instead, vary across the two species. This suggests that the separation of the two species An. daciae and An. messae is a recent event. A consequence of the present study is that previous observations of An. messae in Sweden should be regarded as An. messae s.l. until their ITS2 sequences have been determined. In further studies, methods aiming to look for population differences at a broader scale, such as double digest restriction associated DNA sequencing could be used to further investigate the relationship between An. daciae and An. messeae.