Genetic Diversity and Population Structure of Llamas (Lama glama) from the Camelid Germplasm Bank—Quimsachata

Llamas (Lama glama) are invaluable resources of Peru. Despite their importance, their population is decreasing. The Camelid Germplasm Bank—Quimsachata was created as a guardian of this South American camelid (SAC) species and established a bank of llamas from their two types, Ch’aku and Q’ara. However, these populations need to present high genetic diversity to be considered suitable conservation stocks. Thus, in the present study, 13 microsatellites specific for the SAC were used to assess the current genetic variability and differentiation of the llama population from the Bank. The global population showed high genetic diversity with a total of 157 different alleles, with an average of 12.08 alleles per microsatellite, an expected and observed heterozygosity of 0.758 and 0.707, respectively, and an average polymorphic information content (PIC) of 0.723. Although considered as two different breeds and managed separately, the genetic differentiation between Ch’aku and Q’ara was low (FST = 0.01). Accordingly, the gene flow value was high (Nm = 30.5). Overall, our results indicate the existence of high genetic variation among individuals, and thus, this llama population could be considered a suitable genetic stock for their conservation and for sustainability programs. Additionally, the 13 microsatellites can be used to study other Peruvian llama populations and monitor the genetic variability of llamas from the Camelid Germplasm Bank—Quimsachata.


Introduction
Llamas (Lama glama) are the largest South American camelids (SACs) and the best adapted to a wide range of environmental conditions [1]. They can be mainly found in the Andean region between 2300 and 4000 m above sea level in Peru, Bolivia, Ecuador, northwest Argentina, and central Chile [2]. After Bolivia, Peru is the country with the second-largest population of llamas in the world [3], being primarily found in the Department of Puno, with almost 35% of the total Peruvian llama population [4]. For the local economy, llamas are seen as multipurpose animals; given their low-fat and low-cholesterol but high-protein content, their meat is consumed [2], and occasionally their intestines are used to make string and drums, while their excrement is used as fuel. Moreover, they are pack animals and their fiber is often used for clothing [5]. Two main types of llamas are recognized:

Sample Collection
Blood samples from 251 adult llamas of types Ch'aku (n = 92, 67 females and 25 males) and Q'ara (n = 159, 109 females and 50 males) were collected at the Camelid Germplasm Bank-Quimsachata, located in the district of Santa Lucia in the Department of Puno, Peru ( Figure 1). The protocol used for the blood sample collection agreed with the requirements of the National Law No. 30407 "Ley de Protección y Bienestar Animal (Animal Protection and Welfare Law)". To measure the genetic differentiation between both phenotypes, Ch'aku and Q'ara llamas were considered as two different subpopulations. Llamas were selected after analyzing their pedigree records, and only unrelated animals were sampled. The difference in the number of animals between the two populations and in the number of males and females that were enrolled in the study were due to the fact that the Camelid Germplasm Bank-Quimsachata has more Q'ara than Ch'aku phenotypes and a male to female ratio of approximately 30 to 70. Thus, the sampling was designed to keep this ratio and to have a representative sample of the population. Total genomic DNA was extracted using a standard phenol-chloroform and ethanol precipitation protocol [39]. The DNA pellet was resuspended in TE buffer and stored at −20 • C.
Genomic DNA was amplified by polymerase chain reaction (PCR) using a Mastercycler Pro S (Eppendorf, USA). Microsatellites with similar PCR conditions were co-amplified using a multiplex PCR, and the fluorescent labeling of the forward primers allowed for the design of multiloading panels. Each PCR run was performed in a total volume of 10 μL, containing 50 ng of DNA. Reaction mixtures contained 1x PCR buffer, 2 mM of MgCl2, 0.2 μM of each primer, 0.2 mM of dNTPs, and 0.5U of Taq polymerase. The amplification conditions included an initial denaturation step of 95 °C for 5 min, followed by 35 cycles of 95 °C for 30 s, 90 s at 56 °C or 58 °C, 72 °C for 1 min, and a final extension at 72 °C for 30 min. The PCR products were separated by capillary electrophoresis in an automatic ABI Prism 3130XL Genetic Analyzer ® (Applied Biosystems, Foster City, CA, USA). Genotyping was performed using GeneMapper v.4.0 software (Applied Biosystems).

Statistical Analysis
The most used genetic diversity parameters, such as the number of alleles per microsatellite, allelic frequencies, mean number of alleles per microsatellite, private alleles, expected and observed heterozygosity (He and Ho, respectively), and polymorphic information content (PIC), were calculated for each of the 13 microsatellite markers, using the CERVUS 3.0.3 [42] and GENETIX 4.0.5 [43] software. Possible deviations from the Hardy-Weinberg equilibrium (HWE), either due to an excess or to a deficit of heterozygous in the total population and within subpopulations, were estimated using Fisher's exact test implemented in the GENEPOP 4.0.11 software [44]. Furthermore,
Genomic DNA was amplified by polymerase chain reaction (PCR) using a Mastercycler Pro S (Eppendorf, Hauppauge, NY, USA). Microsatellites with similar PCR conditions were co-amplified using a multiplex PCR, and the fluorescent labeling of the forward primers allowed for the design of multiloading panels. Each PCR run was performed in a total volume of 10 µL, containing 50 ng of DNA. Reaction mixtures contained 1x PCR buffer, 2 mM of MgCl 2 , 0.2 µM of each primer, 0.2 mM of dNTPs, and 0.5U of Taq polymerase. The amplification conditions included an initial denaturation step of 95 • C for 5 min, followed by 35 cycles of 95 • C for 30 s, 90 s at 56 • C or 58 • C, 72 • C for 1 min, and a final extension at 72 • C for 30 min. The PCR products were separated by capillary electrophoresis in an automatic ABI Prism 3130XL Genetic Analyzer ® (Applied Biosystems, Foster City, CA, USA). Genotyping was performed using GeneMapper v.4.0 software (Applied Biosystems).

Statistical Analysis
The most used genetic diversity parameters, such as the number of alleles per microsatellite, allelic frequencies, mean number of alleles per microsatellite, private alleles, expected and observed heterozygosity (He and Ho, respectively), and polymorphic information content (PIC), were calculated for each of the 13 microsatellite markers, using the CERVUS 3.0.3 [42] and GENETIX 4.0.5 [43] software. Possible deviations from the Hardy-Weinberg equilibrium (HWE), either due to an excess or to a deficit of heterozygous in the total population and within subpopulations, were estimated using Fisher's exact test implemented in the GENEPOP 4.0.11 software [44]. Furthermore, to guarantee the quality of the results, the null allele frequencies were also calculated using Micro-Checker 2.2.3 [45]. The level of genetic differentiation among individuals, within and between subpopulations, was calculated by the analysis of molecular variance (AMOVA) test using the ARLEQUIN 3.1 software [46]. The extent of genetic differentiation between the two llama subpopulations, Ch'aku and Q'ara, was quantified using the F-statistics (F IS , F IT , and F ST ; [47]) using GENEPOP 4.0.11 and corroborated with FSTAT 2.9.3.2 [48]. The effective number of migrants per generation (Nm, [19]) was estimated by GENEPOP 4.0.11.
The genetic structure was determined through a grouping analysis of the individuals in a different number of inferred clusters (K) using an analysis based on the 'admixture' ancestry model implemented in the STRUCTURE 2.3 software [49]. The burn-in period was set to 50,000 followed by 500,000 Markov chain Monte Carlo (MCMC) iterations. Independent runs of K were performed from 1 to 7 clusters and were repeated 4 times to check the consistency of the results. Finally, a factorial correspondence analysis was performed with GENETIX 3.0.3 to further assess the genetic relationships between the llama types, describing the association of qualitative variables in which each individual is represented just once for the value of each modality (microsatellites) and variable (alleles per microsatellite).

Genetic Diversity Assessment
The 13 microsatellites used in this study were polymorphic and revealed a high level of genetic diversity in the population of llamas from the Camelid Germplasm Bank-Quimsachata. The number of alleles, He, Ho, and PIC values for each microsatellite for Ch'aku and Q'ara as separated populations and as a global population are shown in Table 1.
In total, 157 alleles were found across the 251 individuals, with an average of 12.08 alleles per microsatellite (ranging from 8 for LCA54, LCA83, and LCA85 up to 19 for YWLL08 and YWLL59). The He value was higher than the Ho value, with an average of 0.707 in the global population. Additionally, the He and Ho values were always higher in the Ch'aku subpopulation compared to the Q'ara. The PIC value ranged from 0.517 (LCA54) to 0.883 (YWLL08).
For the global population, 7 out of the 13 microsatellites departed from the HWE (p < 0.05). The inbreeding coefficient (F IS ) was estimated for each microsatellite for the global population and for the Ch'aku and Q'ara subpopulations, and they were positive and significantly different from zero (p > 0.05) ( Table 2). Moreover, we observed a deficit of heterozygotes and 7.1% more homozygotes than would be expected under the HWE ( Table 2). The existence of null alleles was assessed with Micro-Checker, and only the microsatellites LCA82A, YWLL59A, LCA85A, and GLM4 showed signs of having them (Table S2). On the other hand, a total of 35 private alleles were observed (15 for Ch'aku and 20 for Q'ara subpopulations), but their frequencies were very low (<0.05) ( Table 3).

Genetic Differentiation between Ch'aku and Q'ara
The genetic differentiation coefficient value (F ST ) between the Ch'aku and Q'ara types showed very low genetic differentiation (F ST = 0.01) ( Table 4), whereas a high gene flow value was observed (Nm = 30.9 migrants per generation). Additionally, the analysis of molecular variance indicated that the highest variance was due to variation within the populations (93.95%), while only 1.02% variance was observed between the two subpopulations (Table S3). The F ST value calculated by Arlequin 3.1 was 0.0102 and coincided with the value given by FSTAT 2.9.3.2 [48]. Bayesian analysis carried out by the STRUCTURE software [49], using seven independent runs (K = 1, 2, 3, 4, 5, 6, 7) and each one repeated four times, did not show any evidence of genetic differentiation or population subdivisions. The highest likelihood was obtained when K = 2 and the individuals of both phenotypes were assigned to two clusters, but there was not a clear separation between the two llama subpopulations (Figure 2). Thus, there was no evidence of population structure in the Camelid Germplasm Bank-Quimsachata. Additionally, the correspondence factorial analysis showed a graphic representation of the genetic relationship between the Ch'aku and Q'ara subpopulations ( Figure 3). This analysis indicated once again the low genetic differentiation between the subpopulations, given that both overlap with each other and do not form clear independent groups. Genes 2020, 11, x FOR PEER REVIEW 9 of 13 individuals of both phenotypes were assigned to two clusters, but there was not a clear separation between the two llama subpopulations (Figure 2). Thus, there was no evidence of population structure in the Camelid Germplasm Bank-Quimsachata. Additionally, the correspondence factorial analysis showed a graphic representation of the genetic relationship between the Ch'aku and Q'ara subpopulations (Figure 3). This analysis indicated once again the low genetic differentiation between the subpopulations, given that both overlap with each other and do not form clear independent groups.

Discussion
There has been an upsurge in attempts to conserve natural resources, since it is well known that the loss of genetic variability diminishes the ability to recover endangered species [50] and decreases the chance to improve the performance of animals involved in production systems [51]. In this scenario, the Camelid Germplasm Bank-Quimsachata plays a crucial role in conserving a natural, economic, cultural, and historical resource of the SAC populations in Peru. The population of llamas in this Bank must have high genetic diversity to be considered an appropriate genetic stock that could contribute to ensuring its conservation, aid in implementing future strategies to face the loss of diversity, and increase the viability and productivity of llama populations from other regions of Peru. This study is the first to report on the genetic diversity and population structure of the llama population in the Camelid Germplasm Bank-Quimsachata using microsatellite markers. Genes 2020, 11, x FOR PEER REVIEW 9 of 13 individuals of both phenotypes were assigned to two clusters, but there was not a clear separation between the two llama subpopulations (Figure 2). Thus, there was no evidence of population structure in the Camelid Germplasm Bank-Quimsachata. Additionally, the correspondence factorial analysis showed a graphic representation of the genetic relationship between the Ch'aku and Q'ara subpopulations (Figure 3). This analysis indicated once again the low genetic differentiation between the subpopulations, given that both overlap with each other and do not form clear independent groups.

Discussion
There has been an upsurge in attempts to conserve natural resources, since it is well known that the loss of genetic variability diminishes the ability to recover endangered species [50] and decreases the chance to improve the performance of animals involved in production systems [51]. In this scenario, the Camelid Germplasm Bank-Quimsachata plays a crucial role in conserving a natural, economic, cultural, and historical resource of the SAC populations in Peru. The population of llamas in this Bank must have high genetic diversity to be considered an appropriate genetic stock that could contribute to ensuring its conservation, aid in implementing future strategies to face the loss of diversity, and increase the viability and productivity of llama populations from other regions of Peru. This study is the first to report on the genetic diversity and population structure of the llama population in the Camelid Germplasm Bank-Quimsachata using microsatellite markers.

Discussion
There has been an upsurge in attempts to conserve natural resources, since it is well known that the loss of genetic variability diminishes the ability to recover endangered species [50] and decreases the chance to improve the performance of animals involved in production systems [51]. In this scenario, the Camelid Germplasm Bank-Quimsachata plays a crucial role in conserving a natural, economic, cultural, and historical resource of the SAC populations in Peru. The population of llamas in this Bank must have high genetic diversity to be considered an appropriate genetic stock that could contribute to ensuring its conservation, aid in implementing future strategies to face the loss of diversity, and increase the viability and productivity of llama populations from other regions of Peru. This study is the first to report on the genetic diversity and population structure of the llama population in the Camelid Germplasm Bank-Quimsachata using microsatellite markers.
A representative sample of llama individuals (n = 251) was analyzed by using 13 microsatellites specific to the SAC. The results showed a high level of genetic diversity in the population of llamas, with an average of 12.08 alleles per microsatellite and an expected heterozygosity of 0.758. The Camelid Germplasm Bank-Quimsachata keeps the birth and mating records of each individual and that information is used to avoid mating between close relatives [13]. Thus, our results might indicate that this management strategy has contributed to maintaining high genetic variability, and it will be complemented with the routine assessment of the genetic variability by means of molecular markers such as microsatellites or single-nucleotide polymorphisms (SNPs). Furthermore, our results are similar to those from Argentinean (7.33-8.33 mean alleles per microsatellite and He of 0.47-0.9; [1,52]) and Bolivian (12.04 alleles per microsatellite and He = 0.68) llama populations [2]. Likewise, in those studies, the He was always higher than the Ho, and a higher value of He was reported in the Quimsachata llama population compared to the Bolivian llama population. Nonetheless, it is important to mention that they used different sets of microsatellite markers (many of which were individually used in this study), and in the case of Barreta et al. (2013), they were not necessarily specific for SACs. Thus, to make a more adequate comparison, we would ultimately have to use the same set of microsatellites used in those studies.
The F IS media value for the total population (0.063) and the microsatellites showing deviations from the HWE were explained by heterozygote deficiency. The excess of homozygotes in a domestic population indicates loss of genetic variability and could be explained by a lack of random mating, which occurs during the artificial selection of herds [1], population subdivisions (Wahlund effect), gene flow, or the existence of null alleles [53]. Indeed, we observed the presence of null alleles in four microsatellite markers (LCA82A, YWLL59A, LCA85A, and GLM4; Table S2), which could have potentially increased the observed homozygosity value. Additionally, other possible explanations for homozygote excess could be due to the evolutionary history of llamas, consisting of polygynous behavior in which herds contain an α male who controls the access of other males to its territory and expels young males while retaining the female individuals. This behavior of a young male's exclusion is still present in managed populations [1,54]. A departure from the HWE could be due not only to a genotyping error (null alleles) in some of the microsatellites but also to the repetitive mating of individuals within the same herd and to the low quantity of breeding males.
Regarding the llama population structure from the Camelid Germplasm Bank-Quimsachata, although we found a large number of private alleles in the Ch'aku and Q'ara subpopulations (15 and 20, respectively; Table 3), the AMOVA analysis showed that most of the variation came from the variability within the Ch'aku and Q'ara subpopulations (93.95%), whereas only 1.02% was due to the variability between the subpopulations (Table S3). Furthermore, their genetic differentiation index (F ST ) was 0.0102 along with a high gene flow value (Nm = 30.5). These values, according to [47], are indicators of low genetic differentiation between the two subpopulations and a weak genetic structure. These results are confirmed by the cluster analysis and the correspondence factorial analysis ( Figure 3) and are similar to the values estimated by [2], among and within regional groups of Bolivian llamas.
Importantly, the breeding management of the Bank has been highly controlled since its beginnings, that is, all Ch'aku and Q'ara males are reared together but kept separately from Ch'aku and Q'ara females (who are also reared together) and the males and females of the same phenotype are only brought together for breeding. Therefore, the low genetic structure is most likely a result of their own ancestral domestication process, which involved frequent exchanges of reproductive males and, hence, crossbreeding between both types. For future studies to elucidate a potential genetic differentiation between these phenotypes, we propose the use of genes related to observed phenotypic traits (e.g., the diameter of fiber), which could help to identify allelic variants related to the studied phenotypes [2], or the use of SNPs, which cover the whole llama genome [12].

Conclusions
This study is the first assessment of the genetic diversity and structure of the Peruvian llama (L. glama) population from the Camelid Germplasm Bank-Quimsachata utilizing microsatellites. The set of microsatellites used was highly polymorphic and, hence, can be utilized to track the genetic variability of these animals to avoid a reduction in the effective population size. Overall, our results show that the llama population of the Camelid Germplasm Bank-Quimsachata presents a high level of genetic diversity; thus, it can be considered as an adequate stock for the conservation of this natural resource. However, it would be important to identify other Peruvian llama populations that could help to increase the diversity of the Bank and enable a higher representation of Peruvian llamas.
On the other hand, although bred separately and managed as two distinct subpopulations, we observed low genetic differentiation between the Ch'aku and Q'ara phenotypes. Therefore, additional analysis with genes related to phenotypic observed differences should be carried out to assess potential genetic differentiation between these two phenotypes.
Supplementary Materials: The following are available online at http://www.mdpi.com/2073-4425/11/5/541/s1, Table S1: Primer sets used for genotyping the microsatellites, Table S2: Existence and Frequencies of null alleles found in each microsatellite marker analyzed with the software Micro Checker 2.2.3 and CERVUS 3.0.3, respectively Table S3: Analysis of molecular variance (AMOVA) between the subpopulations of llamas Ch'aku and Q'ara.