Genetic Diversity and Population Dynamics of Leptobrachium leishanense (Anura: Megophryidae) as Determined by Tetranucleotide Microsatellite Markers Developed from Its Genome

Simple Summary More than 41 percent of amphibians evaluated by International Union for Conservation of Nature are threatened. It is vitally important to establish scientific and effective protection strategies for these organisms. Leishan Spiny Toad is endemic to China and it has a narrow distribution area. Long-term intentional human use and habitat destruction has caused the species to suffer. Here, we developed newly reliable and efficient molecular markers based on its genome to assess its genetic diversity and population history and provided support for conservation of this toad. Our results show that this toad still possesses high genetic diversity, but population decline may increase the possibility of inbreeding, which could work against persisting survival. Recovering the toad’s habitat and strengthening the publicity and education of wildlife protection can be helpful to its sustainability. Abstract Persisting declination of amphibians around the world has resulted in the public attaching importance to the conservation of their biodiversity. Genetic data can be greatly helpful in conservation planning and management, especially in species that are small in size and hard to observe. It is essential to perform genetic assessments for the conservation of Leptobrachium leishanense, an endangered toad and receiving secondary protection on the list of state-protected wildlife in China. However, current molecular markers with low reliability and efficiency hinder studies. Here, we sampled 120 adult toes from the population in the Leishan Mountain, 23 of which were used to develop tetranucleotide microsatellite markers based on one reference L. leishanense genome. After primer optimization, stability detection, and polymorphism detection, we obtained 12 satisfactory microsatellite loci. Then, we used these loci to evaluate the genetic diversity and population dynamics of the 120 individuals. Our results show that there is a low degree of inbreeding in the population, and it has a high genetic diversity. Recently, the population has not experienced population bottlenecks, and the estimated effective population size was 424.3. Accordingly, stabilizing genetic diversity will be key to population sustainability. Recovering its habitat and avoiding intentional human use will be useful for conservation of this species.


Introduction
Amphibians have long been declining on a global scale, and this trend will continue [1]. Furthermore, some amphibians face extinction or have become extinct [2]. There have been reports of massive declines in amphibians in many places, including areas where all species have been actively conserved [3,4]. Although there has been little consensus on the causes of this phenomenon [5], we recognize that amphibian populations are under serious threat and are in desperate need of conservation.
The Leishan Spiny Toad (Leptobrachium leishanense) is an endemic amphibian to China and is mainly restricted in Leishan county of Guizhou Province. This species inhabits broadleaf forests at elevations ranging from 1100-1800 m and breeds in slow-flowing streams via larval development [6]. The toad suffers from significant habitat loss and is often harvested for local consumption [7]. Thus, the population size has declined dramatically. It is listed as an endangered species on the International Union for Conservation of Nature (IUCN) Red List and receives secondary protection on the new list of state-protected wildlife in China. Formulating scientific conservation strategies is necessary for this species.
Genetic assessment is one of the aims of the conservation of biodiversity [8] and an important measure for amphibian population conservation [9][10][11]. Estimating genetic diversity and effective population size are the main goals of genetic assessments [12]. Genetic diversity reflects the adaptive potential of populations for environmental change [13]. When genetic diversity decreases, the extinction risk of populations increases [14]. Moreover, the levels of genetic diversity are related to population size [15]. It is a consensus that determining effective population size is more vital than measuring census size in populations [16]. In theory, small populations are susceptible to genetic depletion through drift and inbreeding, with adverse consequences for viability [17,18]. Therefore, effective population size can be used to assess the viability of populations.
As next-generation sequencing technologies offer new opportunities for conservation genetics [19], microsatellite markers with high mutation rates and genome-wide distributions reveal recent changes in genetic structure and demography critical for population management [20,21]. Although several studies using 10 dinucleotide microsatellites have shown that L. leishanense has high levels of genetic diversity and has not experienced recent bottleneck events [22][23][24], genetic assessments of this species are not nearly sufficient. In addition, dinucleotide microsatellites are considered less efficient and more unreliable than tetranucleotides because of their minimal PCR stutter [25]. Moreover, the traditional methods of microsatellite isolation and characterization are quite involved, costly, and time-consuming [26]. With the publication of a number of genomes, we can obtain sufficient numbers of different types of useful microsatellite loci more efficiently [27]. The genome sequencing project of L. leishanense has provided the opportunity to isolate and characterize microsatellites at the genomic level [28].
Here, we totally sampled 120 adult toes of L. leishanense from the population in the Leishan Moutain, and 23 of them were used to develop tetranucleotide microsatellite markers with polymorphisms based on one reference L. leishanense genome. After that, we analyzed the genetic diversity and population dynamics using the microsatellite loci we identified. The goals of this study were to (1) develop microsatellite loci with high reliability and efficiency, (2) evaluate the genetic diversity of the L. leishanense population, (3) detect if the population is experiencing a population bottleneck, (4) estimate the effective population size, and (5) provide molecular support for L. leishanense conservation planning.

Sampling
In 2012, 2013, 2014, 2015, and 2018, we collected 24 L. leishanense adults per year in Maoping village of Leishan County, Guizhou, China (Figure 1), sampled their toes, fixed the toes in anhydrous ethanol, and stored them in a −20 • C refrigerator. All individuals were released immediately after sampling. All experiments involving animals were approved by the Animal Ethics Committee of the School of Life Sciences, Central China Normal University (CCNU-IACUC-2019-008). We have complied with all relevant ethical regulations for animal testing and research. proved by the Animal Ethics Committee of the School of Life Sciences, Central China Normal University (CCNU-IACUC-2019-008). We have complied with all relevant ethical regulations for animal testing and research.

DNA Extraction and Primer Selection
DNA samples were extracted using the TIANamp DNA kit (Tiangen, Beijing, China) and stored at −20 °C. MicroSatellite identification tool (MISA-web, Gatersleben, Germany) [29] was used to obtain the simple sequence repeats (SSRs) of L. leishanense from its genome. Then, we randomly selected 87 tetrabase repeat microsatellite markers that were repeated more than 10 times and designed 87 pairs of primers according to the flanking sequences at both ends of each primer. With the extracted DNA as a template, we optimized the annealing temperature of the primers and reaction system. Each polymerase chain reaction (PCR) procedure was conducted in a 10 μL volume, in which the premix was 5 μL, each primer was 0.3 μL, template DNA was 0.6 μL, and ddH2O was 3.8 μL. The procedure was performed with initial denaturation at 94 °C for 5 min, followed by 35 cycles of denaturation at 94 °C for 30 s, annealing at temperature Ta for 30 s, extension at 72 °C for 45 s, and extension at 72 °C for 5 min. PCR products were detected by 1% agarose gel electrophoresis. By adjusting the Ta temperature, a product with a clear band was obtained. The Ta temperature corresponding to the product was used as the optimum temperature for PCR amplification. Under the optimal amplification conditions, we used the DNA of three different individuals to detect the stability of primers in different individuals and screened the primers that could be amplified stably.

Polymorphic Microsatellite Verification
The screened primers were used to synthesize 5′ upstream fluorescent primers (FAM, HEX and TEMED, compounded by Tiangen, Beijing, China). DNA amplification was performed on 23 individuals collected in 2012 and 2013 by PCR with fluorescent primers, and

DNA Extraction and Primer Selection
DNA samples were extracted using the TIANamp DNA kit (Tiangen, Beijing, China) and stored at −20 • C. MicroSatellite identification tool (MISA-web, Gatersleben, Germany) [29] was used to obtain the simple sequence repeats (SSRs) of L. leishanense from its genome. Then, we randomly selected 87 tetrabase repeat microsatellite markers that were repeated more than 10 times and designed 87 pairs of primers according to the flanking sequences at both ends of each primer. With the extracted DNA as a template, we optimized the annealing temperature of the primers and reaction system. Each polymerase chain reaction (PCR) procedure was conducted in a 10 µL volume, in which the premix was 5 µL, each primer was 0.3 µL, template DNA was 0.6 µL, and ddH 2 O was 3.8 µL. The procedure was performed with initial denaturation at 94 • C for 5 min, followed by 35 cycles of denaturation at 94 • C for 30 s, annealing at temperature Ta for 30 s, extension at 72 • C for 45 s, and extension at 72 • C for 5 min. PCR products were detected by 1% agarose gel electrophoresis. By adjusting the Ta temperature, a product with a clear band was obtained. The Ta temperature corresponding to the product was used as the optimum temperature for PCR amplification. Under the optimal amplification conditions, we used the DNA of three different individuals to detect the stability of primers in different individuals and screened the primers that could be amplified stably.

Polymorphic Microsatellite Verification
The screened primers were used to synthesize 5 upstream fluorescent primers (FAM, HEX and TEMED, compounded by Tiangen, Beijing, China). DNA amplification was performed on 23 individuals collected in 2012 and 2013 by PCR with fluorescent primers, and the amplified fluorescence PCR products were sent to Tsingke Biological Company, Beijing, China for SSR scanning and sequenced by an ABI 3730xl analyzer. Then, the products were genotyped and calculated, and the evaluation criterion of the polymorphisms was a PIC value higher than 0.5 [30]. We used Genemarker 1.3 software [31] to read the lengths of alleles, genotyped the microsatellite markers, and selected the sites with obvious poly-Animals 2021, 11, 3560 4 of 13 morphisms for the following analysis. The microsatellite genotyping data in Excel were transformed by using the Microsatellite Toolkit [32]. Cervus 3.0 software [33] was used to calculate the number of alleles (Na), polymorphism information content (PIC), expected heterozygosity (He), and observed heterozygosity (Ho). Micro-Checker 2.2.3 [34] was used to check large allele dropout of the microsatellite markers. GenePop 1.2 software [35] was used to detect the Hardy-Weinberg equilibrium (HWE) and linkage disequilibrium (LD) of the screening microsatellite markers with polymorphisms, and the Bonferroni correction was used for correction. The significance level was p < 0.05.

Genetic Diversity Analysis
A total of 120 DNA samples were amplified by PCR with the screening fluorescent primers described above. The PCR products were sent to Qingke Biological Company for SSR scanning. An ABI 3730xl analyzer was used for sequencing. Data were analyzed using GenAlEx 6.502 [36] to calculate the effective number of alleles (Ne), the mean relatedness of the individuals for every year, and the per year genetic differentiation coefficient (Fst), and Cervus 3.0 software was used again to calculate the values described above. Excel and Microsatellite Toolkit v3.1.1 software were used for preliminary genetic data statistics and data format conversion. FSTAT 2.9.3.2 software [37] was used to calculate allelic richness (Ar), allelic diversity (Hs), and the inbreeding coefficient (Fis).

Population Bottleneck Identification
Bottleneck 1.2.02 software [38] was used to test whether the population had experienced population bottlenecks. Sign and Wilcoxon methods were used to test mutations through three mutation models: the infinite allele model (IAM), the stepwise mutation model (SMM), and the two-phased model of mutation (TPM). TPM was set to 95% SMM, with a variance of 30 and 1,000 iterations.

Effective Population Size Calculation
NeEstimator 2.1 [39] was used to calculate the effective population size by selecting the random mating model, and the confidence interval was 95%.

Distribution of SSR in Genome of L. leishanense
A total of 1,454,145 microsatellite markers were obtained from the genome of L. leishanense. Monobase repeat microsatellite markers and dibase repeat microsatellite markers were the most common among all microsatellite markers. There were 874,773 monobase repeat microsatellite markers and 263,927 dibase repeat microsatellite markers, accounting for 60.16% and 18.15% of the total number of microsatellite markers, respectively, followed by 71,167 tribase repeat microsatellite markers and 23,332 tetrabase repeat microsatellite markers, accounting for 4.89% and 1.60% of the total number of microsatellite markers, respectively. The number of pentabase repeat microsatellite markers and hexabase repeat microsatellite markers was the lowest, with 909 pentabase repeat microsatellite markers and 844 hexabase repeat microsatellite markers, accounting for only 0.12% of the total microsatellite markers (Figure 2).

Polymorphism Microsatellite Loci
Eighty-seven pairs of primers randomly chosen from 23,332 tetrabase repeat microsatellite marker. After primer optimization, 64 pairs of primers were successfully amplified. Then stability detection was used, and 46 pairs of primers were obtained. Employing polymorphism detection, we obtained 12 satisfactory microsatellite loci. The Na of these loci ranged from 6-16. PIC values ranged from 0.537-0.904. Ho and He were between 0.609-0.913 and between 0.622-0.931, respectively (Table 1). All 12 loci were not significant with regard to LD (p > 0.05), and there were no loci deviated from HWE

Polymorphism Microsatellite Loci
Eighty-seven pairs of primers randomly chosen from 23,332 tetrabase repeat microsatellite marker. After primer optimization, 64 pairs of primers were successfully amplified. Then stability detection was used, and 46 pairs of primers were obtained. Employing polymorphism detection, we obtained 12 satisfactory microsatellite loci. The Na of these loci ranged from 6-16. PIC values ranged from 0.537-0.904. Ho and He were between 0.609-0.913 and between 0.622-0.931, respectively (Table 1). All 12 loci were not significant with regard to LD (p > 0.05), and there were no loci deviated from HWE (p > 0.05). According to Micro-Checker 2.2.3, there was no large allele dropout of these microsatellite markers and no scoring error caused by the shadow peak.

Population Genetic Diversity
Using the 12 loci screened above, all 120 individuals were used to study the genetic diversity of the population. HWE detection of the L. leishanense population was performed. After Bonferroni correction, the significance level was p < 0.0042. The results are shown in Table 2. When we used all 120 samples for testing, loci LEA23, LEA7, LEA47, LEA2, and LEA53 deviated from HWE significantly, with Fis values as 0.179, 0.112, 0.119, 0.230, and 0.134, respectively (Table 3). When we separated the samples into each year for testing, locus LEA20 deviated from HWE significantly in 2013, with a Fis value of 0.215. Locus Then, we calculated pairwise year Fst valus in L. leishanense ( Table 3). None of these values is greater than 0.05, suggesting the genetic differentiation between these years is negligible [40]. Further, we calculated the mean relatedness of the individuals for every year (Figure 3). In 2013, 2014, and 2015, mean pairwise relatedness within groups was significantly greater than zero, indicating the samples we collected in these three years have relatively close relationships.    Next, we calculated Na, Ne, PIC, Ho, He, Ar, Hs, and Fis of the population ( Table 4). The results indicated that the genetic diversity of the population was still high. The positive value of Fis and that of Ho was lower than that of He, suggesting that there was a low degree of inbreeding in the population.  Next, we calculated Na, Ne, PIC, Ho, He, Ar, Hs, and Fis of the population ( Table 4). The results indicated that the genetic diversity of the population was still high. The positive value of Fis and that of Ho was lower than that of He, suggesting that there was a low degree of inbreeding in the population.

Population Bottleneck
The average expected heterozygosity (Heq) of the population in the IAM, SMM, and TPM models was calculated (Table 5). In the IAM model, there were 11 sites where He was significantly higher than Heq (p < 0.05), among which LEA22, LEA25, LEA20, LEA35, LEA14, LEA23, and LEA24 were extremely significantly higher than Heq (p < 0.01). In the TPM and SMM models, only He at LEA5 was significantly higher than Heq, showing heterozygote surplus. The sign test and Wilcoxon test were used to detect the heterozygosity surplus of the population under the three models of IAM, TPM, and SMM ( Table 6). The mutationdrift balance of the population was detected under the IAM model, both of the sign and Wilcoxon tests showed significant deviations from the mutation-drift balance of the population. Under the TPM and SMM models, both the results of the sign test and Wilcoxon test showed that the population did not deviate from mutation-drift equilibrium.   The analysis of allele frequency distribution in the L. leishanense population showed that the allele frequency was mainly concentrated between 0.0-0.1, which was approximately 84.52% of the total allele frequency. Alleles with a frequency of 0.1-0.2 accounted for 12.69% of the total allele frequency. The proportion of the frequency distribution interval of 0.2-0.3 was 1.80%, while the allele proportions of the frequency distribution intervals of 0.3-0.4 and 0.4-0.5 were 0.5%. The allele frequency showed a typical "L" type distribution (Figure 4), suggesting that the population has not recently experienced a bottleneck effect.

Discussion
We isolated and characterized 12 tetrabase repeat microsatellite markers with polymorphisms from one reference genome of L. leishanense. Then, we used these loci to study the genetic diversity and population dynamics of this species. We found that the genetic diversity of the population was high and that there was a low degree of inbreeding in the population. Moreover, the population has not recently experienced bottleneck effects, and the estimated effective population size is 424.3.

Tetranucleotide Microsatellite Markers
Although the results above are similar to those of Zhang's research [24], which used 10 dibase repeat microsatellite markers, the 12 tetranucleotide microsatellite markers we developed are more polymorphic and suitable for genetic diversity research. During PCR amplification, a biological phenomenon called stutter is generated due to chain slippage, resulting in typing errors, and the stutter product has one or more fewer duplicates than the real allele product [41,42]. In general, tetranucleotide repeats tend to stutter less than trinucleotide and dinucleotide repeats and are much more accurate and reliable [43,44]. Therefore, in different types of microsatellite systems, tetrabase repeat microsatellite markers are more common than dibase or tribase markers. Moreover, the PIC values of all 12 loci were higher than 0.5, suggesting that the loci we developed had higher polymorphism. Stable and reliable microsatellite markers are a necessary prerequisite for population estimation in the wild [45]. Thus, after primer optimization, stability detection, and polymorphism detection, we finally obtained 12 satisfactory tetranucleotide microsatellite loci.

Discussion
We isolated and characterized 12 tetrabase repeat microsatellite markers with polymorphisms from one reference genome of L. leishanense. Then, we used these loci to study the genetic diversity and population dynamics of this species. We found that the genetic diversity of the population was high and that there was a low degree of inbreeding in the population. Moreover, the population has not recently experienced bottleneck effects, and the estimated effective population size is 424.3.

Tetranucleotide Microsatellite Markers
Although the results above are similar to those of Zhang's research [24], which used 10 dibase repeat microsatellite markers, the 12 tetranucleotide microsatellite markers we developed are more polymorphic and suitable for genetic diversity research. During PCR amplification, a biological phenomenon called stutter is generated due to chain slippage, resulting in typing errors, and the stutter product has one or more fewer duplicates than the real allele product [41,42]. In general, tetranucleotide repeats tend to stutter less than trinucleotide and dinucleotide repeats and are much more accurate and reliable [43,44]. Therefore, in different types of microsatellite systems, tetrabase repeat microsatellite markers are more common than dibase or tribase markers. Moreover, the PIC values of all 12 loci were higher than 0.5, suggesting that the loci we developed had higher polymorphism. Stable and reliable microsatellite markers are a necessary prerequisite for population estimation in the wild [45]. Thus, after primer optimization, stability detection, and polymorphism detection, we finally obtained 12 satisfactory tetranucleotide microsatellite loci.

Genetic Diversity
We speculated that several loci deviated from HWE (Table 2) mainly caused by sampling from the same family ( Figure 3). The sharp decline of the population size has increased the possibility in sampling individuals of same family. The genetic diversity of a population is a long-term process, the population of L. leishanense does not have significant genetic differentiation among these five years (Table 3); accordingly, we considered that these deviating loci were still effective in estimating population genetic diversity. Then, a series of indices were used to measure the genetic diversity of the toad, including Na, Ar, Ne, PIC, Ho, He, and Hs. According to our results, the toad still has high genetic diversity. Threatened species usually have small or declining populations and are prone to loss of genetic diversity due to inbreeding or genetic drift [14]. As an endangered and narrowly distributed toad, the population shows the opposite result. Several studies investigating endangered or narrowly distributed species have obtained similar results [45][46][47][48][49], indicating that endangered species or species with a narrow distribution may also have high levels of genetic diversity. When the earth was in an ice age, some areas with a stable ecological environment became the refuge of organisms, and the populations living in the refuge survived and accumulated rich genetic diversity [50]. The Leishan Spiny Toad is a relatively primitive species, and its formation dates back to the Miocene [22]. The toad survived by staying on Leigong Mountain and retained rich genetic diversity when the ice age came. In addition, two additional distribution sites were found by Zheng et al. [51], suggesting that the toad is not strictly a narrowly distributed species. We may have underestimated the genetic diversity of the species.
Although the Fis value of the population is on the low degree, this does not mean that there is no inbreeding between the individuals in the population. According to our year-byyear field work, its population size is declining. This undoubtedly increases the possibility of its inbreeding. Inbreeding has a negative effect on the fitness of the population, including fertility and viability [52], which is not conducive to the long-term development of the population. We could not find more obvious molecular evidence of inbreeding, possibly due to our restricted sampling size and the relatively high number of alleles found. As with high number of alleles, the probability of obtaining homozygote hgenotypes in one locus is very low. Thus, it will influence our detection of inbreeding.

Population Dynamics
Combining the results of model simulation with allele frequency distribution, we find that the population has not recently experienced a bottleneck effect. We tested three models, and IAM was significant both in the sign test and the Wilcoxon test. Both SMM and TPM were not significant in the sign test and Wilcoxon test (Table 5). IAM assumes that there is only one mutation of an allele in a population, and each mutation produces a new allele, which is generally used in isozyme or DNA sequencing data. SMM supposes that alleles can mutate upward or downward into new alleles. TPM is the synthesis of the previous two models, and the probability of occurrence of two kinds of mutations can be determined. The principle of allelic mutation in microsatellite data is the increase or absence of repeating units, which is represented by the change in sequence length. Some studies believe that TPM is more suitable for microsatellite data [53]. Therefore, we accept the result of TPM that there is no significant excess heterozygosity in this population. That is, the rate of heterozygosity decrease is approximately the same as the rate of allele loss in L. leishanense, indicating that the population has not recently experienced a bottleneck.
However, the ability to detect the population bottleneck based on heterozygosity is limited, and the number of alleles is more sensitive to population fluctuation, so it is more reliable to analyze the distribution of allele rates in the case of heterozygous residues to determine whether the population has experienced the bottleneck effect [54]. To enhance the adaptability to environmental changes, species tend to accumulate many rare alleles with low frequency. Therefore, the frequency distribution of alleles in mutation-drift equilibrium shows an "L" shape. If the species recently experienced a genetic bottleneck, the distribution of alleles with low frequency (0.0-0.1) will change to a mid-frequency distribution (0.1-0.2); thus, the allele frequency distribution will deviate from "L" [55]. In this study, the frequencies of the allele were generally in a typical "L" (Figure 3), suggesting that the population did not experience bottleneck effects.
The LD distance between microsatellite markers can be used to estimate effective population size, and this method has been applied to mammals, fish, amphibians, and other animals [56]. Effective population size is a valuable method in population conservation and management research. Maintaining an effective population of sufficient size is a key factor to maintain the rich genetic diversity of the population. Based on the microsatellite loci we developed, the estimated effective population size of L. leishanense is 424.3. Nei et al. [57] deemed that the population size should be 4-10 times of the effective population size to maintain the stability of population genetic diversity. Therefore, to maintain the stability of the population, the number of Leishan Spiny Toads should be 1697.2-4243. However, while LD information is used to estimate the effective population size, the accuracy of the results is significantly correlated with the sample size [58]. More samples may be needed to obtain more reliable and accurate results in L. leishanense.

Conclusions
Our study has provided 12 reliable tetranucleotide microsatellite loci with polymorphisms, enriching the information regarding the genetic diversity and population dynamics of L. leishanense. Although the genetic diversity is still high based on our results, a low degree of inbreeding indicates that the population is declining. Avoiding habitat fragmentation and intentional human use will be key to the conservation of this species. Furthermore, recovering the streams and woodlands where the species once existed abundantly will also help to stabilize its genetic diversity.

Data Availability Statement:
The datasets generated and analyzed during the current study are available from the corresponding author on reasonable request.