Genomic characterization is an important step toward implementing efficient breeding and conservation programs for endangered local breeds. Animal genetic diversity is becoming critical for food security and rural development especially in light of changing conditions such as climate change, new or reviving human or animal disease threats, as well as changes in market and societal needs. Local breeds that are well adapted to local and increasingly drier/warmer conditions should play a more prominent role in the livestock production and food security in both developed and developing countries. Better genetic characterization of local goat breeds is essential for more efficient genetic improvement programs targeting adaptive traits and conservation management strategies.
In this study, the genetic diversity analyses of the Slovenian Drežnica goat, five Austrian goat breeds and one from South Tyrol in north-eastern Italy was performed for the first time. A comparison was done with already published alpine goat breeds, most of which are important and rare local breeds. Apart from the importance of these breeds in food production, these breeds also represent cultural heritage of local societies, but they are often threatened by replacement or crossbreeding with more productive cosmopolitan commercial breeds. This process often leads to a significant decrease in population size, which in turn results in inbreeding depression and lower performance, providing additional reasons for their replacement. The worst-case scenario is extinction of local breeds and consequently the loss of key traits for the survival and management of flocks especially in extensive production systems. Traits such as resistance to local diseases, resilience, adaptation to poor forage and water resources, homing and gregarious behavior are crucial in the harsh alpine environment. Many local alpine goat breeds have retained such characteristics, which can help them overcome challenges related to negative effects of climate change. Specifically, the average temperature in the alpine region has risen in recent decades by nearly 2 °C, which is almost twice as large as the average global increase (42). Genetic characterization studies of rare local breeds such as those studied in here can provide proof of their genetic identity and a basis for maintaining genetic variation, improving their performance and conservation strategies.
Overall, our genetic diversity analysis of goats from the alpine area revealed that they still retained relatively high levels of genetic variability, as was also found in other studies with SNP arrays (7–9, 16). However, we noted some differences that warrant discussion. Since the goat SNP chip used here was designed based on genomic data from cosmopolitan breeds (8), it could be biased against rare and more diverse local breeds, an issue called ascertainment bias (43). To mitigate such bias, we used 4-SNP haplotype blocks as marker units rather than single-SNP alleles. Compared with breeds from the neighboring alpine area, the Drežnica goat tends to have a larger total number of observed alleles, mean number of alleles per block and number of private alleles. These results are encouraging, as the breed went through bottlenecks in the past and the population size is low today. The Drežnica goat has the second highest number of population-specific alleles, supporting its distinct genetic identity as a breed. The high number of private alleles could explain its excellent adaptability to the adverse climate/forage conditions in the Alps and indicate that frequent admixture with other alpine breeds was not common. This is likely why the Drežnica goat breed contributes 6.8% to the synthetic pool of alpine breeds with the maximal total number of alleles (Fig. 3B), which is among the top contributions compared to other alpine breeds. The Drežnica goat also contributed positively to total allelic diversity of the Alp2Step dataset (Fig. 3A). After removing the Appenzell, Toggenburg, Tauern Pied, Valais, and Drežnica goats individually from Alp2Step, the recalculations showed that these breeds accounted for the highest percentage of allelic diversity among the populations. This is in line with the results of the differentiation analysis for these five breeds, whereby they exhibited the highest pairwise distance values from other alpine breeds (Fig. 5). Therefore, based on various genetic diversity analyses, we can conclude that the alpine breeds, including local breeds with small population sizes, generally retained appreciably high levels of genetic variability with a few breeds excelling, such as the Appenzell, Toggenburg, Tauern Pied, Valais, and Drežnica goats, that exhibited the highest percentage of diversity.
The distinct genetic origin of the Drežnica goat was further demonstrated with the neighbor-joining tree (Fig. 6), where the breed formed its own branch. This statement was supported with analysis presented in Fig. 7, where Drežnica goat composed its own cluster that was clearly separated from other alpine breeds. Principal components calculated for all individuals of European breeds (Figure S3A) retained the cluster of Drežnica goat separately and close to other alpine breeds, but surprisingly grouping it together with Landrace goats from Netherlands. This was not so obvious in the neighbor net graph of Alpine (Fig. 6) and European breeds (Figure S1), where the closest node to Drežnica goat branch included three Austrian breeds (Pinzgau, Tauern Pied, and Styrian Pied goats) surrounded by Italian goat breeds. The cluster of Austrian breeds was also the closest on PCA graph of alpine breeds (Fig. 7), followed by a group of Italian breeds and the most distant Swiss breeds. For this reason, we added the third principal component (Figure S3B), which separated the Drežnica goat and Dutch Landrace goat. With the third principal component, Drežnica goat formed a subcluster in the middle of Italian and Austrian breeds. Admixture analysis (Fig. 8) revealed that a small group of Drežnica goat animals contains some admixture signatures of the Styrian Pied, a neighboring Austrian goat breed. Despite that, Drežnica goat is one of 11 breeds (8, 9, 16) with the most uniform population ancestral structure. In summary, Drežnica goat stayed genetically very homogeneous, which could be due to factors such as geographic and/or demographic isolation, bottlenecks, genetic drift, and distinctiveness or a combination of these factors.
Generally, the alpine breeds clustered according to the country of origin (Fig. 7) and geographical proximity, which was expected. One exception was a cluster of cosmopolitan breeds from France, Italy, Switzerland and Austria (FR_CMA, IT_CMA, CH_CHA, and AT_CHA). Moreover, admixture analysis (Fig. 8) revealed that a majority of the genetic ancestry is shared between the pairs IT_CMA-FR_CMA and CH_CHA-AT_CHA suggesting that they are most likely one genetically similar population. This was very likely due to sire semen exchange among these dairy breeds selected for high milk production. Another exception were Italian Passeier and Valpassiria goat positioned in the cluster of Austrian rather than Italian breeds (Fig. 7). The reason for this is probably the geographical location of Passeier and Valpassiria goat on both sides of the state border between North Italy and Southwest Austria. On the neighbor-joining tree (Fig. 6) Italian Passeier goat together with Valpassiria goat clustered together with Austrian Blobe goat. Additionally, the Passeier (this study) and Valpassiria (18) breeds could be regarded as one population based on our admixture analysis (Fig. 8), the later also demonstrating Passeier to Blobe introgression. As these breeds share a relatively small geographic area in the Tyrolean Alps, historical admixture was expected. Similar to the PCA results shown in Burren et al. (11), the Toggenburg and Appenzell goats had the same node, and next to that the Tessin Grey, Nera Verzasca and Peacock goats shared another node (Fig. 6). Tessin Grey, Nera Verzasca, Grisons Striped, and Adamello Blond goats displayed more heterogeneous population structures in our analysis (Fig. 8), compared to studies of Burren et al. (11) and Colli et al. (18).
The positions of populations in our study for the Alpine (Fig. 6), European (Figure S1), and Global datasets (Figure S2) were consistent with those on the PCA, neighbor net and phylogenetic graphs observed in Colli et al. (18). When we expanded the alpine neighbor net with other European breeds, the breeds from northern and northwestern Europe separated from the alpine breeds. Swiss breeds together with cosmopolitan alpine breeds were closer to French and Spanish breeds. On the other hand, the Drežnica goat and Austrian as well as Italian breeds were closer to Turkish breeds. A possible explanation of introduction of some Turkish goat stocks to Austria, Slovenia and further west to Italy could be in geographical closeness with the Turkish (Ottoman) empire. This was a state that encompassed in the 600-year period (14-20th Century) much of Southeastern Europe, Western Asia, and Northern Africa. In 1520, this empire expanded northwest all the way to what is today essentially the eastern border of Slovenia and Austria. The gene flow from the goat domestication center in Mid-East to Europe are thought to occur via two major routes, Danubian and Mediteranean (Amills_Tosser_Klopp 2017). This “Turkish” goat migration route via Balkan during the Ottoman empire could represent another most recent wave of south-Danubian-route introduction of goats to Europe. In the phylogenetic network of the global dataset, the breeds mainly clustered according to the continent. The European dataset ended up grouping with Spanish breeds and breeds from South America followed in the next node from the Spanish breeds. This is in line with historical facts regarding Spanish expansion in the early 16th century that also brought Spanish goats to this continent. We can conclude that most alpine goat breeds, especially local breeds with small population sizes, show relatively high homogeneous genetic structure and strong geographical partitioning, whereas larger population size cosmopolitan alpine breeds exhibit high admixture and geographic spread.
For the purposes of the second objective, we investigated the effect of inclusion or exclusion of outliers from the breed on the genetic diversity and population structure parameters. A composite test encompassing various metrics was used (24) for the detection of admixed outliers. The detection of purebred animals or excluding outliers is important, especially for analyses of small endangered populations of local breeds that are often in danger of crossing with highly selected commercial cosmopolitan breeds. Likewise, outliers in datasets could affect the estimation of genetic diversity and inference of population structure. Consequently, discrepancies in genetic parameters can be large when datasets with or without outliers are compared. For this reason, we formed different datasets using the multivariate outlier test to exclude outliers. If a repeated test still detected “new” outliers, they were again excluded in the second iteration. Finally, related animals were dropped out of the dataset in the last step. Our follow-up comparative analyses clearly demonstrated that datasets with or without outliers could affect the outcome of analyses. Major effects were observed in parameters such as the total number of observed alleles within a breed, number of private alleles, number of semiprivate alleles and mean number of alleles per block. To illustrate this issue for the case of the Drežnica goat, noticeable differences in results (Table 1, Table 3 and Fig. 3) were found when comparing the results of analyses of the Alp1Step and Alp2Step datasets for all the mentioned parameters. It is important to note that these effects did not arise due to differences in the number of animals per dataset, as each dataset included 50 goats, but rather due to differences in representative animals between the two datasets. Interestingly, we noticed also the effect on the diversity and structure parameters of breeds where no outliers or related animals were detected. For example, the Chamois Colored goat (Austria) and Peacock goat maintained the same number and the same animals across the datasets, but the results for these two breeds were also affected due to different dataset constructions of other breeds. In fast changing environment conditions, the level of within population genetic variation is one of the signals for extinction resistance of the breeds (4). This emphasized the need to select the optimal method to form the datasets for evolutionary and conservation genetics analyses. In case of employing the commonly used method in dataset formation without excluding outliers or highly related animals, there is a danger of getting the incorrect results and wrong assessment of actual adaptation capacity for the particular breed.
Changes in estimated genetic diversity parameters were observed in many breeds when different datasets were used. However, it was difficult to ascribe these differences to the compositional differences of the datasets versus the number of animals in the datasets. For this reason, we used the tools that account for differences in sample size among populations (Table 1 and Fig. 3). The relative contributions of the breeds to the total allelic diversity were different when we analyzed the Alp1Step or Alp2Step datasets. The Drežnica goat, Chamois Colored goat from Switzerland, and Camosciata Alpine goat from Italy displayed the largest changes between the datasets (Fig. 3). The differences in allelic diversity within populations were largely responsible for this change. We can explain this with the fact that excluding admixed outliers decreased the allelic diversity in particular breeds, which consequently also affected their contribution to the global allelic diversity. In the case of the Italian Camosciata Alpine goat, there was a loss of total allelic diversity after removing the breed from Alp1Step, but when the program removed the breed from Alp2Step, there was a gain in the total allelic diversity. The Camosciata Alpine goat from Italy showed the largest difference in mAR as well. In Alp2Step, this breed had around 11% lower mAR than in Alp1Step (Table 1). Together with the Swiss Chamois Colored goat, the Italian Camosciata Alpine goat is a member of the alpine breed cluster (Fig. 7) and these two breeds hence share a large proportion of alleles. The removal of admixed animals within individual breed causes a decrease in allelic richness and the average number of private alleles. As shown in Fig. 3A, AT_CHA and FR_CMA positively contributed to the total allelic richness in the Alp2Step dataset. In contrast, CH_CHA and IT_CMA were not among the breeds with a positive contribution. The reason for this is most likely that the input of CH_CHA and IT_CMA to total allelic diversity was mostly covered by AT_CHA and FR_CMA, since all four are cosmopolitan breeds that were generated from the same Alpine breed. This was not distinguished in the analysis with the Alp1Step dataset. Three breeds with the largest set of animals in AlpInit (IT_CMA, CH_CHA, and SI_DRZ) had the highest drop of mAR in Alp2Step. In their case, a larger number of animals consequently also means a larger number of outliers. This clearly demonstrated, how mAR in Alp1Step was overestimated, because of outliers that remained in the dataset. After detecting and removing them in Alp2Step, we got objective estimates of allelic richness for each breed. Considering both optimized datasets of alpine breeds, we can conclude that the exclusion of admixed and related animals reduces AS and AT in breeds with low allelic diversity, but increases AS and AT in breeds with high allelic diversity.
Analyses of the population diversity is often the basis for the criteria used in the management of breeding and conservation programs. Besides that, allelic diversity parameters are good indicators of long-term response to natural or artificial selection (5). Local breeds represent the majority of all datasets in this study and these breeds are particularly under long-term pressure of natural selection in harsh environments due to traditional extensive production systems. Estimates of allelic diversity can be directly linked to rates of adaptation and have a potential to be used as objective conservation criteria of each breed (44). As we demonstrated, it matters which animals we choose as breed representatives for the diversity analyses to obtain more objective results of genetic diversity and structure to form a basis for informed decisions in the breeding and conservation programs.
Differences in how the datasets were prepared also affected diversity parameters within populations, similar to relationship parameters between breeds. In fact, differences in the number of private alleles between datasets of a particular breed increased or reduced AT in Alp2Step, which contributed to better contrast in clustering and more distinct population differentiation. As an example is analysis in Tables S3 and Fig. 4, exhibiting that DEST values between the Drežnica goat and other breeds were lower in Alp1Step than in Alp2Step. This means that removal of admixed Drežnica goats resulted in higher pairwise distances with other populations than removal of closely related Drežnica goats only. Analyses with other breeds showed similar results. However, the largest differences in DEST values between Alp1Step and Alp2Step were found for the Italian Camosciata Alpine, Swiss Chamois Colored, and Drežnica goats most likely because these breeds had the largest sample size in the AlpInit dataset before optimizing Alp1Step and Alp2Step. Although the same or similar number of animals in Alp1Step and Alp2Step represented these breeds, the composition of animals in each breed was different between datasets. The highest difference in DEST values was calculated between IT_CMA and FR_CMA (deviation of 0.073, Table 1). As alluded earlier, semen exchange is common among cosmopolitan breeds, and genetically, these breeds resemble essentially one population (Fig. 8). Regardless, excluding the most admixed animals within each breed sharpened the differences between them, which was clearly demonstrated by the differences in distances based on Alp1Step or Alp2Step dataset. As expected, eliminating admixed animals from Alp2Step on average resulted in an increase in the pairwise distances between most breeds. When we removed significant outliers from each breed in the Alp2Step dataset, we obtained more objective distances between breeds in the alpine area.
Our results therefore clearly demonstrate that the procedure used for post-genotyping dataset optimization could have a significant impact on the outcome of genetic diversity and population structure analyses. Non-optimal dataset optimization could lead to erroneous conclusions about genetic diversity, identity and relatedness to other breeds. Choosing an objective method for exclusion of outliers can lead to more accurate and unbiased estimation of allelic diversity. Considering strong correlation of allelic diversity and long-term adaptation to the new optima (1, 2, 38) the improved estimation of allelic diversity could be considered as important part for improvement of the conservation prioritization. Taking all of our results together, we propose the mvOutlier test to be considered as a statistically proven, objective and effective tool for identifying outliers to allow for more reliable genetic parameter estimations, especially in local breeds with small population sizes. Likewise, this tool could be included in conservation and breeding programs to avoid or reduce breeding admixed animals in critically endangered populations. The vectors included in our mvOutlier analyses, originally proposed in Ramljak et al. (24), could be improved or replaced by more sophisticated vectors, and any optimization in this sense could be of broad interest.