Anopheles maculipennis Complex in The Netherlands: First Record of Anopheles daciae (Diptera: Culicidae)

: Despite their past importance as vectors of indigenous malaria, the species composition and spatial distribution of the members of the Anopheles maculipennis complex have been studied to a limited extent in the Netherlands. Therefore, this investigation focuses on the distribution of the members of this complex in the Netherlands, including Anopheles daciae , which has recently been found in countries bordering the Netherlands. In the framework of a national mosquito surveillance between 2010 and 2021, a total of 541 specimens of An. maculipennis s.l. were analyzed from 161 locations covering the entire territory. In addition, 89 specimens were analyzed from overwintering sites during the winter of 2020/2021. All individual mosquitoes were identiﬁed to species-level using Sanger sequencing of the ribosomal internal transcribed spacer 2. To characterize the habitat of An. maculipennis s.l. in the Netherlands, land cover use data was extracted in a 1 km buffer area around each ﬁnding location. For populations collected in summers between 2010 and 2021, the most frequent species was An. messeae , present in 88.19% of the locations, followed by An. maculipennis s.s. (11.80%), An. atroparvus (3.72%) and An. daciae (3.72%). Anopheles daciae was found in the southern inland areas of the country. Furthermore, An. messeae and An. daciae occurred in sympatry at overwintering sites. This study provides relevant information on the occurrence of species of the Anopheles maculipennis complex in the Netherlands, contributing to a better estimation of the risk of mosquito-borne disease in the country.


Introduction
In the Netherlands, the last published checklist of mosquitoes (Diptera: Culicidae) included 35 indigenous species [1].This list is not static and was recently updated with a new indigenous species, Culiseta longiareolata [2].Four members of the Anopheles maculipennis complex (= s.l.) are reported as present in the published checklist [1], some of which are capable of carrying pathogens of medical importance, including malaria: Anopheles atroparvus van Thiel, 1927, An. messeae Falleroni, 1926, An. maculipennis sensu stricto Meigen, 1818, and An.melanoon Hackett, 1934.However, this latter species is not considered to occur in the Netherlands by [3], and An.maculipennis s.s. is considered uncommon in the country [4].Anopheles atroparvus is reported as the only malaria vector (Plasmodium vivax and Plasmodium malariae) in the coastal areas of the Netherlands.Experimentally, An. atroparvus can also be an effective host and capable of transmitting Plasmodium ovale to humans [5].A comparison of the anopheline species composition between 1935 and 1999 showed a prevalence shift from An. atroparvus to An. messeae in the Delta of the Rivers Rhine and Meuse, coinciding with the disappearance of indigenous malaria [6].In a study of overwintering mosquitoes in several farms in the Netherlands, An. messeae individuals were more frequently found and An.atroparvus was less common [4].A decline of An. atroparvus over the 20th century was recorded in multiple European countries and was assumed to be linked to major ecological changes, such as drainage practices, surface water pollution, loss of suitable resting sites for hibernation, etc. [6][7][8].
Despite their past importance as vectors of indigenous malaria and their potential role in the transmission of imported tropical malaria, leading to the reappearance of autochthonous malaria cases in Europe [9], the species composition and spatial distribution of the members of An. maculipennis s.l. in the Netherlands has been poorly studied.Individuals of An. maculipennis s.l. were found at 144 sampling sites during a nationwide inventory of indigenous mosquitoes, involving natural, rural and urban habitats [10].However, this nationwide inventory did not include DNA-based species identifications to distinguish the members of the complex.Yet, the nuclear ribosomal internal transcribed spacer 2 (ITS2) flanked by portions of the conserved 5.8S and 28S rDNA is useful in this respect [11][12][13].Since the members of the An.maculipennis complex are difficult to discriminate by morphological characteristics, their identification needs to be verified by ITS2 sequencing.
Anopheles daciae Linton, Nicolescu & Harbach, 2004, is a recently described species of the An.maculipennis complex that is distributed throughout continental Europe [14] and has been found over the past years in the countries bordering the Netherlands, including Germany [15,16], the United Kingdom [11], and Belgium [17].Nevertheless, An. daciae has not yet been reported in the Netherlands.Therefore, the aim of this study was 1) to investigate the distribution of An. maculipennis s.l.members present in the Netherlands by applying ITS2 sequencing and 2) to find evidence of the presence of An. daciae in the Netherlands, where it is expected to occur, given its presence in neighboring countries.

Sampling
Mosquito specimens for the study were collected from different surveys.Most of the specimens (n = 531) were collected between the months of May and September using Mosquito Magnet Liberty Plus traps (WoodstreamTM Co., Lititz, PA, USA) using octenol, in the framework of the National Mosquito Survey [10].This included 145 specimens collected in 2011, 74 in 2012, 146 in 2013, 72 in 2014, 88 in 2015, 3 in 2016, 12 in 2017 and one in 2021.Additional specimens were collected during Exotic Mosquito Surveys [18] using a variety of sampling methods such as BG-Sentinel traps (n = 4) or BG-Mosquitaire traps (n = 5) (Biogents AG, Regensburg, Germany) both using BG-Sweetscent, and larval sampling using aquarium nets (n = 1).In total, 541 specimens (540 adults and 1 larva) of An. maculipennis s.l. were collected at 161 locations covering the entire territory of the Netherlands.
In addition to these 541 specimens, a total of 89 An. maculipennis s.l.specimens were collected in February (n = 65) and March (n = 24) 2021 from six bunkers of the New Dutch Waterline, located in the municipality of West-Betuwe, the Netherlands.These bunkers are well-known overwintering sites for several mosquito species [19].
All specimens were transported to the laboratory and were morphologically identified to the Anopheles maculipennis complex level using the key of Becker et al. [20].After identification, specimens were placed in sterile vials and kept frozen at -20 • C until further processing.

DNA-Based Species Identification
Individual DNA was extracted from a leg or a part of abdomen using the NucleoSpin ® Tissue DNA extraction kit (Macherey-Nagel, Düren, Germany), following the manufacturer's protocols, with an elution volume of 70 µL.For the mosquitoes collected from overwintering sites in February and March 2021, DNA was extracted following the ammoniumhydroxideprotocol as described by [21].The ITS2 fragment was amplified using the primers of [16], with thermal cycling conditions, PCR reactions and purification for sequencing following [17], except for the mosquitoes collected in overwintering sites in February and March 2021.For these latter, ITS2 fragment were amplified using MyTaq ® HS Red mix (Bioline, UK) using the primers as described in [16] and the following thermal cycling conditions: 1 min at 95 • C, followed by 35 cycles of 95 • C for 15 s, 53 • C for 15 s, and 72 • C for 10 s.Forward and reverse strands were assembled and corrected with Geneious ® Prime v.2019.2.3 (Biomatters Ltd., Auckland, New Zealand), after which consensus sequences were generated and trimmed to remove the primers and low-quality ends.
ITS2 consensus sequences were used as queries to search for most similar sequences in GenBank (NCBI, National Centre for Biotechnology), using the Basic Local Alignment Search Tool (https://blast.ncbi.nlm.nih.gov/Blast.cgi,accessed on 15 June 2022).To discriminate between An. daciae and An.messeae, aligned consensus sequences were visually checked for the presence of the five species-specific diagnostic sites [12].Species-assigned consensus sequences were then aligned together with sequences of all other Anopheles species occurring in the Netherlands, namely An. algeriensis Theobald, 1903, An. claviger (Meigen, 1804), An. plumbeus Stephens, 1828 [22], and with outgroup sequences (namely An. funestus sensu stricto Giles, 1900, and An.minimus Theobald, 1901), using ClustalW in Geneious ® Prime v.2019.2.3.We also included ITS2 sequences of An. melanoon Hackett, 1934, which is part of An. maculipennis s.l. and was reported to occur in the Netherlands in the past [22] but was not identified in previous reports [3].Conspecific identical sequences were removed from the database to retain unique ITS2 sequences in the final alignment.Using the web application FindModel (http://hiv.lanl.gov/content/sequence/findmodel/findmodel.html,accessed on 15 June 2022), the Kimura 2-parameter was identified as the best evolution model describing our data [23,24].A rooted maximum likelihood tree (ML) was constructed using MEGA X v.10.0.5 (Kumar et al., 2018), with branch support assessed by 1000 bootstrap replicates.Average interspecific K2P distances were calculated with the R v.3.6.1 package Spider v3.6.2 [25,26].

Habitat Characterization
To characterize the habitat of An. maculipennis s.l. in the Netherlands, a 1 km buffer area was created around each finding location (excluding overwintering locations).Using this pre-defined buffer zone in ArcGIS v.10.7. 1 [27] and the 2018 version of the raster file of the Corine Land Cover [28], the areas covered by the five main Land Cover Classes were extracted: artificial or urban areas, agricultural areas, forest and seminatural areas, wetlands, and water bodies.The expected number of specimens per Land Cover Class was calculated using the following formula: (total number of specimens per species × percentage Land Cover in the Netherlands)/100.Habitat association per species was verified by comparing observed number of specimens per species per Land Cover Class with expected number of specimens per species per Land Cover Class using a Fisher's exact test.All statistical analyses were conducted in RStudio v1.4.1717 [25].

Results
The ITS2 fragment was scored in the 541 specimens collected between May and September.Of these, 496 specimens were assigned as An.messeae, 25 as An.maculipennis s.s., 11 as An.atroparvus, and nine as An.daciae.Of the 89 mosquitoes collected from the six bunkers, 82 specimens were identified as An.messeae (92.13%), and seven were identified as An.daciae (8.99%).
No ITS2 sequences were shared between the four species of An. maculipennis s.l.collected in the Netherlands, with each of the four species involving one unique species-specific haplotype (Figure 1).Average interspecific K2P distances ranged from 0.687 to 8.258% (Table 1).Double peaks at two of the five supposedly diagnostic sites discriminating An. messeae from An. daciae were observed in three ITS2 sequences of An. daciae, namely position 218 (A/T) (n = 2) and 220 (C/T) (n = 1) (site numbering following [12]).Such ambiguities were not recorded in the 578 An. messeae ITS2 sequences (214 (T), 218 (T), 220 (C), 416 (G), and 436 (G)).Excluding the sampled overwintering sites (i.e., the six bunkers), the most frequent species captured was An. messeae, being present in 88.19% of the locations, followed by An. maculipennis s.s.(11.80%),An. atroparvus (3.72%) and An.daciae (3.72%) (Figure 2, Table S1).Anopheles messeae was found in a total of 142 locations, being found in sympatry with An. maculipennis s.s. at four locations and with An. atroparvus at five locations.Anopheles daciae was found in sympatry with An. maculipennis s.s. at three out of the six locations where it was identified (Figure 2).Using ITS2, the present investigation provides the first solid evidence of the occurrence of An. daciae in the Netherlands.The 16 identified specimens were captured only in 2011, 2012, 2013 and 2021, in the southern part of the country (Figure 2).At overwintering collection sites, the seven An. daciae were collected in February from five out of six bunkers; no An.daciae were identified in March.Anopheles maculipennis s.s. and An.atroparvus were not found at the overwintering sites.Except for An.maculipennis s.s., preferred Land Cover Class for An.messeae, An. daciae and An.atroparvus are areas with predominant agricultural use (Figure 3).Anopheles maculipennis s.s. was found in Land Cover Classes artificial habitats (e.g., industrial or residential areas) and agricultural areas in similar proportions.Anopheles atroparvus was not found in Land Cover Classes forests or seminatural areas.Anopheles atroparvus was collected in the north of the country, at six locations near the coast, where it is more probable to encounter mixing seawater and fresh water.A significant difference between expected and observed distributions per Land Cover Class was observed for An.messeae (p < 0.001), indicating a significantly higher occurrence at artificial Land Cover Classes than expected.For An. maculipennis s.s.(p = 0.164), An. daciae (p = 1) and An.atroparvus (p = 1), no significant differences between expected and observed distributions were found (Figure 3).

Discussion and Conclusions
In the Netherlands, four species of Anopheles maculipennis s.l. were identified: An. messeae, An. maculipennis s.s., An. atroparvus and An.daciae.Our study presents the first report of An. daciae in the Netherlands using ITS2 sequencing.
Each identified species of the complex involved a single ITS2 haplotype, with smallest average interspecific K2P distances between An. daciae and An.messeae (0.7%).However, ambiguous sites at some of the species-diagnostic positions were observed [17,[29][30][31].Double peaks in chromatograms can result from slight differences among ITS2 copies (heterozygosity) and their regular occurrence in specimens from different countries, surveys and years, suggests that only two sites are diagnostic between An. daciae and An.messeae, namely positions 416 (A/G) and 436 (C/G).For the remainder, phylogenetic DNA sequence analyses of ND4, ND5, COI and Hunchback gene fragments provide no support for the distinction of An. daciae and An.messeae [17,31], while other taxonomically diagnostic features between these nominal species are still poorly investigated (e.g., hybrid incompatibility, morphology, ecology, cytotaxonomy, etc.).Therefore, it has been proposed to regard An. daciae as a species inquirenda (i.e., a species of doubtful identity [32]) [17,33].
Our study shows that An. messeae is the species with the widest distribution in the Netherlands compared to any other species of An. maculipennis s.l.Similarly to a study in Germany [14], An. messeae was the most frequent species in the analyzed samples.In the Netherlands, the species was most commonly found in agricultural areas.Important and ecologically relevant features of the Dutch agricultural landscape are the drainage ditches, which are an important refuge for biodiversity in agricultural landscapes [34] and represent a preferred breeding site for Anopheles mosquitoes.However, when comparing the distribution of An. messeae with the overall distribution of Land Cover Classes in the Netherlands, the species seems to be associated with artificial Land Cover Classes, more than we would expect based on the land cover distribution of the entire country.Interestingly, the prevalence of An. messeae observed in the present study is comparable to that reported by [4], which focused on overwintering mosquitoes.In the latter study, An. maculipennis s.s. was not detected, and the presence of An. daciae was not investigated by ITS2 sequence data.Anopheles atroparvus, a common species breeding in brackish waters, was found at three sites in the coastal areas of the Netherlands during the study of [4].
Anopheles maculipennis s.s. is also widely distributed in the Netherlands, covering large areas from North to South, but it was not found near the coastal areas.Our data show that this species occurs in similar proportions in both urban and agricultural Land Cover Classes.In our study, no preference of An. maculipennis s.s.towards a specific class was identified.However, similar to An. atroparvus and An.daciae, the collected number of An. maculipennis s.s.specimens was too low for adequate analyses.As such, further research is needed to identify habitat preferences.In Belgium, An. maculipennis s.s.appears to be the most frequent and widespread species of the complex [17].However, this observation was based on a survey of artificial breeding sites.Therefore, the higher prevalence of An. maculipennis s.s. is not surprising, since this species seems better adapted to artificial habitats compared to other species of An. maculipennis s.l.[35][36][37].In the present study, most of the specimens were collected using adult traps at randomly generated locations across the country, including urban, rural and natural Land Cover Classes, thus preventing a sampling bias [10].The occurrence of An. maculipennis s.s. in urban areas in our study (>40% of the sampling sites) seems to corroborate its preference for man-made habitats.
Except for An.atroparvus [6], the potential role of An. messeae, An. maculipennis s.s. and An.daciae in the historical transmission of malaria in the Netherlands is unknown.For An. daciae, the present highlighted species distribution does not fit with the historical areas where the malaria parasite occurred [6].The main vector of malaria was An. atroparvus [6], a species dependent on brackish water.Nowadays, in comparison with the other members of the species complex, the species distribution of An. atroparvus indicates the presence of the species in scarce locations nearby coastal areas.
This study also shows that An. messeae and An.daciae live in sympatry during winter.It remains unclear, however, how the species within the An.maculipennis complex differ in overwintering strategies.Earlier studies have shown that An. atroparvus enters diapause in early winter, but occasionally continues its blood-feeding behaviour to maintain fat reserves, in contrast to An. messeae, which remains inactive throughout winter [38,39].To our knowledge, this is the first study to find An.daciae overwintering in artificial shelters together with other mosquito species that are in diapause [19].Furthermore, An. daciae was only found in the southern inland areas and occurring in areas with Land Cover Classes associated with agricultural activities.
This study provides accurate and unbiased information on the occurrence of the species of the Anopheles maculipennis complex in the Netherlands and shows that An. messeae is the most frequent species in the country during summers and can occur in sympatry with An. daciae in winters.Unbiased occurrence data are needed to develop mosquito species distribution models that can contribute to a better estimation of the risk of mosquito-borne diseases in the country.In addition, while the ITS2 gene fragment is an adequate tool for the identification of the species of the Anopheles maculipennis complex, exploring the species' whole genomes will further help elucidate the phylogenetic relationships between the complex members and support their taxonomic status.

Figure 1 .
Figure 1.Condensed ITS2 ML-tree (K2P model) of five members of Anopheles maculipennis s.l.(An.maculipennis s.s.; An. messeae; An. atroparvus; An. daciae; An. melanoon), four of which were collected in the Netherlands in the present study, including An. plumbeus, An. claviger and An.algeriensis occurring in the Netherlands, and An.funestus sensu stricto and An.minimus as outgroups (GenBank accession numbers: KP298399, KP298400, OK570292, OK570315).Duplicate sequences per species were excluded, with the ITS2 databases for An.messeae and An.daciae including a few sequences displaying ambiguous sites.Numbers at nodes are bootstrap support values.* species reported in the Netherlands.

Table 1 .
Descriptive statistics of the genetic diversity of ITS2 within Anopheles maculipennis s.l. in the Netherlands, including the average interspecific K2P distances among sequences (excluding conspecific identical sequences).: number of haplotypes.N P : number of polymorphic nucleotide sites.* ambiguities recorded at two of the five species-diagnostic sites.Numbers in bold are specimens collected at overwintering sites (2020/2021).
n: sample size.N H