Ethics statement

Blood sampling for GRF goats was conducted by veterinarians and no invasive procedures were applied. Thus, in accordance to the 2010/63/EU guide and the adoption of the Law D.L. 04/03/2014, n.26 by the Italian Government, an ethical approval was not required for our study.

Genomic data

Five millilitres of blood were collected from each of 48 female GRF goats (between 2 and 9 years old; all registered in the herdbook) from a pool of 269 females in the Garfagnina district (Media Valle del Serchio, Lucca, Central Italy). Animals were genotyped with the Illumina GoatSNP50 BeadChip (Illumina Inc., San Diego, CA) containing 53,347 SNPs [17]. Genomic data of nine Italian autochthonous goat breeds, namely Argentata dell’Etna (ARG), Bionda dell’Adamello (BIO), Ciociara Grigia (CCG), Di Teramo (DIT), Garganica (GAR), Girgentana (GGT), Orobica (ORO), Valdostana (VAL) and Valpassiria (VSS) were downloaded from the online repository (https://datadryad.org/stash/dataset/doi:10.5061/dryad.v8g21pt) of the ADAPTmap project [3, 18]. The breeds were selected based on breed abbreviation in the PLINK sample information file (.fam) downloaded from the repository and breed description (code and country) reported in Table 1 of [18]. The two datasets were merged and quality control was conducted in PLINK v1.9 [19, 20] based on the following criteria: i) only autosomes were kept, ii) call rate per SNP >95% and iii) missing values per sample <10%. After editing, 260 samples and 48,716 SNP were retained (Table 1). The distribution of SNP per chromosome (CHR) is presented in S1 Fig.

Table 1

Names of breeds, breed code and number of animals analyzed before (pre-QC) and after (post-QC) quality control per breed.

Breed name	Breed code	No. pre-QC	No. post-QC
Argentata dell’Etna	ARG	25	24
Bionda dell’Adamello	BIO	24	24
Ciociara Grigia	CCG	19	19
Di Teramo	DIT	24	24
Garganica	GAR	20	20
Girgentana	GGT	30	30
Garfagnina	GRF	48	48
Orobica	ORO	24	23
Valdostana	VAL	24	24
Valpassiria	VSS	24	24

Runs of homozygosity

Analysis of ROH was conducted in R (v. 3.5.0) using the package detectRUNS v. 0.9.5 [21, 22]. The consecutive method [23] that runs under the main function consecutiveRUNS.run was adopted. The required parameters were set to: i) minimum number of 15 SNPs/ROH, ii) 1 Mbp minimum length of ROH and iii) allow one heterozygous SNP within an ROH (to account for genotyping errors). In addition, ROH lengths were split into five classes (0–2, 2–4, 4–8, 8–16 and >16 Mbp). For each of the class and breed, descriptive statistics of ROH per breed, per chromosome, per SNP and per length class were estimated. Principal component analysis (PCA) was used to identify (dis)similarities among breeds, relative to the average number of ROH identified per chromosome. PCA was applied via a singular value decomposition (prcomp function in R [21]). In addition, genomic inbreeding (F_ROH) was calculated per breed (defined as the length of ROH over the total autosomal length per goat and summed over all goats per breed). Regions with a high frequency of ROH (≥45%) were detected and genes located within ± 1 Mbp were annotated by using the Capra hircus ARS1 (http://www.ensembl.org/index.html) and the variant effect predictor (https://www.ensembl.org/Tools/VEP) Ensembl databases.

Population stratification and ancestry

PCA and admixture analysis were used to infer the presence of distinct populations based on the SNP data. PCA was conducted on the matrix of genotypes. The proportion of mixed ancestry in the breeds was assessed by the ADMIXTURE 1.22 software [24, 25]. The number of ancestries (K) to be retained in the admixture analysis (K = 2 to 10) was evaluated via a 10-fold cross-validation (CV). The final selection on the number of ancestries was done by inspecting the CV error.

Discriminant analysis of principal components

Discriminant analysis was applied to assess the breed traceability of the GRF goats using SNP data. To achieve this, the methodology of discriminant analysis of principal components (DAPC) [26] implemented in the adegenet R package [21, 27, 28] was adopted. In brief, DAPC is a 2-step approach: firstly, a PCA on the matrix of the genotypes is conducted and then, a small number of selected PCs (instead of the original SNP genotypes) is used as an input for the linear discriminant analysis (LDA). The selection on the optimal number of PCs to be further used in the LDA is done via cross-validation (CV) where the data is split in training and validation sets. The following criteria were implemented for selection of PCs: i) 10-fold CV with 30 repetitions, ii) a maximum number of 259 PCs were tested, and iii) the number of PCs to be retained was based on number of PCs associated with the highest mean success. Three different scenarios of DAPC were applied as described below:

Scenario 1 (supervised learning). The full dataset was analyzed simultaneously. In this scenario, all available data were used for model training and the discriminant functions were extracted based on all animals. This is not, however, a real case scenario, since the discriminant functions were developed utilizing the entire data set. The objective for a practical application is to identify an external individual membership to a group (i.e., external validation). Hence, two more scenarios were developed adopting a CV scheme also for the discriminant function.

Scenario 2 (semi-supervised learning). Assessment of correct assignment of GRF goats was done via a semi-supervised CV (CV_SS). Five GRF goats were sampled representing the testing set of the DAPC analysis. The training set (TRN) was constituted by the remaining of 43 GRF goats plus all the goats from the other breeds. The five GRF samples were classified in one of the 10 breeds presented in the TRN set via the function predict.dapc. The procedure was repeated 10 times and results were averaged over the 10 repetitions.

Scenario 3 (unsupervised learning). Assignment of GRF goats in a breed but without the presence of any GRF goats in the TRN set (unsupervised CV; CV_US). This scenario could also be viewed as a method to assess the genomic similarity of the GRF with the rest of the breeds (i.e., type of clustering). The approach was similar to Scenario 2 other than the testing population consisted of the entire GRF set and GRF samples had to be classified in one or more of the other 9 breeds. To increase the number of the tested samples in each round of the CV, 80% of GRF breed was sampled. Moreover, to test for the effect of the size in model training for the assignment of GRF, different proportions of TRN set were sampled (20, 30, …, 90%) 10 times each, and results were averaged over 10 repetitions. In other words, the size of the TRN set varied between 42 to 191 goats. All nine breeds were present in each scenario and all GRF goats were used in this scenario.

The terms (semi/un)-supervised should not be confused with the terminology used in machine learning. The introduction of these terms has been used in the manuscript to distinguish among the three approaches that were used in the DAPC analysis, and, although they are, up to a point, analogous with the same terms used in machine learning they are not identical.

Runs of homozygosity

Summary results of the detected ROH regions as either total counts or averaged based on the number of samples per breed are presented in Table 2 and Fig 1a, respectively. A relatively high number of ROH was detected for GRF (n = 2,450), which was the third largest among the breeds in the study. The greatest number of ROH among all breeds was for GGT (n = 2,762) followed by ORO (= 2,693), while the smallest was found for ARG (n = 465). For GRF, the number of ROH per Capra hircus chromosome (CHI) varied from 35 (CHI25) to 158 (CHI1). The maximum length of ROH per chromosome was found on CHI1 (568,887,711bp) and the minimum on CHI23 (113,323,345bp). In general, the total length of ROH per CHR followed the same pattern of the total ROH number per CHR (Fig 1b).

Fig 1

Summary results of runs of homozygosity (ROH): (a) average number detected per breed, (b) number and length of ROH per chromosome (CHR) in Garfagnina breed, (c) frequency distribution of the number of ROH in the breeds analyzed per length class and (d) frequency distribution of the number of ROH in different length classes per breed.

ARG: Argentata dell’Etna; BIO: Bionda dell’Adamello; CCG: Ciociara Grigia; DIT: Di Teramo; GAR: Garganica; GGT: Girgentana; GRF: Garfagnina; ORO: Orobica; VAL: Valdostana and VSS: Valpassiria; In (a) horizontal bars within each boxplot represent the median, and red rhombus the mean.

Table 2

Total number of runs of homozygosity (ROH) detected per breed.

Breed	No. ROH
ARG	465
BIO	814
CCG	496
DIT	813
GAR	730
GGT	2,762
GRF	2,450
ORO	2,693
VAL	1,613
VSS	768

For all breeds analyzed except CCG and DIT, the number of ROH relative to the length on the genome decreased with an increase in length (Fig 1c and 1d). In DIT samples, ROH were more frequent in length classes of 4–8 and 8–16 Mbp compared to 2-4Mbp. The percentage of ROH with a length >16Mbp per breed varied between 0.89% to 13.53% for ORO and DIT, respectively. For small ROH length (<2Mbp) the proportion over the total number detected reached ~78% in ARG, while only 35% of ROH was observed for DIT. The pattern of ROH length class was similar among GRF, GGT and ORO with ~50% of ROH having a length < 2Mbp, ~25% between 2-4Mbp, ~13% between 4-8Mbp, ~5% and ~2% between 8-16Mbp and > 16Mbp, respectively (S1 Table).

For GRF, common ROH (more than 45% in the GRF samples analyzed) were detected on CHI12, between ~34.6–35.3 Mbp (Fig 2, Table 3). In total, 14 SNP were contained in this region. The same area was also present in DIT, while a broader region (~33,9–36.5 Mbp) was shared among ARG, CCG, and GGT. To identify similarities among breeds relative to the number of ROH per chromosome, a PCA was conducted on the average number of ROH identified per chromosome. In addition, a heatmap on the actual number of ROH per chromosome was produced. Both approaches placed GRF closer to ORO and GGT with respect to the rest of the breeds (S2 Fig).

Fig 2

Number of times (%) each SNP was detected inside a run of homozygosity (ROH) in Garfagnina (GRF) goats.

Table 3

Most common (≥45% in each breed) runs of homozygosity (ROH) detected per breed on Capra hircus chromosome 12, with the start-end regions and number of SNP per ROH.

Breed	Start-SNP	End-SNP	No. SNP	Start-region (bp)	End-region (bp)
ARG	snp30406-scaffold335-807385/ rs268262652	snp30399-scaffold335-501928/ rs268262645	6	34,949,980	35,255,437
CCG	snp30421-scaffold335-1558038/ rs268262666	snp30391-scaffold335-171300/ rs268262637	27	34,199,327	35,586,065
DIT	snp30413-scaffold335-1113038/ rs268262659	snp30397-scaffold335-418126/ rs268262643	14	34,644,327	35,339,239
GGT	snp30428-scaffold335-1839052/ rs268262673	snp11142-scaffold140-760668/ rs268243983	53	33,918,313	36,518,132
GGT	snp17454-scaffold1805-39262/ rs268250096	snp55342-scaffold855-361968/ rs268286936	26	40,620,364	42,257,561
GRF	snp30413-scaffold335-1113038/ rs268262659	snp30397-scaffold335-418126/ rs268262643	14	34,644,327	35,339,239
VAL	snp50169-scaffold717-4207960/ rs268281880	snp3193-scaffold1095-4352995/ rs268236233	87	24,592,901	28,744,348

Genomic inbreeding coefficients (F_ROH) were found to be intermediate for GRF compared to the rest of the breeds analyzed, with a mean value of 0.069 (Fig 3a, S2 Table). The highest values were observed for GGT (0.143) and ORO (0.137). Moreover, the distribution of F_ROH calculated per CHR was similar, with some high values (>0.5) observed for CHI7, 9, 16, 22 and 25 (Fig 3b).

Fig 3

Summary of the genomic inbreeding coefficients (a) per breed and (b) of the Garfagnina breed per chromosome (Chr).

ARG: Argentata dell’Etna; BIO: Bionda dell’Adamello; CCG: Ciociara Grigia; DIT: Di Teramo; GAR: Garganica; GGT: Girgentana; GRF: Garfagnina; ORO: Orobica; VAL: Valdostana and VSS: Valpassiria.

Population stratification and ancestry

As a first step, a PCA was conducted on the complete data to visualize the general structure and relationships among breeds. The first axis distinguished GRF goats from ARG, CCG, DIT, GAR and GGT, while the second axis further separated GRF from the rest of the breeds (Fig 4). An inspection of all the pairwise comparison between the first 10 axes (PCs) was carried out. In general, GRF was clearly separated from the rest of the breeds. Interestingly, a small GRF subgroup consisting of six goats was observed on the 9th axis (S3 Fig).

Fig 4

Scatterplot of the first two principal components¹.

¹Singular value decomposition was applied on the matrix of genotypes. ARG: Argentata dell’Etna; BIO: Bionda dell’Adamello; CCG: Ciociara Grigia; DIT: Di Teramo; GAR: Garganica; GGT: Girgentana; GRF: Garfagnina; ORO: Orobica; VAL: Valdostana and VSS: Valpassiria.

An admixture analysis was conducted to complement the PCA results. A varying number of group ancestries was investigated, ranging from K = 2 up to 10. The model with the minimum CV error was the one with eight group ancestries (S4 Fig). In general, the admixture results were in agreement with PCA, depicting the uniqueness of the GRF gene pool (Fig 5a and 5b). At K = 8 (number of ancestries selected after CV) GRF had almost a breed-specific ancestry, with a small percentage of the GRF goats sharing co ancestry with i) BIO and VSS and ii) ARG and CCG and to a small extent with ORO, VAL, DIT and GGT. It should be noted that apart from GRF, group-specific ancestries, at least to a great extent, existed almost for all breeds except ARG, CCG, GAR and VSS.

Fig 5

Admixture analysis with K = 8 coancestry groups (a) per individual and (b) averaged per breed.

ARG: Argentata dell’Etna; BIO: Bionda dell’Adamello; CCG: Ciociara Grigia; DIT: Di Teramo; GAR: Garganica; GGT: Girgentana; GRF: Garfagnina; ORO: Orobica; VAL: Valdostana and VSS: Valpassiria.

Discriminant analysis of principal components

In the first scenario of DAPC, all data were used. The first 40 PCs, explaining ~35.75% of the total variability in the SNP data (S5 Fig), were used in the final DAPC model, resulting in an assignment success rate of 100% for all GRF goats to their breed of origin. The pattern of the genetic diversity based on DAPC is presented in Fig 6, where a clear genetic distance of GRF from the rest of the breeds can be observed.

Fig 6

Scatterplot of the first two discriminant components of the DAPC.

Breeds are presented by different colors and symbols. ARG: Argentata dell’Etna; BIO: Bionda dell’Adamello; CCG: Ciociara Grigia; DIT: Di Teramo; GAR: Garganica; GGT: Girgentana; GRF: Garfagnina; ORO: Orobica; VAL: Valdostana and VSS: Valpassiria.

An external validation scenario, which better reflects a practical application of the discriminant model, was further assessed. In the first analysis (CV_SS scenario), GRF had representative animals in the TRN set where the DAPC model was developed. Also, in this case, a 100% correct classification of GRF goats was observed (S3 Table). Interestingly, the classification of GRF was invariant to the number of PCs selected (ranged between 10 to 70) in DAPC (S4 Table). In the second scenario (CV_US) there were no representative GRF samples in training the model of DAPC. In that case the majority of animals were classified as CCG while few were assigned to DIT in some of the CV replicates (S5 Table). Similar results were obtained with an increased size of TRN set. In the majority of the scenarios, GRF goats were classified either as CCG or DIT, where there were a few cases in which GRF goats were also assigned as VSS, GAR or BIO, but in none of the cases as ORO, GGT or VAL (Fig 7).

Fig 7

Percentage of assignment of the GRF goats in the CVUS scenario.

Results were averaged over 10 replicates (CV) in each data subset (from 20 to 90%). CVUS: unsupervised CV, where GRF breed had no representative goats in model training, hence GRF goats had to be classified in one of the other 8 breeds; ARG: Argentata dell’Etna; BIO: Bionda dell’Adamello; CCG: Ciociara Grigia; DIT: Di Teramo; GAR: Garganica; GGT: Girgentana; GRF: Garfagnina; ORO: Orobica; VAL: Valdostana and VSS: Valpassiria.

Runs of Ηomozygosity

As it has been highlighted by Bertolini et al. [35], crossbred goat populations tend to have smaller total ROH length and number compared to purebred populations. The same pattern has been observed when comparing unselected vs. selected populations undergoing breeding programs. Nevertheless, as it has been pointed out by different colleagues [35, 36] the 50k chip is not considered efficient for the accurate detection of small ROH regions, resulting in underestimation of small ROH hits. Despite this, our analysis was based on purebred goats, and thereby a smaller bias is expected. Results of ROH and F_ROH for the 9 Italian breeds of the AdaptMap project were in agreement with results previously reported [35].

The general pattern of ROH (i.e. in terms of total—and by chromosome—number and length of ROH) for GRF was similar to GGT and ORO. Moreover, common ROH were found for GRF (more than 45% in GRF samples analyzed) on CHR 12 at, roughly, between 34.6–35.3 Mbp (Table 3). The same region was also detected in DIT breed, while the broader region (~33,9–36.5 Mbp) was shared among ARG, CCG, and GGT. Further, a search of genes presented in the top ROH region identified for GRF (~50.25–50.94Mbp, updated on the ARS1 assembly) and 1Mbp up-downstream (~49.25–51.94Mbp, updated on the ARS1 assembly) was carried out. Interestingly, the region ~49-52Mbp has been previously reported in goats [35, 37]. It is worth noting that within this region lay the genes of the general gap junction protein family GJA3 (gap junction protein alpha 3; ~50.642–50.644Mbp), GJB2 (gap junction protein beta 2; ~50.675–50.676Mbp), GJB6 (gap junction protein beta 6; ~50.694–50.695Mbp). The GJB2 and GJB6 are associated with the nervous system, hearing functions and ectodermal processes [38, 39]. Moreover, the SAP18 (Sin3A associated protein 18; ~51.136–51.141Mbp) that is related to gonad development [40], was also mapped in this region.

The narrow region of the detected top ROH runs for GRF on CHI12 spanned the CENPJ (centromere protein J; 50.23–50.27Mb) and the IL17D (interleukin 17D; 50.91–50.93 Mb). More precisely, the snp30397-scaffold335-418126 was found to be an intron of the CENPJ gene, while the snp30413-scaffold335-1113038 was downstream from IL17D. There is a series of studies that have linked the CENPJ with primary microcephaly in humans and mice [41–44]. Moreover, CENPJ has been found to regulate the neurogenesis and the cilia disassembly in the developing cortex in mice [45]. Also in mouse, disruption of the CENPJ can cause the Seckel Syndrome [46]. Possible associations with these physiological functions and phenotypic differences among the breeds studied are not obvious.

Apart from the CENPJ and IL17D genes, 9 more genes are found within this genomic area, namely PARP4 (poly(ADP-ribose) polymerase family member 4; ~50.29–50.35Mbp), MPHOSPH8 (M-phase phosphoprotein 8; ~50.36–50.41Mbp), PSPC1 (paraspeckle component 1; ~50.44–50.48Mbp), ZMYM2 (zinc finger MYM-type containing 2; ~50.56–50.63Mbp), as well as the CRYL1 (crystallin lambda 1; ~50.77–50.83Mbp) and the IFT88 (intraflagellar transport 88; ~ 50.84–50.91Mbp).

Downstream this region, in 1Mbp expansion, the genes ATP12A (ATPase H+/K+ transporting non-gastric alpha2 subunit; ~50.08–50.11Mbp) and RNF17 (ring finger protein 17; ~50.11–50.23Mbp) were located. Moreover, upstream the region there were also mapped the EEF1AKMT1 (EEF1A lysine methyltransferase 1; ~50.94–50.95Mbp), LATS2 (large tumor suppressor kinase 2; ~51.07–51.09Mbp), ZDHHC20 (zinc finger DHHC-type containing 20; ~51.16–51.22Mbp), MICU2 (mitochondrial calcium uptake 2; ~51.23–51.28Mbp), and the FGF9 (fibroblast growth factor 9; ~51.34–51.37Mbp).

Breed diversity, ancestry and discrimination

Two approaches, complementary to each other, have been used to infer the GRF relationships with 9 native Italian breeds, namely principal component and admixture analysis. Both analyses confirmed the distinct and unique genetic background of the GRF breed and results were in general agreement with each other. Overall, GRF was placed closer to VSS, BIO and CCG, with a small percentage of ancestries shared with all 9 breeds. Interestingly, both approaches identified a subgroup of 6 goats. In admixture, the six goats showed a unique ancestry (Fig 5a), while this subgroup was identified in PCA on the 9th axis (S3 Fig). Nevertheless, the percentage of variance explained by this axis was <1%.

The DAPC model was used to assess breed traceability based on SNP and was able to classify with 100% success GRF goats to their breed of origin (S3 Table). Moreover, an unsupervised learning was applied, where GRF had no representative samples in TRN set. Results were consistent with the PCA and admixture and assigned the majority of GRF goats to the CCG breed with a small number of goats (varied between 4 to 6) assigned as DIT (S5 Table).

As mentioned in the material and methods section, the primary step of the DAPC analysis is to select the number of PCs to be used in the discriminant model. Hence, a basic question emerged on how robust the DAPC could be considered relative to the number of PCs used. Our analysis showed that, although the assignment success was invariant to the number of PCs in the semi-supervised DAPC analysis (number of PCs varied between 10–70), a pattern was found in the case of the unsupervised model. More precisely, when the DAPC contained 40 PCs some of the goats were classified as DIT. In the rest of the cases, where 60 or 20 PCs were used, all GRF goats were assigned as CCG.

PCA analysis seems to detect common ancestries between GRF goats and Alpine Arc goat breeds whereas the DAPC approach identifies similarities between GRF and the goat breeds from Central Italy. Both hypotheses are consistent with the history of Tuscan goat populations that experienced migratory flows from both northern and central Italy.

Overall, the genomic analysis confirmed the hypothesis that GRF is a result of crosses among goats from the Alpine Arc and Tuscan-Emilian Apennines regions. Nevertheless, based on the genomic information analyzed here, GRF represents a unique genetic pool and was genetically distinguished from 9 different native Italian breeds. This, in turn, resulted in breed traceability with a 100% success rate after CV. To sum up, our analysis complements previous work on various zootechnical and adaptive characteristics [7–8] of GRF and provides with a more complete description of the breed.

Suggestions for monitoring and conservation

One of the strategic priorities of the Global Plan of Action for Animal Genetic Resources of FAO is breed conservation, either in vivo or in vitro or both [29, 32]. The population status of GRF could be described and summarized in the following points: GRF i) has a distinct genetic pool, ii) exist in a small number of farms and animals, iii) is concentrated in a small geographical region, iv) lacks a formal breeding strategy, v) is reared in semi-extensive and family type farming systems, vi) is not subject to reproductive technologies such as artificial insemination, vii) has low pedigree quality (or almost absent) and viii) lacks of an organized action plan for monitoring and conservation. These points contribute to the critical status of the breed and highlight the need of conservation policies. For conservation in farms, farmers need to be financially supported by local and international authorities. A direct link between the breed and the production of niche products could add an extra economic value. Moreover, the cultural and the environmental components should be further investigated and quantified in the near future.

In addition, storage of biological material in gene banks offers an extra level of security. As recently been reported by Zomerdijk et al [29], gene banks greatly vary among countries and species and goats come second in order in gene banks after cattle. Even though gene banks, are cost-efficient [47, 48], they face various challenges, such as economic losses, loss of genetic material due to lack of liquid nitrogen and risks from diseases and catastrophic events such as floods and earthquakes [29]. Hence, the storage of biological material in a gene bank should be seen as a complementary rather a substitution to in vivo conservation.

Moreover, the idea of the foundation of a GRF breeding association should be considered. Animals are used by farmers if they provide them with profit and profitability is strongly associated to productivity. In this regard, research institutes should provide GRF farmers and breeders with vital scientific support. For example, the development of new tools for breed traceability and recognition might increase farmers’ interest. A breeding nucleus, where phenotypic and genetic variability are studied, should be considered. The startup of a breeding scheme, however, should be carefully reflected, since in the Italian legislation, no financial support is provided for livestock conservation if the breed undergoes directional selection, even if the population is at risk of extinction.

Furthermore, genotyping of a large number of animals (males and females) is highly recommended. To reduce costs, genotyping could be targeted, at a first step, only to animals that will contribute to the next generation. In the near future, genomic information could be further utilized for the development of a genomic-based breeding program that would help to select animals early in life and boost the genetic progress. Besides, optimal contribution selection [49, 50], penalization of the number of offspring per male (in the absence of a breeding objective) and a rotational mating scheme, with a controlled buck exchange among farms [51], are some ideas to be considered.

In the three pillars for a successful conservation scheme described in [34] (“discover, secure and sustain”), GRF is still under investigation. Various studies have been carried out, while our work presented a first whole-genome scan on the breed’s genetic diversity.