Updated distribution of anopheline mosquitoes in Hokkaido, Japan, and the first evidence of Anopheles belenrae in Japan

Background In Hokkaido, northern island of Japan, at least seven cases of falciparum malaria were reported by 1951. A survey conducted at that time was unsuccessful in implicating any mosquito species as the possible vector. Although active anopheline mosquito surveillance continued until the middle of the 1980s, there is very limited information on their current status and distribution in Japan. Therefore, this study is an update on the current status and distribution of anopheline mosquitoes in Hokkaido based on a 15-year entomological surveillance between 2001 and 2015. Methods A survey of mosquitoes was conducted at 22 sites in Hokkaido, Japan, from 2001 to 2015. Adult mosquitoes were collected from cowsheds, lakesides, shrubs, and habitats ranging from open grassland to coniferous forest using a Centers for Disease Control and Prevention (CDC) miniature light trap enhanced with dry ice, aspirators, and sweeping nets. Larvae were collected from lakes, ponds, swamps, stagnant and flowing rivers, and paddy fields. All specimens were morphologically identified and subjected to polymerase chain reaction (PCR)-based sequence analysis of the internal transcribed spacer 2 ( ITS2) region of rDNA. Phylogenetic trees were reconstructed using the neighbor-joining method with the Kimura 2-parameter model on MEGA X version 10.2.2. Results A total of 46 anopheline specimens were used for the phylogenetic analysis. During the survey, a new member of the Anopheles hyrcanus group, An. belenrae, was discovered in eastern Hokkaido in 2004. Anopheles belenrae has since then been consistently found and confirmed to inhabit only this area of Japan. Four members of the An. hyrcanus group, namely An. belenrae, An. engarensis, An. lesteri, and An. sineroides, have been found in Hokkaido. The results also suggest that An. sinensis, formerly a dominant species throughout Japan, has become a rarely found species, at least currently in Hokkaido. Conclusion The updated distribution of anopheline mosquitoes in Hokkaido, Japan, showed considerable differences from that observed in previous surveys conducted from 1969 to 1984. In particular, areas where An. sinensis was previously distributed may have been greatly reduced in Hokkaido. The phylogenetic analysis revealed a novel An. hyrcanus group member identified as An. belenrae, described in South Korea in 2005. It is interesting that An. belenrae was confirmed to inhabit only eastern Hokkaido, Japan. Graphical abstract

Background Malaria cases reported in Japan reached 28,000 annually in 1945 and 1946, with over 7000 cases of vivax malaria up to the end of the 1950s. Surprisingly, at least seven cases of falciparum malaria were reported between 1947 and 1951 in Rubeshibe, Hokkaido (43.78 N, 143.61 E), located in the north of Japan. Although a survey was conducted to determine the vector mosquitoes involved at that time, no suspected species were found [1]. Endemic malaria was considered eliminated by 1960. The number of malaria cases has decreased drastically since then, with less than 80 imported cases annually in the past 10 years, and 20 imported cases in 2020 [2].
Anopheline species contain the most important malaria vector species. Among those recorded in Japan, Anopheles sinensis is the most widespread and common anopheline species. This species is considered the major vector of vivax malaria in Korea and China. Previous surveys conducted in Japan from 1970 to 1986 revealed that An. sinensis was the dominant anopheline species in Japan, including Hokkaido; An. lesteri was commonly found in Hokkaido, with only a few An. sineroides [3][4][5][6]. These surveys also found a new member of this group, An. engarensis [3][4][5]. Thus, several malaria vector species, including An. sinensis, An. engarensis, and An. lesteri, continue to inhabit Japan. Despite the need for a nationwide survey to systematically assess these species, very little information is available, mostly gathered in the 1980s. Recently, several DNA barcoding projects have been conducted on mosquitoes in Japan, and a small number of genomic data on anopheline mosquitoes were included [7][8][9]. However, these studies were not specific to malaria vector mosquitoes.
At the onset of this survey, the presence of five species of the An. hyrcanus group, namely An. sinensis, An. sineroides, An. lesteri, An. engarensis, and An. yatsushiroensis, had been confirmed in Japan. Moreover, of these five species, only An. yatsushiroensis has never been reported in Hokkaido [10][11][12][13], the region of interest in this study. Nonetheless, the highly similar morphological features of the members of this group, particularly An. engarensis and An. sinensis, make it difficult to distinguish between species morphologically. Therefore, the frequency of clasper movements in males, hybridization studies, and chromosomal studies were used to distinguish An. engarensis from the Japanese population of An. sinensis [3][4][5]. They have recently been effectively identified using polymerase chain reaction (PCR) and sequence analysis. Among the molecular markers used for mosquito taxonomy, the cytochrome oxidase subunit 1 (COI) sequences of the DNA barcoding region [14][15][16] and the internal transcribed spacer 2 (ITS2) region of rDNA are the most efficient. ITS2 in particular is very efficient in distinguishing between closely related species such as the An. maculipennis complex, An. quadrimaculatus complex, An. culicifacies complex, and An. gambiae complex [17][18][19][20]. ITS2 has also been used to address taxonomic issues in the An. hyrcanus group [21][22][23][24][25][26].
For about 20 years after the last survey in 1984 [6], very few surveys of malaria vector mosquitoes were conducted in Japan. We therefore conducted nationwide surveys between 2001 and 2015 to determine the current status and distribution of anopheline mosquitoes in Japan. In the present study, species identification and determination of genetic distances between specimens was carried out by analyzing the ITS2 region. Special attention was given to determining the distribution of An. (Anopheles) belenrae described in South Korea in 2005 [24] in Japan. Finally, we updated the information from previous surveys [3][4][5][6] on the current distribution of the anopheline mosquitoes in Hokkaido.

Mosquito sampling
Hokkaido, the study site of the present survey, is generally divided into four areas: Donan (southern Hokkaido), Doo (central Hokkaido), Doto (eastern Hokkaido), and Dohoku (northern Hokkaido). Table 1 summarizes the species and numbers of mosquitoes collected from 2001 to 2015 at 22 sites in Hokkaido classified according to the above areas. The Dohoku area was not included in this survey because it is extremely mountainous, has a harsh biological environment, and does not have very convenient transportation. Anopheline mosquitoes are the most active in Hokkaido during the summer season from July to August. This period is recognized as a very short but active season for anopheline mosquitoes. Therefore, our mosquito surveys were conducted once a year in late July or August at approximately the same location.
Adult mosquitoes were collected in cowsheds and around their habitats such as lakesides, shrubs, and open grassland to coniferous forest throughout the day using a Centers for Disease Control and Prevention (CDC) miniature light trap enhanced with dry ice [27], aspirators, and sweeping nets for approximately 3 h after sunset. Collected adult mosquitoes were frozen and transported in an icebox to the National Institute of Infectious Diseases (NIID), Tokyo, Japan. Larval mosquitoes were Keywords: Anopheles hyrcanus group mosquitoes, Hokkaido Japan, ITS2 sequence collected from paddy fields, swamps, stagnant and flowing rivers, lakes, and ponds using dippers. Larvae were transported alive to NIID and reared to adults under laboratory conditions of 25 °C, 60-70% relative humidity, and a photoperiod of 16:8 (L:D) h. Morphological identification was performed on all adult individuals using taxonomic keys [11,28]. All classified mosquito specimens were transferred individually into 1.8 ml microtubes (Eppendorf, Hamburg, Germany), and stored at −80 °C until subsequent analyses by ITS2 sequencing. The collectors always wore long-sleeved shirts, long pants, and hats, and applied repellent to the bare skin of their hands and faces to prevent mosquito bites.
In this study, eight specimens collected in domestic areas outside Hokkaido were used as a reference specimen in phylogenetic analysis. Seven of the eight specimens were from Japan, and the last was from Vietnam.  Table 2.

DNA extraction, ITS2 amplification, and sequencing
Total genomic DNA was extracted from individual samples using the REDExtract-N-Amp Tissue PCR Kit (Sigma Chemical Co., St. Louis, MO, USA) according to the manufacturer's protocol. Extracted mosquito DNA was subjected to PCR-based sequence analysis and phylogenetic analysis using primers of the ribosomal DNA ITS2 region (forward primer, 5′-TGT GAA CTG CAG GAC ACA-3′; reverse primer, 5′-TAT GCT TAA ATT CAG GGG GT-3′) [29]. Amplification conditions were as follows: initial denaturation at 95 °C for 2 min, followed by 35 cycles of 95 °C for 30 s, 55 °C for 30 s, and 72 °C for 1 min, and a 4 min final extension at 72 °C using the Veriti ™ 96-Well Thermal Cycler (Thermo Fisher Scientific Inc., Waltham, MA, USA).
All visible PCR-amplified DNA fragments were purified using the QIAquick PCR Purification Kit (QIAGEN, Venlo, Netherlands) or extracted using MonoFas (GL Sciences Inc., Tokyo, Japan) from a 2% low-melting-point agarose gel (SeaPlaque GTG Agarose, Cambrex Corp., East Rutherford, NJ, USA) after preparative gel electrophoresis and visualization with ethidium bromide. Each purified double-stranded PCR product was directly cycle-sequenced from both ends using the BigDye Terminator Cycle Sequencing FS Ready Reaction Kit v3.1 (Thermo Fisher Scientific) and the PCR primers [29]. The thermal profile used was 25 cycles of 96 °C for 10 s, 55 °C for 5 s, and 60 °C for 4 min using a thermal cycler, and an ABI PRISM 3730 Genetic Analyzer (Thermo Fisher Scientific Co.).
Sequence analysis was performed using GENETYX software version 14 (Genetyx Corp., Tokyo, Japan). Sequences of the PCR-amplified DNA fragments were then used to perform BLAST searches on the GenBank nucleic acid database of the National Center for Biotechnology Information website (http:// www. ncbi. nlm. nih. gov/ BLAST/) for species identification.

Phylogenetic analysis
Multiple alignment of the ITS2 sequences with those of related species available in the GenBank library was performed using the CLUSTALW program [30]. Phylogenetic trees were produced using the neighbor-joining (NJ) program with the Kimura 2-parameter model [31] on MEGA X version 10.2.2 [32]. The statistical significance of the resulting NJ trees was evaluated using a bootstrap test with 1000 replications.

Distribution map
A map of Hokkaido in the geodatabase (ArcGIS data collection standard pack 2014, ESRI Japan, Tokyo, Japan) was used to map the collection sites. The geographical positions of the collection sites were obtained from both previous studies [3][4][5][6] and the current study. The geographical positions obtained in this study were recorded using a geographical positioning system (GPS: GPS-MAP64, Garmin, USA) ( Table 1). These collection sites were plotted on the map using ArcGIS 10.0 (ESRI Inc., Redlands, CA, USA).

Phylogenetic analysis
A total of 248 specimens (181 adults and 67 larvae) were collected in Hokkaido between 2001 and 2015. The collected specimens were classified into four anopheline species of the An. hyrcanus group: An. belenrae, An. engarensis, An. lesteri, and An. sineroides (Table 1). Interestingly, An. sinensis was not collected from Hokkaido during our survey. Phylogenetic analysis was performed using the 485-base-pair (bp) ITS2 sequence of 38 specimens, collected from different sites (30 from Hokkaido, seven from Japanese regions outside Hokkaido, and one from Vietnam) in different years, and eight reference sequences from the GenBank database ( Table 2).
The NJ phylogenetic trees revealed five robust clades, consisting of the four species listed above and An. sinensis (Figs. 1, 2). Unfortunately, An. sinensis was not detected in Hokkaido in this study. Therefore, we analyzed An. sinensis collected from areas outside Hokkaido. The cluster of An. sinensis showed small differences (Figs. 1, 2). Nonetheless, there were no differences among the four Japanese specimens of An. sinensis (Yokohama08, Echizen379, Kaifu353, and Misawa391) and one specimen from South Korea (isolate 1).
In 2004, two larvae morphologically identified as An. sinensis were confirmed to be An. belenrae using the ITS2 sequence, marking the first record of An. belenrae in Japan (specimen Akan44). Subsequent phylogenetic analysis showed that An. belenrae was the closest related species to An. sinensis, followed by An. engarensis (Figs. 1, 2).

Intra-and interspecific ITS2 variation
The levels of nucleotide variation detected between pairs of specimens in the An. hyrcanus group are presented in Table 3. There were no genetic differences between the 10 Japanese specimens of An. belenrae (Akan44, Kushiro10, Kushiro201, Kushiro313, Kushiro418, Kawakami60, Akan712, Kushiro503, Nakagawa807, and Nakagawa26) and the Korean strain (isolate 3), with 0% pairwise divergence. This suggests that the Japanese An. belenrae and the Korean An. belenrae are the same, at least based on the ITS2 sequences. Among the An. sinensis strains, the Vietnamese strain (GLVN59) showed slight differences from the other strains, with 0.26% pairwise divergence.  An. engarensis  Regarding the clusters of An. engarensis, one specimen (Daisen55) collected in Akita Prefecture, an area outside Hokkaido, was slightly different from the nine specimens collected in Hokkaido (Abashiri02, Akashuri18, Yubari25, Yufutsu20, Kameda40, Yufutsu115, L1532, A14-22, and L1468) (Figs. 1, 2). No genetic differences were observed among the An. engarensis strains, except for the specimen Daisen55, with 4.63% pairwise divergence from the others ( Table 3). As mentioned above, the intraspecific variation between these three species was very low. The pairwise divergence was 0% for An. belenrae, 0.26% for An. sinensis, and 4.63% for An. engarensis (Table 3). Both An. sinensis and An. engarensis indicated a few differences based on the collection areas. Regarding the interspecific variation, the pairwise divergence was 13.25-13.57% between An. belenrae and An. sinensis, 13.27-14.26% between An. belenrae and An. engarensis, and 13.86-15.19% between An. sinensis and An. engarensis ( Table 3). The values among these three species indicate high levels of genetic differentiation. The NJ phylogenetic trees also showed that one specimen of An. kleini from South Korea was located closer to An. engarensis than to An. belenrae and An. sinensis (Figs. 1, 2). A detailed study of the genetic background of these species will be necessary.
In An. sineroides, no differences were found between the two specimens from Hokkaido (Fukagawa13 and Kushiro343) and the one from South Korea (SINEK02), with 0% pairwise divergence (Table 3). However, the two specimens from areas outside Hokkaido (Hida386 and Toyama80) showed a few differences from the three mentioned above, with 0.26 and 0.8% pairwise divergence (Table 3). In contrast, a large intraspecific variation was observed in An. lesteri (Figs. 1, 2). Pairwise divergence was in the range of 0.26-2.14% among the nine specimens from Hokkaido (Abashiri01, Aba-shiri17, Abashiri42, Hakodate31, Kamiiso35, Kam-eda38, Kushiro317, Kushiro505, and A2651) ( Table 3). The two Korean strains (specimen 1 and specimen 2 classified as type B and type C of An. lesteri, respectively) were quite distant from the other An. lesteri strains (Figs. 1, 2). Pairwise divergence among the 11 An. lesteri specimens was 0.26-8.96%, indicating that An. lesteri appeared to form a highly divergent population. The cluster of An. lesteri revealed low pairwise divergence, ranging from 0 to 2.14%, between the nine An. lesteri specimens from Japan and the Chinese strain of An. anthropophagus (SMMU-FK1) ( Table 3), suggesting that they may belong to the same species.

The first record of An. belenrae in Japan
Our surveys from 2001 to 2015 revealed a significant change in the distribution range of the An. hyrcanus group in Hokkaido reported until the 1980s [3][4][5][6], including the first record of An. belenrae in Japan. Two larvae collected in the Kushiro Wetland in 2004 were tentatively named An. sinensis Kushiro strain, based solely on the morphological characteristics of the emerged adults. However, phylogenetic trees constructed using ITS2 sequence revealed that this An. sinensis Kushiro strain formed a robust clade that was clearly different from the clades of An. sinensis and other Anopheles species. Interestingly, the ITS2 sequence of the Kushiro strain was not identical to that of the An. sinensis strain collected in southern Japan, outside Hokkaido, but to that of An. belenrae, a new strain reported in South Korea in 2005 [24]. The Kushiro strain could confidently be included in the An. belenrae cluster because of the absence of intraspecific divergence as mentioned above. This species was consistently found in the Kushiro Wetland after the first detection in 2004. In contrast, An. belenrae was not found outside Hokkaido in our 15-year nationwide survey. Thus, we concluded that this species is restricted to the Kushiro Wetland in Hokkaido.
The Kushiro Wetland is the largest marshland/wetland in Japan and is located in the Kushiro Plain. The Kushiro Wetland has been the focus of nature conservation efforts since before World War II, was registered as a Ramsar site in 1980, and designated as a national park in 1987. It is also famous for being the breeding ground for Japanese cranes, Grus japonensis, and many other wild birds and a protected area for natural monuments, birds, and animals; thus, land development is strictly regulated. In South Korea, An. belenrae is found in the northern part of the country near the border with North Korea [24,38,39]. In China, An. belenrae is reportedly distributed in Shandong and Liaoning Provinces in northeastern China, facing the Korean Peninsula [40]. These areas are not only geographically close to Japan, but may also have similarities in climate, vegetation, and some environmental factors with the Kushiro Wetland. However, further investigation is needed to compare the morphological characteristics of Japanese and Korean An. belenrae, and to determine the distribution of this species in locations outside Hokkaido in Japan. We hope that ecological and evolutionary factors impacting the emergence of An. belenrae will be elucidated with the development of molecular biological technology.

No information for An. sinensis from Hokkaido
The next noteworthy finding was the disappearance of An. sinensis from Hokkaido. In previous surveys, An. sinensis was generally distributed throughout Hokkaido [3][4][5][6] (Fig. 3). Although it is often found in the same larval habitat as An. lesteri, it is thought to occur more frequently in developed paddy fields and swamps [41]. In the 2000s, we did not find any An. sinensis in the habitat of An. lesteri, nor did we find any new sources or habitats (Fig. 4). It is possible that the larval habitat of An. sinensis changed drastically during the 20-year period between the previous studies [3][4][5][6] and this current study. For example, in the 1949 [1] and 1976 [6] surveys, four members of the An. hyrcanus group were detected in northeastern Hokkaido, around Rubeshibe (Fig. 3). At that time, there were paddy fields all over the district, and forestry and horse-logging were the main industries. In recent times, however, the horse-logging industry has An   Table 1 declined dramatically, and the paddy fields have been replaced with upland crops. Furthermore, neither An. sinensis nor An. sineroides was found in this area, around Ozora, during our survey (Fig. 4). It is highly likely that the changes in vegetation and industry have affected the distribution of anopheline mosquitoes.
There may be other reasons for the disappearance of An. sinensis from Hokkaido. The classification of organisms was mainly based on morphological keys until the 1990s. Although adults of An. belenrae can be separated morphologically from those of An. lesteri, An. sinensis and other species [24], in fact, it is likely that An. belenrae and An. sinensis could not be differentiated morphologically. Therefore, it should be noted that An. belenrae may have been classified as An. sinensis. The results from the ITS2 sequences in this study revealed that these two species were genetically the most closely related. In addition, the pairwise interspecific distance in mitochondrial genomes calculated by each fragment showed minor or no differences between An. sinensis, An. belenrae, and An. kleini [40]. Phylogenetic analysis of COI indicated that ancient hybridizations probably occurred among these three closely related species [42] Fig. 2 Phylogenetic relationships among members of the Anopheles hyrcanus group. Neighbor-joining (NJ) phylogenetic tree was constructed based on partial ITS2 region nucleotide sequences. A distance matrix was calculated using the Kimura 2parameter evolutionary model and the tree was constructed using the NJ approach in MEGA X version 10.2.2. This phylogenetic tree is a rootless and radiation tree with An. pullus (AY186792) as an outgroup sequence to determine more closely related species. The scale bar indicates the proportion of sites changing along each branch. Bootstrap values are not indicated in this figure. The species names in this group are in bold. The members in the same species are shaded. Abbreviations of specimen and sequence accession numbers of specimens used in this study are listed in Table 1  Table 3 Percentage of pairwise divergence among five members of Anopheles hyrcanus group mosquitoes in Japan

An. belenrae
The number of base substitutions per site from between sequences are shown. Analyses were conducted using the Kimura 2-parameter model [25]. This analysis involved 46 nucleotide sequences. All ambiguous positions were removed for each sequence pair (pairwise deletion option). There were a total of 485 positions in the final dataset Evolutionary analyses were conducted in MEGA X [26] Specimens collected in this study are labeled with the species name, specimen name, and country listed in Table 1. Specimens found in GenBank are labeled with species name, specimen name, and country Values denoted in italic and bold indicate ITS2 intraspecific differences   28 49.18 this problem, we tried to extract DNA and decipher the nucleotide sequence from age-old, dried specimens previously classified as An. sinensis collected in Hokkaido [3][4][5][6]. However, no new information could be obtained from these specimens. We hope that techniques for genetic analysis using age-old specimens will soon be available.

Anopheles engarensis in Japan
Anopheles engarensis is also a species whose distributional range in Hokkaido has decreased. This species, first described in Engaro-cho (northeastern Hokkaido) in 1977 [3], was also found in Monbetsu, Kushiro, and Obihiro until 1984 [6], suggesting a wide distribution in Hokkaido [3][4][5][6] (Fig. 3). However, our surveillance found this species to be restricted to western and southern Hokkaido (Fig. 4). In addition, the species was collected in northern Tohoku, Akita Prefecture, suggesting a southward shift presumably due to environmental changes, including the climate of larval habitats. In terms of classification, An. engarensis was recognized as a new species in the An. hyrcanus group only after its chromosomal structure was determined to be different from that of An. sinensis [4]. This was because of the high morphological similarity between the two species. Indeed, the only distinguishing feature was the unique number of clasper movements of An. engarensis males during artificial mating, a common method for laboratory maintenance of anopheline mosquitoes [5].
In general, ITS2 is known to have high interspecific and low intraspecific variability; however, extensive intraspecific variations have been reported in anopheline mosquitoes. For instance, ITS2 intraspecific variations ranged from 0.2 to 19.0% for the Latin American anophelines [43]. In the An. hyrcanus group, the average intraspecific distance was 0.3%, but no intraspecific variations were observed in An. belenrae [42]. These results suggest that the ITS2 spacer is a good marker for differentiating between members of the An. hyrcanus group. In this study, there were no intraspecific variations in the An. belenrae, An. engarensis, and An. sineroides strains from Hokkaido. However, there was significant intraspecific variation among the nine An. engarensis strains from Hokkaido and the specimen Daisen55 from Akita This map was created from references [3][4][5][6]. This map shows mosquito species collected from 29 sites. Black, blue, red, and orange circles indicate collection sites of An. sinensis, An. engarensis, An. sineroides, and An. lesteri, respectively Prefecture. The genetic distance of 4.7% was considerably greater than the 0.22% intraspecific variation in the An. sinensis strains from Vietnam, South Korea, and Japan. We inferred that the Daisen55 An. engarensis strain was not introduced from Hokkaido but inhabited the Tohoku region independently. On the other hand, in species groups consisting of recently diverged members, such as the An. gambiae complex, the interspecific differences in ITS2 were reported to be minor, ranging from 0.4 to 1.6% [20]. It is possible that An. engarensis is a recently diverged lineage.

Anopheles lesteri in Hokkaido
In Rubeshibe Hokkaido, at least seven cases of falciparum malaria were recorded between 1946 and 1947. A survey conducted to determine the vector mosquitoes involved in the transmission was unsuccessful, although An. sinensis and An. sineroides were collected [1]. During the falciparum malaria epidemic in the vicinity of Guangdong City, China, around 1942, the transmission was inferred to have involved An. lesteri and not An. sinensis.
This inference was based on results from field investigations and subsequent infection experiments with Plasmodium falciparum [44]. Based on this inference, it was suggested but never confirmed that An. lesteri may have been involved in the outbreak of falciparum malaria in Rubeshibe, Hokkaido. In terms of distribution, An. lesteri, which was initially thought to be restricted to western islands of Japan such as the Kyushu Island [44], was also found in various areas of Honshu in the mainland of Japan [12], Hokkaido [45], Okinawa Island, and Yaeyama Islands [46]. The present survey confirmed that An. lesteri is still widely distributed in Hokkaido (Figs. 3, 4). At the start of our survey in 2001, we observed female mosquitoes collected in Ozora, Hokkaido, to have an intense affinity for human blood. These female mosquitoes were therefore considered, and subsequently confirmed, to be An. lesteri based on the reported high anthropophilic nature of An. lesteri relative to An. sinensis and other members of the An. hyrcanus group [12,47]. Consequently, we expected to easily collect An. lesteri in subsequent surveys in Hokkaido. In our study, the ITS2 intraspecific divergence in An. lesteri was 0-9.44%. These values suggest that An. lesteri is a highly divergent species when An. lesteri type B (specimen 1) and type C (specimen 2) from South Korea are included in this species. Since the ITS2 distance of this species varies even within Hokkaido, there is a possibility that An. lesteri includes crypto-species. In defining this species, it is necessary to analyze both the COI barcoding region and the ITS2 region. Moreover, a large number of specimens collected outside Hokkaido will be necessary. In a previous study, a short interspecific divergence of 7.2% was observed between An. kleini and An. engarensis [42]. We obtained similar ITS2 divergence of 6.47% and 8.67% between An. kleini and our An. engarensis specimens from Hokkaido and An. kleini and An. engarensis specimen Daisen55, respectively. Although these results may provide validation that An. kleini is a synonym of An. engarensis, further analysis is required. We also presented evidence that An. anthropophagus and An. lesteri were conspecific, based on the ITS2 divergence between them. Our results based on interspecific comparisons of ITS2 divergence may also support previous reports that An. belenrae and An. sinensis are genetically distinct [24,25], and An. anthropophagus is a conspecific species of An. lesteri [34,48].

Changes in the distribution of the An. hyrcanus group in Hokkaido
In previous studies, An. sinensis, An. lesteri, and An. yatsushiroensis in Japan were reported to preferentially invade livestock barns and houses [12]. These species are considered to be more endophilic. However, there is very limited information on the behavior of other members of the An. hyrcanus group. Therefore, in order to collect as many mosquitoes as possible in this study, we attempted to collect both adults and larvae of mosquitoes from all areas using several methods regardless of mosquito behavior. In addition, the surveys were conducted in late July and August, when the mosquitoes are the most active in Hokkaido. In the previous studies [3][4][5][6], most specimens were obtained in August as well. The selection of the present survey sites was also based on these precedents. As we did not plan to conduct regular annual surveys, some negative factors such as the small number of surveys and collection sites may be considered as limitations. However, at least at the time of our survey, An. sinensis may not have been distributed in Hokkaido or may no longer have been present in sufficient numbers to be collected. In fact, the An. sinensis specimens used in this study were collected easily in areas south of Hokkaido, such as Kanagawa, Aomori, Fukui, and Tokushima Prefectures, using the same method as that followed in Hokkaido. Therefore, it is unlikely that the sampling methodology or the timing of the survey influenced our results. Similarly, our results indicate a significant change in the distribution range of other members of the An. hyrcanus group in Hokkaido from that reported until the 1980s [3][4][5][6]. However, although multiple surveys were conducted in the Doo and Doto areas, no further information is available for the Donan area, except for the results of a survey conducted in 2003. The possibility that the distribution of An. sinensis was reconfirmed after 2003 cannot be ruled out. To address these questions, continued mosquito surveys are needed in the future.

Conclusions
ITS2 sequence divergence revealed the current distribution of the An. hyrcanus group mosquitoes in Hokkaido, demonstrating great differences from surveys conducted between 1969 and 1984. In particular, the area inhabited by An. sinensis has greatly diminished, and the newly discovered An. belenrae was confirmed to inhabit only eastern Hokkaido. In summary, this study showed that Hokkaido harbored four members of the An. hyrcanus group, namely An. engarensis, An. belenrae, An. sineroides, and An. lesteri.
Our research has revealed that two anopheline species reported as malaria vectors, An. lesteri and An. belenrae, are present in Hokkaido today. Although the malaria vector capacity of the Japanese strain of An. belenrae has not yet been evaluated, the Korean strain is considered to be a vector or potential vector of Plasmodium vivax [24,49]. Fortunately, all recently reported cases of malaria in Japan have been imported. However, emergence of potential autochthonous malaria epidemics should always be of concern, because multiple malaria vector species still remain in Japan, as confirmed in this study.