# Phylogeography andWolbachia Infections Reveal Postglacial Recolonization Routes of the Parthenogenetic Plant Louse Cacopsylla myrtilli (W. Wagner 1947), (Hemiptera, Psylloidea)

Authors: Christina Nokkala; Valentina G. Kuznetsova; Nazar A. Shapoval; Seppo Nokkala
Journal: Journal of Zoological Systematics and Evolutionary Research
 (2022)

---

## Abstract

To reveal the phylogeographic pattern of the parthenogenetic psyllidCacopsylla myrtilli (W. Wagner 1947) (Hemiptera, Psylloidea), we sequenced a 638 bp fragment of the mitochondrial COI gene from 962 individuals. These insects originated from 46 sampling sites, which cover a significant part of the northern Palearctic distribution range of the species. The sequence data revealed 40 haplotypes, with three main (H1, H2, and H3) and 37 derived ones. The main haplotypes H1 or H2 or both were present at all sampling sites. The star-like shape of the haplotype networks indicated recent population expansion. In most cases, the derived haplotypes were specific for each country, suggesting that the main haplotypes H1 and H2 are of refugial origin, and the derived haplotypes have emerged after the postglacial recolonization process. Based on the haplotype sequences, we suggest H3 to represent the ancestral haplotype from which H1 and H2 have evolved. We suggest that the main haplotype H3 together with its derived haplotypes represents bisexual C. myrtilli, which shows a limited distribution on both sides of the border between Finland and Russia in northern Fennoscandia. The genetic diversity was the highest in Sjoa in southern Norway and also high in the White Sea region in northwest Russia. Higher diversity in Sjoa was attributed to both earlier recolonizations compared to that of the White Sea region and the absence of Wolbachia infection. We suggest that these sites were colonized from different Pleistocene refugia, i.e., from western and eastern refugia, respectively. From the White Sea region, recolonization continued eastwards to Ural Mountains and westwards to Finland and further north to Kola Peninsula. From northern Finland, recolonization continued to Finnmark, Norway, and further to Sweden and finally reached a secondary contact zone with colonizers from Norway in Central Sweden. The Caucasus and Siberian/Manchurian refugial regions have played an important role in the origin of C. myrtilli populations in Siberia and the Russian Far East.


---

## Body

## 1. Introduction

Jumping plant lice or Psylloidea, a hemipteran superfamily, consist of mostly monophagous or oligophagous species feeding on plant sap, being so-called phloem feeders [1–4]. The life cycle of psyllids is very tightly coupled with the hostplant. Mating occurs on the plant, and females lay their eggs on the leaves and shoots of the plant. In addition, nymphs develop on the plant, and adults feed on the plant. In previous studies, all psyllids were considered to reproduce bisexually [5], though suggestions of parthenogenetic reproduction were made for species where no males had been encountered [6–9]. Cacopsylla myrtilli (W. Wagner, 1947) is a triploid species with an obligatory apomictic parthenogenetic reproduction [10]. In some populations, rare nonfunctional males with abnormal meiosis, as well as rare diploid females, occur with a low frequency (<5-10%) [11, 12]. However, in a geographically restricted area, in northern Finland and northwest Russia, males show regular meiosis leading to normal sperm production [11] and are thus functional. Cytogenetic analysis, together with the mitochondrial COI markers, demonstrates that triploid parthenogenetic females produce rare males and diploid females in every generation [12]. Recently, apomictic parthenogenetic reproduction has also been found in the triploid C. ledi (Flor, 1861) [13] and in the pentaploid C. borealis Nokkala & Nokkala, 2019 [14]. The apomictic parthenogenesis is a strictly clonal reproduction system. Consequently, mutations are the only way to generate genetic variation and hence new clones.The boreoalpine speciesC. myrtilli lives on two closely allied host plants, Vaccinium myrtillus L. and V. uliginosum L. [15]. The present-day distribution of the species covers the entire subarctic and alpine areas in the northern hemisphere. In the western Palearctic, the species forms extremely dense populations in northernmost Fennoscandia and at high altitudes but is much sparser at altitudes near sea level. It is very rare in southwest Finland (pers. observ.), has not been recorded from southern Sweden or Denmark [15], and appears in northern Italy, Bulgaria, Czech Republic, Germany, Switzerland, Austria, northern and eastern Russia, and Carpathians at high altitudes [15]. This kind of patched distribution at high latitudes and altitudes is thought to reflect adaptation to harsh environments (see [16]), and C. myrtilli seems to fit the model very well.It is well known that parthenogenetic species are most abundant within marginal areas, such as deserts, islands, high altitudes, high latitudes, and areas that were formerly covered by ice or exposed to extreme aridity, a phenomenon known as geographical parthenogenesis (reviewed by [17]). It is speculated that if there were both a bisexual and a parthenogenetic form of the species in the refugia during the glacial period, it was the parthenogenetic one that first recolonized areas freed from ice cover. After ice retreat, the postglacial colonization of northern Europe is considered to follow two alternative routes. The western route led directly from the south when Sweden was physically connected to the European continent via Denmark, while the eastern route ran across Russian Karelia and northern Finland, while large areas of present-day Sweden and Finland were still covered with ice or water [18–20]. The weight of the ice caused the Earth’s crust to be pressed down, and when the ice melted, the crust began to rise slowly. The influence of the former ice cover is evident even today, as the uplifting of land is continuing along the coasts of the Gulf of Bothnia, which carried the thickest ice load [21, 22].According to the traditional concept, the southern refugial peninsulas, Iberia, Italy, and the Balkans, together with possible additional refugia near the Caucasus and the Caspian Sea, have been the major regions where species have survived during the last glacial period that ended ca. 10 000 BP. At that time, most of the European continent was covered with ice [18–20]. Depending on various biogeographical, abiotic, and biotic factors, the recolonization of northern Europe has followed various routes for different species but with some general features. According to Hewitt [19], most species recolonized northern Europe from the Pleistocene refugia in Iberia or Balkans as the Alps slowed down recolonization from Italy. The expansion patterns produced hybrid zones when lineages originating from different refugia met.More recently, the importance of so-called cryptic or extra-Mediterranean refugia has gained considerable support (reviewed by [16, 23–25]). It has been suggested that smaller refugial regions have existed in several locations in central Europe and even as far north as the Norwegian coast [26, 27]. Stewart and Lister [23] have summarized the existing evidence for extra-Mediterranean refugia and presented how combined recolonization from southern peninsular refugia with recolonization from smaller northern refugia may enhance greater complexity of the recolonization process and the observed genetic structure of present-day populations. However, common to all views of postglacial recolonization is the presence of multiple routes, originating from different refugia. Although general patterns of recolonization are known, the details of the recolonization of a particular species require to be individually determined.In the present study, we sampledC. myrtilli populations throughout its distribution range in Fennoscandia and Northern Russia as well as from the eastern part of the range in the Sea of Okhotsk region and analyzed the distribution of cytochrome c oxidase subunit I (COI) haplotypes. We aimed to describe the phylogeographic population structure of this species to bring new insights into the routes of recolonization of an apomictic parthenogenetic species after the last glacial period. We also addressed the question whether the present-day distribution of populations originates from a single distinct refugium or several refugia and sampled three central European populations from Czech Republic and Bulgaria. In addition, we aimed to recognize the original refugial haplotypes and possible new local haplotypes that emerged after postglacial recolonization.


## 2. Materials and Methods

### 2.1. Samples

The speciesC. myrtilli was identified using the key and description by Ossiannilsson [15]. The specimens were collected on V. myrtillus and V. uliginosum in 46 locations, 6 in Norway, 6 in Sweden, 14 in Finland, 17 in Russia, 2 in the Czech Republic, and 1 in Bulgaria during 2007-2020 (Table 1 and Figure 1) and stored immediately in 96% ethanol in the field. Partial COI gene was sequenced from 964 specimens. Because of heteroplasmy (Supplementary Figure S1), two of these were not included in the study (Table 1).Table 1
Number ofCacopsylla myrtilli individuals sequenced in samples from Norway, Sweden, Finland, Russia, Czech Republic, and Bulgaria. Altitudes below 200 m a.s.l. not shown.
Sampling locationNLatitudeLongitudeAltitude a.s.l./mDateNorway1Finnmark, Šuoššjávri3269°22′11″24°18′20″Aug 10, 20122Sjoa, Rudihøe3661°46′27″9°17′16″1000Aug 1, 20093Sjoa, Kvernbrusætrin2461°42′27″9°19′25″950Aug 1, 20094Sjoa, Stålane2461°41′15″9°14′27″1000Aug 1, 20095Sjoa, Kringlothaugen2261°43′06″9°22′40″700Aug 1, 20096Sjoa, Rindhovda 10803061°43′05″9°05′12″1080Aug 16, 2010Sweden7Abisko, Björkliden Fjället2268°24′32″18°39′55″500Aug 3, 20098Abisko, Lapporten 5702868°19′14″18°51′05″570Aug 10, 20129Soppero2168°00′39″21°39′25″450Aug 10, 201210Jokkmokk2566°35′36″19°49′20″300Aug 8, 201211Sorsele2665°29′07″17°33′36″430Aug 8, 201212Storuman2163°08′23″17°15′33″400Aug 8, 2012Finland13Utsjoki, Ailigas4569°53′51″27°03′32″320July 30, 201514Utsjoki, Hietala2469°51′06″27°00′34″Aug 11, 201215Utsjoki, Karigasniemi1769°25′14″25°58′32″450Aug 16, 201716Kilpisjärvi, Pikku-Malla1769°08′50″20°44′20″620Aug 6, 201617Sodankylä, Puisuvanto1867°46′52″26°46′09″200Aug 3, 201818Sodankylä1267°22′49″26°49′35″Aug 12, 201219Salla, Tuntsa1467°39′40″29°45′53″July 21, 201020Kuusamo1866°13′57″29°08′48″Aug 12, 201221Liminka2864°43′53″25°23′04″Aug 4, 200922Paltamo5064°33′28″27°43′41″July 16, 200923Juuka1263°15′04″29°11′12″Aug 12, 201224Suonenjoki662°40′56″27°03′48″July 17, 200925Tohmajärvi1862°23′12″30°19′12″Aug 12, 201226Forssa460°58′37″23°31′00″Aug 13, 2012Russia27Kildin Island2769°19′11″34°20′55″Aug 14, 202028Murmansk, North Nagornoje2068°55′03″33°04′01″July 16, 201229Murmansk, Kola1868°52′35″33°01′11″July 28, 200730Murmansk, Laplandskii Reserve3268°04′12″32°16′12″Aug 1, 201931Apatity2067°34′29″33°23′27″Aug 6, 201232Vorkuta3067°30′64°02′Aug 4, 201333Island Malyi Gorelyi (White Sea)1666°17′58″33°37′18″Aug 4, 201234Island Bolshoy Gorelyi (White Sea)366°17′96″33°37′30″Aug 4, 201235White Sea, mainland3366°17′40″33°39′06″Aug 18, 201636Island Sredniy (White Sea)4166°17′17″33°38′01″July 13, 201237Kartesh (White Sea)366°20′40″33°38′23″Aug 6, 201238Arkhangelsk2664°42′43°22′Aug 22, 201339Kostamuksha264°40′30°40′August, 200740Syktyvkar1261°39′51″50°49′01″Sep 2, 202041Tomsk260°00′05″84°13′32″1560Aug 31, 202042Magadan3059°34′24″150°46′04″July 22, 202043Kemerovo2755°11′24″86°29′24″July 10, 2019Czech Republic44Krušne Hory 11250°24′56″12°40′30″925July 28, 202045Krušne Hory 21050°23′33″12°27′48″915July 28, 2020Bulgaria46Rila Mountains442°08′20″23°27′54″1950Aug 29, 2020Figure 1
Map showing the sampling localities ofCacopsylla myrtilli. The localities corresponding to these numbers are given in Table 1.


### 2.2. DNA Extraction, PCR, and Sequencing

DNeasy Blood and Tissue Kit (QIAGEN) was used to extract total genomic DNA from the material as described in our previous studies [12–14]. Complete adult individuals were homogenized in extraction buffer, and the procedure was carried out following the manufacturer’s instructions. When the yield was below 20 ng/μl, DNA was concentrated according to the standard precipitation procedure as described in Nokkala et al. [13]. A fragment of cytochrome c oxidase subunit I (COI) gene was amplified using the primer pair HybCacoCO/HybHCOMod, and PCR was carried out as earlier [12–14]. PCR products were purified with QIAquick PCR Purification Kit (QIAGEN), and the purified products were sent to Macrogen Europe (Amsterdam, the Netherlands) for Sanger sequencing. Obtained sequences were trimmed to span a 638 bp fragment of the gene and aligned using BioEdit 7.2.0 software [28]. Alignments were checked by visual inspection. New haplotype sequences found in the present study are deposited in GeneBank under the accession numbers OK174423-OK174455.


### 2.3. Data Analysis and Genealogical Relations between Haplotypes

The number of polymorphic and parsimony informative sites in the alignment of total sequences and haplotype (h) and nucleotide (π) diversities for each geographical region were calculated with DnaSP v 6.1 [29]. DnaSP was also used to investigate the demographic history of populations in each study region with Tajima’s D [30] and Fu and Li’s F [31] statistical tests for neutrality and the R2 test [32] for population size change. Significance of neutrality tests was obtained based on 1000 coalescent simulations under a model of selective neutrality. To study the genealogical relations between the haplotypes, we reconstructed a haplotype network using the median-joining algorithm [33] implemented in PopART [34] with epsilon set to 0.


## 2.1. Samples

The speciesC. myrtilli was identified using the key and description by Ossiannilsson [15]. The specimens were collected on V. myrtillus and V. uliginosum in 46 locations, 6 in Norway, 6 in Sweden, 14 in Finland, 17 in Russia, 2 in the Czech Republic, and 1 in Bulgaria during 2007-2020 (Table 1 and Figure 1) and stored immediately in 96% ethanol in the field. Partial COI gene was sequenced from 964 specimens. Because of heteroplasmy (Supplementary Figure S1), two of these were not included in the study (Table 1).Table 1
Number ofCacopsylla myrtilli individuals sequenced in samples from Norway, Sweden, Finland, Russia, Czech Republic, and Bulgaria. Altitudes below 200 m a.s.l. not shown.
Sampling locationNLatitudeLongitudeAltitude a.s.l./mDateNorway1Finnmark, Šuoššjávri3269°22′11″24°18′20″Aug 10, 20122Sjoa, Rudihøe3661°46′27″9°17′16″1000Aug 1, 20093Sjoa, Kvernbrusætrin2461°42′27″9°19′25″950Aug 1, 20094Sjoa, Stålane2461°41′15″9°14′27″1000Aug 1, 20095Sjoa, Kringlothaugen2261°43′06″9°22′40″700Aug 1, 20096Sjoa, Rindhovda 10803061°43′05″9°05′12″1080Aug 16, 2010Sweden7Abisko, Björkliden Fjället2268°24′32″18°39′55″500Aug 3, 20098Abisko, Lapporten 5702868°19′14″18°51′05″570Aug 10, 20129Soppero2168°00′39″21°39′25″450Aug 10, 201210Jokkmokk2566°35′36″19°49′20″300Aug 8, 201211Sorsele2665°29′07″17°33′36″430Aug 8, 201212Storuman2163°08′23″17°15′33″400Aug 8, 2012Finland13Utsjoki, Ailigas4569°53′51″27°03′32″320July 30, 201514Utsjoki, Hietala2469°51′06″27°00′34″Aug 11, 201215Utsjoki, Karigasniemi1769°25′14″25°58′32″450Aug 16, 201716Kilpisjärvi, Pikku-Malla1769°08′50″20°44′20″620Aug 6, 201617Sodankylä, Puisuvanto1867°46′52″26°46′09″200Aug 3, 201818Sodankylä1267°22′49″26°49′35″Aug 12, 201219Salla, Tuntsa1467°39′40″29°45′53″July 21, 201020Kuusamo1866°13′57″29°08′48″Aug 12, 201221Liminka2864°43′53″25°23′04″Aug 4, 200922Paltamo5064°33′28″27°43′41″July 16, 200923Juuka1263°15′04″29°11′12″Aug 12, 201224Suonenjoki662°40′56″27°03′48″July 17, 200925Tohmajärvi1862°23′12″30°19′12″Aug 12, 201226Forssa460°58′37″23°31′00″Aug 13, 2012Russia27Kildin Island2769°19′11″34°20′55″Aug 14, 202028Murmansk, North Nagornoje2068°55′03″33°04′01″July 16, 201229Murmansk, Kola1868°52′35″33°01′11″July 28, 200730Murmansk, Laplandskii Reserve3268°04′12″32°16′12″Aug 1, 201931Apatity2067°34′29″33°23′27″Aug 6, 201232Vorkuta3067°30′64°02′Aug 4, 201333Island Malyi Gorelyi (White Sea)1666°17′58″33°37′18″Aug 4, 201234Island Bolshoy Gorelyi (White Sea)366°17′96″33°37′30″Aug 4, 201235White Sea, mainland3366°17′40″33°39′06″Aug 18, 201636Island Sredniy (White Sea)4166°17′17″33°38′01″July 13, 201237Kartesh (White Sea)366°20′40″33°38′23″Aug 6, 201238Arkhangelsk2664°42′43°22′Aug 22, 201339Kostamuksha264°40′30°40′August, 200740Syktyvkar1261°39′51″50°49′01″Sep 2, 202041Tomsk260°00′05″84°13′32″1560Aug 31, 202042Magadan3059°34′24″150°46′04″July 22, 202043Kemerovo2755°11′24″86°29′24″July 10, 2019Czech Republic44Krušne Hory 11250°24′56″12°40′30″925July 28, 202045Krušne Hory 21050°23′33″12°27′48″915July 28, 2020Bulgaria46Rila Mountains442°08′20″23°27′54″1950Aug 29, 2020Figure 1
Map showing the sampling localities ofCacopsylla myrtilli. The localities corresponding to these numbers are given in Table 1.


## 2.2. DNA Extraction, PCR, and Sequencing

DNeasy Blood and Tissue Kit (QIAGEN) was used to extract total genomic DNA from the material as described in our previous studies [12–14]. Complete adult individuals were homogenized in extraction buffer, and the procedure was carried out following the manufacturer’s instructions. When the yield was below 20 ng/μl, DNA was concentrated according to the standard precipitation procedure as described in Nokkala et al. [13]. A fragment of cytochrome c oxidase subunit I (COI) gene was amplified using the primer pair HybCacoCO/HybHCOMod, and PCR was carried out as earlier [12–14]. PCR products were purified with QIAquick PCR Purification Kit (QIAGEN), and the purified products were sent to Macrogen Europe (Amsterdam, the Netherlands) for Sanger sequencing. Obtained sequences were trimmed to span a 638 bp fragment of the gene and aligned using BioEdit 7.2.0 software [28]. Alignments were checked by visual inspection. New haplotype sequences found in the present study are deposited in GeneBank under the accession numbers OK174423-OK174455.


## 2.3. Data Analysis and Genealogical Relations between Haplotypes

The number of polymorphic and parsimony informative sites in the alignment of total sequences and haplotype (h) and nucleotide (π) diversities for each geographical region were calculated with DnaSP v 6.1 [29]. DnaSP was also used to investigate the demographic history of populations in each study region with Tajima’s D [30] and Fu and Li’s F [31] statistical tests for neutrality and the R2 test [32] for population size change. Significance of neutrality tests was obtained based on 1000 coalescent simulations under a model of selective neutrality. To study the genealogical relations between the haplotypes, we reconstructed a haplotype network using the median-joining algorithm [33] implemented in PopART [34] with epsilon set to 0.


## 3. Results

Sequencing the 638 bp fragment of the mitochondrialCOI gene from the 962 individuals revealed 40 haplotypes (Table 2 and Figure 2). All nucleotide substitutions observed were confirmed with visual inspection of the corresponding electropherogram (Supplementary Figure S2). Of the nucleotide substitutions present in the data set, 35 were polymorphic, and 22 of these were parsimony informative. No indel mutations or stop codons were found. The transition/transversion (ti/tv) ratio was calculated from the base substitutions observed (Table 2) and was found to be 4.4. The sequence analyzed is coding for a polypeptide of 212 amino acids within the COI protein. Thirty-two haplotypes produce identical polypeptides, i.e., the underlying nucleotide substitutions were synonymous. Only eight nucleotide substitutions were nonsynonymous (haplotypes H11, H21, H22, H27, H28, H34, H36, H39) causing one amino acid change each in the polypeptide chain (Supplementary Alignment S1).Table 2
Haplotypes H1-H40, the distribution of derived haplotypes, haplotype-specific substitution(s), and accession numbers of each haplotype. Each derived haplotype is named according to one of the main haplotypes, Rindhovda (H1), Rudihoe (H2), or Utsjoki (H3) from which it is derived. Nor: Norway; Swe: Sweden; Fin: Finland; Rus: Russia; Cze: Czech Republic.
HaplotypeDistributionSubstitution(s)Accession numberReferenceH1KF494327.1Nokkala et al. 2015 [12]H2KF494326.1Nokkala et al. 2015 [12]H3KF494329.1Nokkala et al. 2019 [14]H4NorH1_92A > GOK174423This studyH5Nor, SweH2_191T > GOK174424This studyH6NorH2_191T > G +599C > TOK174425This studyH7NorH1_320G > AOK174426This studyH8NorH2_347C > TOK174427This studyH9Nor, RusH1_458C > TOK174428This studyH10NorH2_566A > TOK174429This studyH11NorH2_566A > T +319 T > COK174430This studyH12NorH2_581A > GOK174431This studyH13NorH1_605A > GKF494330.1Nokkala et al. 2019 [14]H14NorH1_605A > G +218C > TOK174432This studyH15NorH1_608G > AOK174433This studyH16SweH2_11G > AOK174434This studyH17FinH1_98C > TOK174435This studyH18FinH3_146C > TKF494328.1Nokkala et al. 2019 [14]H19FinH3_146C > T +11G > AOK174436This studyH20FinH1_182A > GOK174437This studyH21FinH1_183G > AOK174438This studyH22FinH1_501T > AOK174439This studyH23FinH1_563A > GOK174440This studyH24FinH2_569T > GOK174441This studyH25Fin, RusH2_636C > TOK174442This studyH26RusH2_17A > GOK174443This studyH27RusH2_17A > G +478 T > COK174444This studyH28RusH2_63C > GOK174445This studyH29RusH2_83G > AOK174446This studyH30RusH2_185C > TKF494331.1Nokkala et al. 2019 [14]H31RusH2_185C > T +311A > GKF494332.1Nokkala et al. 2019 [14]H32RusH3_224G > A +254C > TOK174447This studyH33RusH2_299A > GOK174448This studyH34RusH1_319T > COK174449This studyH35RusH1_323C > TOK174450This studyH36RusH3_372G > AOK174451This studyH37RusH1_461A > GOK174452This studyH38RusH2_461A > GOK174453This studyH39RusH1_466C > AOK174454This studyH40CzeH2_461A > COK174455This studyFigure 2
(a) Median-joining haplotype network based on theCOI haplotypes of C. myrtilli. Mutational steps are shown by the number of hatch marks. Circle sizes are proportional to haplotype frequency. (b) Distribution of the COI haplotypes among the populations of C. myrtilli studied.
(a)(b)Derived haplotypes are formed when a mutation occurs in one of the main haplotypes. If there is a mutation in H1 haplotype at, e.g., site 92 changing A to G, the derived haplotype is labeled as H1_92A > G (Table2). Most of the derived haplotypes were caused by a single nucleotide substitution in one of the main haplotypes. Seven haplotypes were found, where a second substitution occurred along with a known single substitution. In one case with two substitutions (H32, derived from H3), no intermediate haplotype (with a single substitution) was found. A maximum difference of five nucleotide substitutions was found between the most divergent haplotypes. The genetic diversity observed in haplotype networks (Figure 2(a)) within each country is summarized in Figure 2(b) and Table 3. As a rule, each derived haplotype represented a private haplotype found in only one country, indicating separate genealogical histories in different geographical regions. The only exceptions are H5 near Norwegian-Swedish and H25 near Finnish-Russian borders in these countries, representing geographically neighboring localities. A quite different situation is H9, frequent in Norway but also found in single specimens in Laplandskii Reserve and Srednyi Island populations in northwest Russia. Two regions stand out from the data, showing relatively high haplotype diversity, the Sjoa region in Norway with 14 and the White Sea mainland, in Russia with 8 haplotypes. In most populations, only 1-2 derived haplotypes, if any, were present, in addition to the main haplotype(s). Detailed data on haplotype diversity in the populations studied is given in Supplementary material (Table S1). All 130 cases of derived haplotypes included in the data set were confirmed by checking the corresponding electropherograms. In two cases (0.21%), it was evident that the change in nucleotide sequence was not genuine but caused by heteroplasmy. These cases were excluded from further analysis.Table 3
Summary of haplotype diversity for the sampling regions. The number of individuals sequenced (N), the number of segregating sites (S), the number of different haplotypes (H), haplotype (h), and nucleotide (π) diversities with standard deviation (SD) are given. The number of individuals with main haplotypes is shown in separate columns, and unique derived haplotypes are listed for each locality (number of individuals, if >1). Tajima’s D, Fu’s F, and R2 of Ramos-Onsins and Rozas values are shown.
RegionNSHMain haplotypesDerived haplotypesh (SD)π (SD)Tajima’sDFu’s FsR2H1H2H3Norway16814144966H4, H5 (18), H6, H7, H8 (8), H9 (11), H10 (3), H11, H13 (4), H15 (2), H12 (2), H140.746 (0.022)0.00251 (0.00012)-0.89178-4.0480.0560Sweden1434431109H5 (2), H160.382 (0.040)0.00116 (0.00012)0.049140.834Finland28311121071261H17 (19), H18, H19, H20, H21 (19), H22, H23 (2), H24, H25 (4)0.652 (0.017)0.00209 (0.00006)-0.56195-2.6720.0604Russia3421718144174H9 (2), H25 (2), H26, H27, H28 (2), H29, H30, H31, H32, H33, H34 H35 (4), H36 (2), H37 (2), H38, H390.565 (0.014)0.0018 (0.00005)-1.40125-9.6010.0339Czech Republic221221H400.091 (0.081)0.00014 (0.00013)-1.1624-0.957Bulgaria41400N.A.N.A.Since the sample locations covered an extensive region from Norway in the west to Magadan in the Russian Far East, it was considered important to study if there are differences in haplotype distribution in various geographical regions, and a haplotype network was constructed (Figure2(a)). The network revealed three main haplotypes in C. myrtilli, the Rindhovda (H1), and Rudihoe (H2) haplotypes already reported previously [12]and a third main haplotype in Utsjoki (H3). These main haplotypes are characterized by specific nucleotides at positions 170 and 311 in the fragment sequenced. In H1, these are C and G, respectively, in H2, T and A, and in H3, C and A (Table 4). The network shows star-like topologies with main haplotypes forming the center, which indicates recent population expansion. The main haplotypes H1 and H2 are extremely common and geographically widespread. Either one or both occurred in high frequency in all sampling locations studied. The H2 haplotype was more common than the H1 haplotype (54.0% and 34.4%, respectively, of all 962 specimens studied).Table 4
Nucleotides at positions 170 and 311 characterizing three main haplotypes H1, H2, and H3.
Nucleotide position170311H1 (Rindhovda)CGH2 (Rudihoe)TAH3 (Utsjoki)CAIn populations from Northern Europe (Fennoscandia and Russia), haplotype diversity was relatively high, varying from 0.382 to 0.746, and nucleotide diversity was low, ranging from 0.0116 to 0.00251 (Table3). In these populations, haplotype diversity was the highest in Norway and lowest in Sweden. The lowest haplotype diversity was observed in populations from Sorsele in Sweden, Utsjoki-Hietala, Kuusamo, Suonenjoki, Tohmajärvi, and Forssa in Finland, Vorkuta, Bolshoy Gorely Island, and Tomsk in Russia, where only one haplotype was found. In Central Europe, in a population from Krušne Hory mountains in Czech Republic, one private derived haplotype was found in addition to H2 haplotype, whereas in Southern Europe in samples from Rila Mountains in Bulgaria, only the main haplotype H2 was found. Although the star-like topology of the haplotype networks and negative values of Tajima’s D and Fu’s S neutrality tests suggested recent population expansion and an excess of singleton haplotypes in Norway Finland, and Russia, the values were significant only for the Russian populations. The R2 values indicative of population size change were not significant but close to it.The data set presented in Table2 allowed us to estimate the nucleotide substitution rate in parthenogenetic C. myrtilli. If we leave out the main haplotypes and haplotypes derived from the H3 haplotype, 33 derived substitutions have appeared in the 10 000 years passed since the ending of the last glacial period in the 638 bp fragment of the COI gene sequenced in 956 individuals. As C. myrtilli has only one generation per year, we can estimate a rate of 5.4×10−9 substitutions per site per generation.


## 4. Discussion

### 4.1. Geographical Variation in Haplotype Diversity

Of the three main haplotypes found inC. myrtilli, H1 (Rindhovda) and H2 (Rudihøe) are extremely common throughout the Palearctic localities studied. Either H1 or H2 is the most common haplotype in all geographical regions studied. There is, however, a significant difference in the distributions of H1 and H2. The H1 haplotype is the prevailing one in high altitude and high latitude populations, as well as in the harshest environments, e.g., Utsjoki populations in Finland, Abisko in Sweden, Finnmark, and sampling sites in the Sjoa region in Norway, Syktyvkar, and Magadan in Russia. The H1 haplotype is clearly associated with adaptation to the most extreme environments. The H2 haplotype, in turn, is dominant in lower altitudes and latitudes and is, in fact, the most common haplotype overall. The proportion of H2 haplotype increases towards the south and (with its derivatives) is the only haplotype found in several populations (Figure 2(b), Supplementary Table S1).The third main haplotype H3 (Utsjoki) and its derivatives showed a restricted distribution in northern Finland and northwest Russia. A single male in Utsjoki displayed the main haplotype H3. The derivative haplotypes were carried by two males from Paltamo, Finland (H18 or H19), one male from White Sea Mainland, Russia (H32), and two females from Arkhangelsk, Russia (H36). The H3 haplotype is highly interesting as a transition type of nucleotide substitution in position 311 would yield H1 and a similar type of substitution in position 170 would yield H2. The H3 haplotype is most probably the ancestral one inC. myrtilli as the very same nucleotides appear in corresponding positions in the closely allied parthenogenetic species C. borealis and C. ledi, as well as in the bisexual species C. lapponica [14]. Nokkala et al. [14] also showed the basal position of the H3 haplotype within C. myrtilli phylogeny. In addition, males carrying the H3 haplotype or its derivatives display normal spermatogenesis, i.e., they are functional males [11]. All these attributes allow us to suggest that these psyllids form a distinct lineage well separated from the triploid parthenogenetic linage and actually represent bisexually reproducing C myrtilli that participated in generating the triploid parthenogenetic C. myrtilli. However, it remains unclear how infrequent diploids can persist within parthenogenetic populations. It is possible that the existence of mixed parthenogenetic and bisexual populations is restricted to northern Europe and exclusively bisexual populations would presumably be found in southern refugial regions where they could have survived during the Pleistocene glacial period.Except for the haplotype H3 and its derivatives, all remaining new haplotypes are derived from either H1 or H2 main haplotypes, by one or two nucleotide substitution(s). These new haplotypes, private or shared with close neighboring populations, appear to be quite specific to each country. Different distribution patterns of the main and derived haplotypes strongly suggest that the main haplotypes are of refugial origin, whereas the derived haplotypes have emerged after the recolonization process when the lineages have already reached their present-day distributions in Holocene. This is also supported by the star-like haplotype networks and statistical analyses. The only exception is provided by H9 that was found both in Norway and northwest Russia (see below). The discovery of heteroplasmic specimens is a sign of ongoing evolution. Originally, each nucleotide substitution has appeared as a heteroplasmic mutation and has become later fixed producing a new nucleotide substitution and a new haplotype. Xu et al. [35] have studied the evolution of the mitochondrial genome and estimated that the fixation of a new mutation takes several generations. During this time, the two types of mitochondria coexist in cells and cause heteroplasmy.It is evident that genetic diversity is the highest in the Sjoa district, Norway, where samples were collected in five high altitude sampling locations situated in a quite small area of9×16km. In addition to the two main haplotypes, H1 and H2, 12 derived haplotypes of which three included two nucleotide substitutions were encountered in Sjoa. Also, in the Russian White Sea mainland population, diversity was moderately high with six derived haplotypes in addition to the main haplotypes H1 and H2. On the other hand, the lowest genetic diversity in northern Europe was found in Sweden where only one derived haplotype not found in any other location was encountered (Figure 2(b) and Table 3), suggesting that the populations in the Sjoa and White Sea mainland are of considerably more ancient origin than populations in Sweden.


### 4.2. Recolonization Routes ofC. myrtilli

The explanation for the existence of the highest genetic diversity in the Sjoa region is probably tightly coupled with the recolonization history of the species and influenced by three cooperating factors. Firstly, it seems apparent thatC. myrtilli has recolonized southern Norway much earlier than northwest Russia, including the Kola Peninsula and Finland. The last region, which C. myrtilli has recolonized, is Sweden, where genetic diversity is extremely low. It is clear that these lines of evidence reflect the phases of retreating ice cover and evolutionary phases of the Baltic Sea. Most likely, colonization from a western refugium through Denmark along the land connection and through an ice-free corridor to southern Norway could occur as early as during the Baltic Ice Lake phase (ca. 13000 years BP), when central and northern Sweden were still covered with thick ice cover. The suggestion parallels the colonizing history of the freshwater amphipod Gammarus lacustris in northern Europe [36]. This colonization route is the oldest one and is reflected in the high genetic diversity in the Sjoa region.Secondly, the highest genetic diversity in Norway could imply that the derived haplotypes present in Sjoa emerged already in a refugium with the main haplotypes H1 and H2. To some extent, this seems to be true. The derived haplotype H9 from Sjoa was also found in populations from northwest Russia, the White Sea mainland, and Laplandskii Reserve near Murmansk. It seems plausible that the derived haplotype H9 in both Sjoa and the White Sea region in Russia originated from the same western refugium. At present, we have no means to determine if any other haplotypes present in Sjoa have their origin in the same refugium.Thirdly, we recently surveyed the presence of theWolbachia bacterial endosymbiont in five Cacopsylla species [37]. The study area covered the same regions as in the present work, and several C. myrtilli populations surveyed were the very same as in the present study. Our analysis showed that all sampling sites in Sjoa were Wolbachia-free whereas samples from Sweden, Finland, and Russia were moderately infected [37]. This suggests that C. myrtilli in western refugia were Wolbachia-free. This refugium could be located somewhere in Western Europe, e.g., in the Belgian Ardennes, a region demonstrated to be one of the cryptic northern refugia (reviewed by [23]) or more south, e.g., in the Iberian Peninsula.It is evident that there must also be separate eastern refugial regions from whichC. myrtilli clones carrying Wolbachia infections presently distributed in Russia, Finland, and Sweden have originated. Shapoval et al. [37] found that C. fraudatrix, a species closely related to C. myrtilli [38], carried similar Wolbachia strains as C. myrtilli. Samples of this species were collected from Biezczady Mountains in southeast Poland in the western part of the large Carpathian Mountain area. Therefore, we consider it plausible to suggest that the Carpathian Mountains represent one eastern refugium from which Wolbachia positive C. myrtilli originate. Provan and Bennett [24] have reviewed evidence for the existence of a cryptic northern refugium in this area. However, there must also exist an additional eastern Wolbachia-free refugium or refugia, since Syktyvkar, Czech Republic (Krušne Hory), Bulgarian (Rila Mountains), and Kemerovo (western Siberia) populations were not infected by Wolbachia. The Balkans could have served as a refugium, but only very few C. myrtilli females from Rila Mountains, Bulgaria, were included in the present study, thus not allowing firm conclusions. In addition to the Balkans, also Caucasus, shown to have been a refugium during the glacial period (see [39]), could have served as an eastern Wolbachia-free refugium for C. myrtilli.It is important to note that in populations infected byWolbachia, the bacterial infection can modify genetic diversity by causing selective sweep [40–44]. That is why our estimation for the nucleotide substitution rate in the parthenogenetic form of C. myrtilli could be somewhat biased, as in populations in Sweden, Finland, and Russia the prevalence of Wolbachia infections is approximately 23% [37]. More accurate rate can be determined only after the distribution of haplotypes within refugial regions have been analyzed.The eastern recolonization route first reached the White Sea region from where psyllids were distributed further north to the Kola Peninsula and northern Finland, towards the west to Finland and later towards southern Finland as the ice retreated. Northern Norway, Finnmark, was recolonized via northern Finland and northern and Central Sweden from the north. This colonization route is strongly supported by the similarity ofWolbachia infections in northern Finland, Finnmark, and Sweden, as well as the extremely shallow genetic diversity in these populations. The western and eastern recolonization lineages came into secondary contact in Central Sweden forming a hybrid zone. The region from central to northern Sweden is known to harbor several hybrid zones [19, 45–48], as it was the last region in Europe to dispose of ice cover ca. 9000 BP. However, the presence of a shared COI haplotype H5 in Sjoa and Storuman populations in Central Sweden is evidence that the western lineage originated from Norway. A visual summary of the suggested recolonization routes of C. myrtilli is presented in Figure 3. Hewitt [19, 20] introduced three broad patterns of recolonization described as paradigms, named as “grasshopper,” “hedgehog,” and “bear,” according to model case studies, and gave several examples of species’ recolonization histories falling in one of those paradigms. Our scenario for C. myrtilli parallels quite well with the “bear” paradigm; although, we suggest the western recolonization route of C. myrtilli to have taken an even more western path than outlined in Hewitt’s [19, 20] “bear” paradigm.Figure 3
Hypothetical recolonization routes ofC. myrtilli in Europe after the last glacial maximum. Western routes are shown in red and eastern routes in blue.While considering the phylogeographic history of the eastern Palearctic populations of Kemerovo, Tomsk, and Magadan, we must consider alternative refugia. It is evident that Magadan and most probably Kemerovo and Tomsk have been recolonized from the Siberian and/or Manchurian refugia. However, we cannot totally rule out the possibility that recolonization from a more western refugium or refugia, e.g., Caucasus, would have reached Kemerovo and Tomsk since there are examples such as the case of the adderVipera berus, where one of the postglacial recolonization routes is assumed to have originated from an extra-Mediterranean refugium and led to Siberia [49].


## 4.1. Geographical Variation in Haplotype Diversity

Of the three main haplotypes found inC. myrtilli, H1 (Rindhovda) and H2 (Rudihøe) are extremely common throughout the Palearctic localities studied. Either H1 or H2 is the most common haplotype in all geographical regions studied. There is, however, a significant difference in the distributions of H1 and H2. The H1 haplotype is the prevailing one in high altitude and high latitude populations, as well as in the harshest environments, e.g., Utsjoki populations in Finland, Abisko in Sweden, Finnmark, and sampling sites in the Sjoa region in Norway, Syktyvkar, and Magadan in Russia. The H1 haplotype is clearly associated with adaptation to the most extreme environments. The H2 haplotype, in turn, is dominant in lower altitudes and latitudes and is, in fact, the most common haplotype overall. The proportion of H2 haplotype increases towards the south and (with its derivatives) is the only haplotype found in several populations (Figure 2(b), Supplementary Table S1).The third main haplotype H3 (Utsjoki) and its derivatives showed a restricted distribution in northern Finland and northwest Russia. A single male in Utsjoki displayed the main haplotype H3. The derivative haplotypes were carried by two males from Paltamo, Finland (H18 or H19), one male from White Sea Mainland, Russia (H32), and two females from Arkhangelsk, Russia (H36). The H3 haplotype is highly interesting as a transition type of nucleotide substitution in position 311 would yield H1 and a similar type of substitution in position 170 would yield H2. The H3 haplotype is most probably the ancestral one inC. myrtilli as the very same nucleotides appear in corresponding positions in the closely allied parthenogenetic species C. borealis and C. ledi, as well as in the bisexual species C. lapponica [14]. Nokkala et al. [14] also showed the basal position of the H3 haplotype within C. myrtilli phylogeny. In addition, males carrying the H3 haplotype or its derivatives display normal spermatogenesis, i.e., they are functional males [11]. All these attributes allow us to suggest that these psyllids form a distinct lineage well separated from the triploid parthenogenetic linage and actually represent bisexually reproducing C myrtilli that participated in generating the triploid parthenogenetic C. myrtilli. However, it remains unclear how infrequent diploids can persist within parthenogenetic populations. It is possible that the existence of mixed parthenogenetic and bisexual populations is restricted to northern Europe and exclusively bisexual populations would presumably be found in southern refugial regions where they could have survived during the Pleistocene glacial period.Except for the haplotype H3 and its derivatives, all remaining new haplotypes are derived from either H1 or H2 main haplotypes, by one or two nucleotide substitution(s). These new haplotypes, private or shared with close neighboring populations, appear to be quite specific to each country. Different distribution patterns of the main and derived haplotypes strongly suggest that the main haplotypes are of refugial origin, whereas the derived haplotypes have emerged after the recolonization process when the lineages have already reached their present-day distributions in Holocene. This is also supported by the star-like haplotype networks and statistical analyses. The only exception is provided by H9 that was found both in Norway and northwest Russia (see below). The discovery of heteroplasmic specimens is a sign of ongoing evolution. Originally, each nucleotide substitution has appeared as a heteroplasmic mutation and has become later fixed producing a new nucleotide substitution and a new haplotype. Xu et al. [35] have studied the evolution of the mitochondrial genome and estimated that the fixation of a new mutation takes several generations. During this time, the two types of mitochondria coexist in cells and cause heteroplasmy.It is evident that genetic diversity is the highest in the Sjoa district, Norway, where samples were collected in five high altitude sampling locations situated in a quite small area of9×16km. In addition to the two main haplotypes, H1 and H2, 12 derived haplotypes of which three included two nucleotide substitutions were encountered in Sjoa. Also, in the Russian White Sea mainland population, diversity was moderately high with six derived haplotypes in addition to the main haplotypes H1 and H2. On the other hand, the lowest genetic diversity in northern Europe was found in Sweden where only one derived haplotype not found in any other location was encountered (Figure 2(b) and Table 3), suggesting that the populations in the Sjoa and White Sea mainland are of considerably more ancient origin than populations in Sweden.


## 4.2. Recolonization Routes ofC. myrtilli

The explanation for the existence of the highest genetic diversity in the Sjoa region is probably tightly coupled with the recolonization history of the species and influenced by three cooperating factors. Firstly, it seems apparent thatC. myrtilli has recolonized southern Norway much earlier than northwest Russia, including the Kola Peninsula and Finland. The last region, which C. myrtilli has recolonized, is Sweden, where genetic diversity is extremely low. It is clear that these lines of evidence reflect the phases of retreating ice cover and evolutionary phases of the Baltic Sea. Most likely, colonization from a western refugium through Denmark along the land connection and through an ice-free corridor to southern Norway could occur as early as during the Baltic Ice Lake phase (ca. 13000 years BP), when central and northern Sweden were still covered with thick ice cover. The suggestion parallels the colonizing history of the freshwater amphipod Gammarus lacustris in northern Europe [36]. This colonization route is the oldest one and is reflected in the high genetic diversity in the Sjoa region.Secondly, the highest genetic diversity in Norway could imply that the derived haplotypes present in Sjoa emerged already in a refugium with the main haplotypes H1 and H2. To some extent, this seems to be true. The derived haplotype H9 from Sjoa was also found in populations from northwest Russia, the White Sea mainland, and Laplandskii Reserve near Murmansk. It seems plausible that the derived haplotype H9 in both Sjoa and the White Sea region in Russia originated from the same western refugium. At present, we have no means to determine if any other haplotypes present in Sjoa have their origin in the same refugium.Thirdly, we recently surveyed the presence of theWolbachia bacterial endosymbiont in five Cacopsylla species [37]. The study area covered the same regions as in the present work, and several C. myrtilli populations surveyed were the very same as in the present study. Our analysis showed that all sampling sites in Sjoa were Wolbachia-free whereas samples from Sweden, Finland, and Russia were moderately infected [37]. This suggests that C. myrtilli in western refugia were Wolbachia-free. This refugium could be located somewhere in Western Europe, e.g., in the Belgian Ardennes, a region demonstrated to be one of the cryptic northern refugia (reviewed by [23]) or more south, e.g., in the Iberian Peninsula.It is evident that there must also be separate eastern refugial regions from whichC. myrtilli clones carrying Wolbachia infections presently distributed in Russia, Finland, and Sweden have originated. Shapoval et al. [37] found that C. fraudatrix, a species closely related to C. myrtilli [38], carried similar Wolbachia strains as C. myrtilli. Samples of this species were collected from Biezczady Mountains in southeast Poland in the western part of the large Carpathian Mountain area. Therefore, we consider it plausible to suggest that the Carpathian Mountains represent one eastern refugium from which Wolbachia positive C. myrtilli originate. Provan and Bennett [24] have reviewed evidence for the existence of a cryptic northern refugium in this area. However, there must also exist an additional eastern Wolbachia-free refugium or refugia, since Syktyvkar, Czech Republic (Krušne Hory), Bulgarian (Rila Mountains), and Kemerovo (western Siberia) populations were not infected by Wolbachia. The Balkans could have served as a refugium, but only very few C. myrtilli females from Rila Mountains, Bulgaria, were included in the present study, thus not allowing firm conclusions. In addition to the Balkans, also Caucasus, shown to have been a refugium during the glacial period (see [39]), could have served as an eastern Wolbachia-free refugium for C. myrtilli.It is important to note that in populations infected byWolbachia, the bacterial infection can modify genetic diversity by causing selective sweep [40–44]. That is why our estimation for the nucleotide substitution rate in the parthenogenetic form of C. myrtilli could be somewhat biased, as in populations in Sweden, Finland, and Russia the prevalence of Wolbachia infections is approximately 23% [37]. More accurate rate can be determined only after the distribution of haplotypes within refugial regions have been analyzed.The eastern recolonization route first reached the White Sea region from where psyllids were distributed further north to the Kola Peninsula and northern Finland, towards the west to Finland and later towards southern Finland as the ice retreated. Northern Norway, Finnmark, was recolonized via northern Finland and northern and Central Sweden from the north. This colonization route is strongly supported by the similarity ofWolbachia infections in northern Finland, Finnmark, and Sweden, as well as the extremely shallow genetic diversity in these populations. The western and eastern recolonization lineages came into secondary contact in Central Sweden forming a hybrid zone. The region from central to northern Sweden is known to harbor several hybrid zones [19, 45–48], as it was the last region in Europe to dispose of ice cover ca. 9000 BP. However, the presence of a shared COI haplotype H5 in Sjoa and Storuman populations in Central Sweden is evidence that the western lineage originated from Norway. A visual summary of the suggested recolonization routes of C. myrtilli is presented in Figure 3. Hewitt [19, 20] introduced three broad patterns of recolonization described as paradigms, named as “grasshopper,” “hedgehog,” and “bear,” according to model case studies, and gave several examples of species’ recolonization histories falling in one of those paradigms. Our scenario for C. myrtilli parallels quite well with the “bear” paradigm; although, we suggest the western recolonization route of C. myrtilli to have taken an even more western path than outlined in Hewitt’s [19, 20] “bear” paradigm.Figure 3
Hypothetical recolonization routes ofC. myrtilli in Europe after the last glacial maximum. Western routes are shown in red and eastern routes in blue.While considering the phylogeographic history of the eastern Palearctic populations of Kemerovo, Tomsk, and Magadan, we must consider alternative refugia. It is evident that Magadan and most probably Kemerovo and Tomsk have been recolonized from the Siberian and/or Manchurian refugia. However, we cannot totally rule out the possibility that recolonization from a more western refugium or refugia, e.g., Caucasus, would have reached Kemerovo and Tomsk since there are examples such as the case of the adderVipera berus, where one of the postglacial recolonization routes is assumed to have originated from an extra-Mediterranean refugium and led to Siberia [49].


## 5. Conclusions

The present-day haplotype diversity alone does not allow concluding whether the western and eastern recolonizers originally expanded from a single refugium each or if two or more refugia were involved. However, the absence ofWolbachia infection in sampling locations in the Sjoa region in Norway, together with moderate levels of infection in other geographical regions, Sweden, Finland, and Russia, suggests that the C. myrtilli in Sjoa has originated from a separate, presumably Wolbachia-free western refugium or refugia. However, some recolonization from the western refugia has also reached the White Sea region in Russia as evidenced by the presence of the H9 derived haplotype. It is also evident that the Sjoa district has been the first to be recolonized, as coastal Norway was deglaciated already ca. 13 000 BP, creating an ice-free corridor for early recolonization. At the same time, western Russia south of the White Sea was also deglaciated and applicable for recolonization. The White Sea region was probably also recolonized early, possibly in multiple recolonizing events and certainly from multiple refugia. In addition, we cannot rule out that C. myrtilli might have inhabited a cryptic refugium on the Norwegian coast, as evidence of its existence are presented [23, 50, 51], and hence also Norway would have been recolonized from multiple sources. To achieve a more comprehensive understanding of the postglacial recolonization process of parthenogenetic C. myrtilli, a systematic study of haplotype structures in putative refugial regions is needed.



---

*Source: 5458633-2022-07-08.xml*