Microsatellite variation and population genetic structure of Anatolian mountain frogs

,

Genetic diversity is one of the most notable types of biodiversity, because it demonstrates the population variation response to environmental changes (Frankham 1995). This type of biodiversity is indispensable for assessing population trends and for species conservation. Resource-rich environments generally incorporate larger populations with high levels of genetic diversity and species richness. Contrarily, habitat loss and isolation are essential features in determining the sizes of smaller populations, which leads to a genetic decrease in these populations (Dolgener et al. 2012). The past evolutionary process also has an important impact on the population genetic structure of organisms. Primarily, the glacial and inter-glacial ages played an important role in biodiversity at mid and high latitudes and altitudes (Hewitt 2000;Taberlet et al. 1998). Due to their relatively low dispersal capacity, amphibians are good representative species for studying the genetic differentiation of isolated geographical populations (Stuart 2006). Anatolia has a unique biodiversity level in terms of amphibians. Anatolian mountain frogs, which are members of the amphibian fauna of Anatolia, have distinctive features such as an adaption to cold and high mountains. Although researchers have formerly described four amphibian species of Anatolian mountain frogs (Rana macrocnemis, Rana camerani, Rana holtzi and Rana tavasensis), recent studies have shown that this group includes only two species: Rana macrocnemis and Rana tavasensis (Veith et al. 2003;Ergül Kalayci et al. 2017).
Rana macrocnemis Boulenger, 1885 was first found in Bursa (Türkiye) and is widely distributed population genetic structure of the Anatolian mountain frogs using six different microsatellite markers from 31 different locations in Anatolia and demonstrated the degree of sub-population division.

Sampling
We collected 138 Anatolian mountain frogs (129 for R. macrocnemis and 9 for R. tavasensis) from 31 different localities in Türkiye (Fig.1, Table 1). The sampling was carried out with the permission of the Local Ethics Committee for Animal Experiments (Approval Reference Number: 2014/36). Muscle tissue samples were obtained by toe-clipping and were stored in 95% alcohol for further analyses. We pooled our populations for the population genetic analysis based on a Structure analysis, which constructs clusters according to the genetic relatedness among those groups (Fig. 1). Based on the Structure analysis, individuals from the Black Sea, Central Anatolia and Eastern regions of Türkiye were categorised as the East Anatolia cluster; while individuals from the Aegean region, Marmara region and some localities (Bolu and Kastamonu) from the west side of the Black Sea regions were categorised as the West Anatolia cluster. The Central Taurus Mountain, including its localities, was categorised as Central Taurus, which is also known as terra typica of R. holtzi (Table 1).

DNA Isolation and PCR
We extracted genomic DNA using the DNA Purification Kit (Promega, Madison, USA) according to the manufacturer's instructions. The following in the forest and subalpine belts of the Caucasus and Asia Minor, to Iran and Southwestern Turkmenia. Meanwhile, Rana tavasensis Baran and Atatür, 1986 is endemic to West Anatolia, and its distribution is restricted to the area around Tavas (Denizli), Girdev Lake (Muğla) and the adjacent regions. R. tavasensis is categorised as endangered in the International Union for Conservation of Nature (IUCN) Red List database.
Despite some molecular methodologies being used for the purpose of highlighting the genetic relationships among Anatolian mountain frog populations (Veith et al. 2003;Ergül Kalayci et al. 2017;Picariello et al. 1999Picariello et al. , 2016Picariello et al. , 2018Najibzadeh et al. 2017aNajibzadeh et al. , 2021, no study has yet been conducted on the genetic relationships among the members of this group using high polymorphic markers such as microsatellites. Microsatellite markers are codominant nuclear DNA markers that are useful for inferring demographic processes, which consist of repeating mono-, di-, tri-and tetranucleotide units distributed throughout the genomes of most eukaryotic species (Ellegren 2004). We need highly polymorphic markers to clarify the genetic structure among Anatolian mountain frog populations, for determining the biodiversity and conservation status. Moreover, Najibzadeh et al. (2021) highlighted the lack of population genetic studies about Anatolian mountain frogs.
The aims of the present study are to: (1) assess the genetic diversity and population structures of Anatolian mountain frogs; (2) reveal the inbreeding level of each cluster and the degree of genetic differentiation within and among its populations; and (3) detect whether a possible bottleneck event has occurred in each cluster. Following our aims, we compared the Table 1 Sampling information and constructed groups. The number of specimens collected from the exact locality is given in parentheses following the geographical name; numbers in square brackets indicate the locality number in Fig. 1; N   C -Estimation of the population using lnP(D)-derived delta K with the cluster number (K) ranging from one to ten.

Population Genetic Analyses
One hundred thirty-eight individuals from 31 locations were genotyped for six microsatellite loci. The population genetic structure was evaluated by an analysis of molecular variance (AMOVA) using ARLEQUIN v. 3.5.2.2 (Excoffier & Lischer 2010). We estimated the polymorphism information content (PIC value) using the CERVUS 3.0. program (Kalinowski et al. 2007). Estimates of the diversity, number, frequency of alleles and the number of private alleles in the populations were determined using GENALEX 6.5 (Peakall & Smouse 2012). We also used ARLEQUIN v. 3.5.2.2 to test the deviation from the Hardy-Weinberg equilibrium (HWE). The expected and observed heterozygosity and allelic richness were calculated for each cluster using FSTAT 2.9.4 (Goudet 2003).
The genetic differentiation (F ST ) and linkage disequilibrium between pairs of loci was calculated using GENEPOP (Rousset 2008). The inbreeding coefficient (F IS ) varied between −1 and 1. F IS = 0 indicated no inbreeding, while values greater than '0' indicated non-random mating because of inbreeding. We calculated the F IS per loci using FSTAT 2.9.4.
We performed a Principal Coordinate Analysis (PCoA) to determine the population structure based on the F ST genetic distance for individuals using GE-NALEX (Peakall & Smouse 2012).
Every forward primer of each locus was fluorescently labelled with FAM, NED, VIC and PET dyes ( Table 2). The multiplex polymerase chain reaction (PCR) was run using the Type-it multiplex PCR Kit (Qiagen, Hilden, Germany) in 10 µl of the reaction mixture. The PCR reaction mixtures contained 3.0 µl of Master Mix, 0.4 µl of each fluorescently labelled primer (10 µm), and 1 µl of DNA (10-20 ng/µl), while the rest of the mixture consisted of RNase-free water. It was run with the following PCR conditions: initial denaturation at 95°C for 5 min, 32 cycles at 95°C for 30 s, 48-60°C (varied according to the primer annealing temperature) for 90 s, 72°C for 30 s, and a final holding at 60°C for 30 min. The microsatellite analysis was performed using the ABI3730XL capillary analyser. The alleles were scored using the GENEMARKER program (Soft Genetics LLC). were observed in the Western (2) and Eastern Anatolia (7) clusters. We found that the allelic richness of each locus ranged from 2.81 to 11.94, the mean observed heterozygosity was 0.55±0.12, 0.53±0.09 and 0.58±0.06, and the expected heterozygosity was 0.63±0.13, 0.62±0.07 and 0.72±0.08 in the Western, Central and Eastern Anatolia clusters, respectively. The mean number of alleles ranged from 7.00±1.67, 5.17±1.30 and 9.33±2.12 in all the clusters. The observed values of heterozygotes in most of the populations were lower than those expected, according to HWE. The BFG095, BFG134 and RadalE8 locus showed a significant deviation from the HWE in Western and Central Anatolia, and all the loci significantly deviated from the HWE in Eastern Anatolia (Table 3). We found the mean F IS values among all the clusters to be 0.15, 0.18 and 0.20 (Table 3). The low F ST values among all the clusters ranged from 0.078 to 0.105.
Population genetic variations were determined hierarchically at four levels: within individuals, among populations within groups (14.77%), among populations within individuals (5.16%) and among groups (8.56%). The test showed that the within individuals variation accounted for 71.51% of the total variations in the dataset. No significant structuring was found among the groups by the AMOVA test.
We detected significant heterozygotes excess for all clusters under the IAM model. However, no significant heterozygotes excess was found for all clusters under the SMM model. Members of the East Anatolia cluster were the only ones that showed significant heterozygotes excess under the TPM model. These results were inconsistent with the normal L-shaped distribution of allele frequency, indicating no genetic bottleneck in any of the three clusters in the recent past (except the Central Taurus group) (Table 4).

Discussion
Anatolian mountain frogs are categorised into cold-adapted and recently-diverged groups. Notwithstanding the different researchers (Veith et al. 2003;Ergül Kalayci et al. 2017;Najibzadeh et al. 2017aNajibzadeh et al. , 2021 who have explicated the phylogenetic relationship of the R. macrocnemis group, the elaborate population genetic structure of this group still remains a mystery. the populations were used, and three independent runs for K between 1 and 10 with 15,000,000 iterations (with 500,000 removed as a burn-in) were performed. We detected genetically distinct clusters (K) among the groups based on the similarities between the genotypes of the animals. The best-fit K value was determined using ad hoc statistic-based approach implemented in the STRUCTURE HAR-VESTER v0.6 (Earl & VonHoldt 2012;Evanno et al. 2005) software program. We used CLUMPAK (Kopelman et al. 2015) to generate aligned individual assignment bar graphs.
We determined the recent demographic bottleneck effect on the populations by examining the significant heterozygosity excess or deficit using the BOT-TLENECK program (Piry et al. 1999). We used Wilcoxon methods to test mutations through the three mutation models (infinite allele model -IAM, stepwise mutation model -SMM and two-phase model -TPM). Furthermore, a qualitative test of the modal shift was performed to evaluate the frequency distribution of alleles at different microsatellite loci.
No significant presence of null alleles, scoring errors and allelic dropout was observed. However, we found linkage disequilibria (p<0.05) in several loci (BFG095-BFG134, BFG095-RadalF5, RadalG11-RadalF5, BFG143-RadalF5), possibly due to a low frequency of null alleles or minor sibship sampling effects. Because no consistent pattern was observed for specific loci or populations, all six loci were used in the remaining analyses. The alleles per locus ranged from 3 to 16 (Table 2). Private alleles 3 to 20, and the average expected heterozygosities within populations ranged from 0.25 to 0.87 (Monsen & Blouin 2004). In addition, Monsen & Blouin (2004) indicated that the number of alleles per locus ranged from 1 to 11 per population, while the average expected heterozygosity per population ranged from 0.45 to 0.73, and the average observed heterozygosities per population ranged from 0.45 to 0.78 in Rana cascadae. The observed heterozygosity was found to be lower than in R. cascadae. This low observed heterozygosity could be related to Wahlund effects (reduced heterozygosity in populations due to the subpopulation structure), null alleles and deviation from the HWE (Freeland 2005). At least three loci were significantly deviant from HWE in all clusters of Anatolian mountain frogs (Table 3).
The different features of Anatolian mountain frogs make them vital for determining the genetic structure and gene flow among populations. The glacial ages also had an essential effect on the Anatolian mountain frogs, as researchers have reported that all the populations evolved from R. macrocnemis during the glacial and interglacial ages (Veith et al. 2003). This study is the first report on the population genetic structure of Anatolian mountain frogs using highly polymorphic, informative and high specificity (microsatellite) molecular markers.
The genetic variability in the Anatolian mountain frog group is similar to that of Rana cascadae with regard to the number of alleles per microsatellite. In montane Rana species (R. cascadae), the total number of alleles per microsatellite locus ranged from Table 3 Genetic diversity values for each group. Number of individuals (N), Number of effective alleles (Ne), observed heterozygosity (HO), Average expected heterozygosity (HE), Allellic richness (Ar), Average number of alleles per locus (A), inbreeding coefficient (FIS), significant deviation from the Hardy-Weinberg equilibrium (p) (the averages are not shown for some genetic parameters) Private alleles were found only in a single population among a broader collection of populations. There was also a negative correlation between private alleles and the migration rate (Slatkin 1985). Private alleles were observed in the Western (2) and Eastern Anatolia (7) clusters, but we found no private alleles for Central Taurus. Central Taurus is a corridor for the Anatolian mountain frogs between the East and West populations (Bilgin 2011). Therefore, a continuous gene flow over Central Taurus in the glacial ages could be the reason for the deficiency of private alleles. This feature (a high number of private alleles) also confirms that Eastern Anatolia is the source population of Anatolian mountain frogs.
The population expansion of the R. macrocnemis group started from the Caucasus to East Anatolia, and continued to inner and West Anatolia (Veith et al. 2003;Najibzadeh et al. 2021). East Anatolia had the highest values among the structure clusters regarding the number of alleles per locus, expected and observed heterozygosity, allelic richness and the number of effective alleles. Because East Anatolia constitutes the main population of Anatolian moun- The colours correspond to the three groups detected using STRUCTURE (Fig.1): East Anatolia group (orange), Central Taurus group (purple) and West Anatolia (blue). frog Hyla sarda, which is endemic to Corsica, Sardinia and their neighbouring islands (Bisconti et al. 2011). Neutral genetic divergence may also occur as the indirect result of local adaptation, in which natural selection against immigrants limits the gene flow (Wang & Bradburd 2014). Structure Harvester gave three clusters for the Anatolian mountain frog group. The Structure analysis also confirmed that the primary genetic variation appertained to the Eastern populations. The Structure clusters resembled the glacial expansion routes for the Anatolian mountain frogs. The East side of the Black Sea coast is one of the expansion routes for the ancestral populations of R. macrocnemis from the Caucasians. The individuals were distributed widely throughout the Western side of Anatolia after following this route. The allele sharing among the Bayburt, Artvin and Trabzon populations, and those of the Western side of Anatolia, supports this idea. Shared alleles and no significant population differentiation signify population intermixing in the glacial and interglacial periods for the Anatolian mountain frogs. A tree population sub-structure was given in the Structure analysis, in which all localities were grouped with the adjacent regions. This is further evidence that members with the same topography and climatic conditions are genetically relevant to each other.
The low proportion of the variance explained by the first two axes of the PCoA indicated that the planar graph might need to represent many variables more efficiently. Although the groups were slightly differentiated from each other, no strict distinction existed among the three clusters.
Anytime soon, human-mediated climate change and habitat loss may drastically reduce the level of biodiversity, with expected effects on many amphibian lineages (Souza et al. 2023). Researchers have stated that summer will witness the most significant temperature increases, affecting Türkiye's South East, Central Anatolia, Aegean and the Mediterranean regions (Bozoglu et al. 2019). According to the current results, the low population and genetic diversity in the Western population of frogs could be a disadvantage in the future, especially for those populations confronted with a problem like the global climate crisis. Although we have reflected on the genetic structuring of Anatolian mountain frogs in detail, we still need to know how local adaptation, phenotypic plasticity and genetic diversity interact to increase or decrease the fitness in Anatolian mountain frogs. Picariello et al. (1999Picariello et al. ( , 2016Picariello et al. ( and 2018 repeated their conclusion and declared identical S1 satellite DNA profiles for R. macrocnemis and R. tavasensis. Despite the high level of the mtDNA tain frogs, the genetic variability was higher than other in clusters. Furthermore, a hot environment and lowlands are unsuitable for Anatolian mountain frogs (Najibzadeh et al. 2017b), and when we compared the distributed region, the coldest area for the R. macrocnemis group was in East Anatolia. Young populations in the postglacial expansion areas exhibited low levels of genetic differentiation. Western Anatolia was the last area in Anatolia to be recolonised after the last glaciation events by Anatolian mountain frogs.
Many studies employing microsatellite data use the TPM mutation model to detect bottleneck events, because it is an intermediate of the IAM and SMM (Wang et al. 2010). According to the TPM mutation model, we observed a significant bottleneck event in only the East Anatolia cluster (Table 4). Populations from expansion areas often have overtones of recent bottlenecks, followed by solid demographic growth (Lessa et al. 2003;Hewitt 2004). While expanding the range into a new suitable habitat, these bottleneck events may have occurred during the retreat after the glacial periods. Frankham et al. (2002) said that a 0.25 level of inbreeding would be high enough to directly reduce fitness. We found a higher level for the inbreeding coefficient in the East Anatolia cluster. However, this result was not approved by the mode-shifted test. The smaller a bottleneck is, the more likely it will be detectable by the modeshift test (Luikart et al. 1998). We found one shifted cluster, for Central Anatolia, and the IAM mutation model supported this result. This result indicates that the Central Anatolia group had a bottleneck event in the recent past. Furthermore, East Anatolia has also undergone a genetic bottleneck event, but the effective population size has remained relatively high.
Shared genes indicate that the Plio-Pleistocene epoch is fit for Anatolian mountain frogs. Despite the big mountain chain existing among the members of this group, we found a moderate genetic differentiation value (0.078-0.105). The moderate F ST values could result from a lineage admixture during the last glacial period, where a fixed gene flow could have occurred between the diverging populations (Hey 2006). In the meantime, the F IS value was more significant than the zero sign of suffering inbreeding. This is also one of the signs that no current gene flow exists among the Anatolian mountain frogs.
Montane species are more likely to expand their range to the lowlands during glacial periods (Hewitt 2011;Stewart et al. 2010). A glacial phase demographic expansion is also considered to underlie the intraspecific genetic structure of the Tyrrhenian tree diversification rate, we could not find a significant genetic differentiation between R. macrocnemis and R. tavasensis. The mtDNA and microsatellite data would show similar divergence patterns in cases where an admixture between the lineages did not occur. In addition, if the differentiation between the mtDNA and microsatellites resulted from the admixture, the microsatellite and nuclear DNA sequence data sets would display a similar divergence pattern (Yang et al. 2016). Our result showed a similar divergence pattern to the study conducted by Ergül Kalaycı et al. (2017) based on nuclear DNA within the R. macrocnemis group. Therefore, we must carefully consider the species status of R. tavasensis under these circumstances.

Author Contributions
Research concept and design: T.E.K; Collection and/or assembly of the data: T.E.K; Data analysis and interpretation: T.E.K.; Writing the article: T.E.K.; Critical revision of the article: N.Ö., T.E.K; Final approval of the article: T.E.K., N.Ö.