Assessing major genes allele frequencies and the genetic diversity of the native Aosta cattle female population

Abstract

The Aosta cattle breeds have a key role in the economy of the Aosta Valley. In addition to the meat and milk production these autochthonous breeds are important for their cultural value, and for their role in the maintenance of the mountain environment. The knowledge of their genetic makeup represents a fundamental asset to managing the reproduction of the population in order to maintain the existing genetic diversity and, as a possible input, to apply genomic selection in a small population. A total of 3195 Aosta cows were genotyped with the GeneSeek Genomic Profiler® (GGP) Bovine 100K by Neogen in the framework of the DUALBREEDING-2 project. The Aosta Black-Chestnut and Chestnut-Herèn resulted to be a unique population. The Aosta female population had shorter ROH than bulls and shared ROH_islands that harbour adaptative and functional genes. Allele frequencies of major genes highlight the possibility for selection for both milk and meat quality variants and that the Aosta cattle population is free from the known Mendelian inheritance diseases found in cosmopolitan breeds (e.g. BLAD, CVM, HCD).

HIGHLIGHTS

• The autochthonous Aosta cattle breeds are free from the Mendelian inheritance disease variants found in cosmopolitan cattle.

• Allele frequencies revealed the effort in selection for variants related to cheese yield and a viable selection for variants related to meat quality.

• The use of genomic information in Aosta cattle can help maintain the biodiversity of these local breeds and their peculiarities.

Keywords:

• Aosta cattle breeds
• autochthonous cattle
• genetic diversity
• ROH
• Mendelian traits

Introduction

The autochthonous Aosta dual-purpose cattle is raised in the northern-west alpine region of Italy, Valle d’Aosta. The Aosta Red Pied breed (ARP) most likely derives from the Northern Europe pied breeds, brought into the valley around the fifth century by the Burgundians migrations, while the Aosta Black Pied-Chestnut (ABC) with the Swiss Hérens are most likely native cattle breeds sharing a common genetic background (Vezzani 1929; Del Bo et al. 2001; Strillacci et al. 2020; Signer‐Hasler et al. 2023).

According to the herd book regulations approved in 2020 (https://anaborava.it/lg_disciplinare.html), the Aosta cattle population genetic improvement program has two sections: one for the ARP and one for the ABC cattle and the ABC subjects with a Hérens ascendant, coded as Aosta Chestnut - Hérens (ACH).

Since the establishment of the National Breeders Association (Associazione Nazionale Allevatori Bovini di Razza Valdostana - ANABORAVA) in 1980 the selection criteria of Aosta cattle includes milk and meat production. In addition to production traits, farmers always paid attention to the ability of cows to adapt to the farming system based on summer grazing. Summer pasture occurs in the mountains surrounding the Aosta valley up to 2500 m.a.s.l. that is characterised by rough climates and permanent meadows. Even if the grazing activity in such harsh environmental conditions may challenge the cows’ functionality, summer pasture is practised as it represents a great economic value for the valley preserving territories and landscapes, making them accessible to tourism and thus closely linked to the local economy and traditions.

For what concerns productivity, in 2022 the average 305-day milk yield recorded by the Italian breeders Association (AIA) were: 3334 kg and 2459 for ABC and ACH respectively and 3781 kg for ARP (http://bollettino.aia.it/Contenuti.aspx?CD_GruppoStampe=TB&CD_Specie=C4). The milk from the Aosta cattle is mainly used to produce the Fontina cheese which is a Protected Designations of Origin (PDO) and its regulation envisages that it must be produced in the Aosta Valley, exclusively with Aosta cattle milk.

Regarding the cultural and social value of the breeds, a historical tournament, called “Battailes des Reines”, takes place yearly in the Aosta Valley. During this folkloristic event cows challenge each other in a harmless match to establish a hierarchy, thus determining the dominant one. This traditional contest takes place in the span of time from March to October, with the final event at the “Arena Croix Noir” to proclaim the Queen of the Valley (Sartori and Mantovani 2010). Because of this tournament, the combativeness trait has been included in the official selection index of ABC and ACH since 2014, the Cheese yield, Muscularity and Combativeness Index - IRCMC (Sartori et al. 2014). Whereas the selection index for the ARP is the IRCM, based only on Cheese yield and Muscularity.

Since the introduction of SNP genotyping technology, the genomic information has been used to investigate the Aosta breeds’ genetic variation (Strillacci et al. 2020; Signer‐Hasler et al. 2023) as a step forward from the studies based on SSR markers (Del Bo et al. 2001).

Among the indicators of genetic variability, Runs of Homozygosity (ROH) are long portions of DNA in homozygous state that may be shared across individuals and populations (Peripolli et al. 2017), i.e. genomic regions identical by descent (IBD). ROH are investigated in populations as they may occur as a consequence of directional selection (Falchi et al. 2023) in livestock. They are generally used to investigate how the population's history, structure and demography evolved along generations (Peripolli et al. 2017).

The ROH can also be used to calculate the genomic inbreeding coefficient (F_ROH) which is the proportion of ROH-covered genome to the total length of the autosomal genome (Saura et al. 2015). Another commonly used inbreeding coefficient, based on the difference between the observed and expected homozygous genotypes is the F_HOM that can be calculated for each individual and across all markers (Keller et al. 2011).

For the conservation of animal genetic resources, in addition to a precise and accurate phenotype description, it is essential to frame the genetic variability that exists in autochthonous populations. The objectives of this study were: (i) to identify the genetic diversity among the Aosta breeds, (ii) to identify Runs of Homozygosity in Aosta cattle population; (iii) to estimate the genomic inbreeding of individuals; (iv) to determine allele frequencies at loci that influence milk and meat characteristics, mendelian diseases, and fertility haplotypes.

Materials and methods

Ethics statement

For this study the Animal Care and Welfare Committee approval was not required as genotypes were made available by ANABORAVA from their genomic database.

Sampling and genotyping

The 3195 female genotypes, obtained with the GGP Bovine 100K (GeneSeek®) by Neogen, were provided by the Associazione Nazionale Allevatori Bovini di Razza Valdostana (ANABORAVA) 929 - ABC, 158 - ACH and 2108 - ARP. The proportion of breeds in the genotypes well represents the population composition: in fact, in 2020, the Aosta breeders’ association registered a total of 5,598 ABC-ACH and 13,997 ARP cows in the herd book (https://anaborava.it/lg_consistenza.html). Cows were sampled across all farms with individuals recorded in the herd book and represent as such the actual population.

The genotypes part of the ANABORAVA genomic database have been produced by the breeder association within the national project DUALBREEDING - Fase2, funded by EU EAFRD and the Italian Ministry of Agriculture. The aims of the project are to preserve cattle populations’ biodiversity and to integrate the information of the DUALBREEDING project for a better genotypic and phenotypic characterisation.

The available marker dataset consisted of 89,762 SNP markers mapped on autosomes with known chromosomal positions according to the ARS-UCD1.2 bovine reference genome. This SNP dataset does not include redundant SNPs to detect the same mutation. All individuals had a call rate value >95%.

Population structure and genetic diversity analysis

Genotypes were subjected to quality control: SNP with a call rate lower than 95% were filtered out. To analyse the population structure and genetic variation within and among breeds the following approaches have been used:

• The principal component analysis (PCA) has been performed in SNP & Variation Suite (SVS) v8.9 (Golden Helix Inc., Bozeman, MT, USA) using only the SNP genotypes filtered for MAF >0.05; leaving a total of 84,078 SNPs for the analysis. Two analyses have been performed one for the three breeds together and one for only the ABC and ACH breeds. The graphical representation of the results has been obtained using the R package “ggplot2” (Wickham 2011).

• The pairwise fixation index (i.e., Wright’s F-statistic - F_ST) was estimated using the dedicated module implemented in SVS v8.9. The F_ST was calculated for all possible pairs of breed combinations.

• The observed and expected heterozygosity has been obtained by PLINK v1.9 software using the “–hardy” command (Purcell et al. 2007).

• The admixture analysis has been performed with the ADMIXTURE software (Alexander et al. 2009). To correct for individual relatedness the pairwise LD pruning (using the PLINK v1.9 software) has been performed on a sliding window of 50 SNPs with a step of 5 SNPs and an r² of 0.5. After the pruning a total of 58,836 SNPs were available for the admixture analysis. The Admixture analysis has been carried out using from 2 to 7 K and the one with the lower cross-validation error (CV) has been selected (K = 2, CV = 0.58851).

Effective population size

The estimation of effective population size (Ne) was obtained using the SNePv1.1 program as described in Barbato et al. (2015). The approach is based on the known relationship between the variance in linkage disequilibrium (calculated using allele frequencies) and the effective population size. The software can estimate the historical Ne based on the relationship of the LD spread, directly calculated on the SNP genotypes, and on the recombination rate. The recombination rate between a pair of SNPs is inferred by SNePv1.1 considering the relationship between the physical and genetic distances, with a default value of 1 Mb = 1 cm.

Detection of runs of homozygosity

The ROH detection was carried out with the “detectRUNS” package of the R software (Biscarini et al. 2018), using the “consecutiveRUNS” method with the following parameters: i) the minimum length of the ROH’s was set to 1 Mb to avoid the detection of short and common ROH across the genome due to linkage disequilibrium (LD), because no LD-based pruning was performed; ii) a minimum of 30 homozygous consecutive SNPs; iii) no heterozygote nor missing genotypes were allowed; iv) a maximum gap of 1 Mb between the SNPs, to ensure that the SNP density did not affect the ROH.

The top 1% SNPs of the SNP distribution in ROH were used to identify the ROH_islands. This defines a threshold of 26% and 51% of SNPs in ROH, i.e. the proportion of individuals with uninterrupted sequences of the top 1% SNP in a ROH. The QTL annotation of the genes inside the ROH_islands, identified with the gene set (Bos taurus: Annotation Release 105) downloaded from NCBI, has been obtained using the Animal QTL database for Cattle (https://www.animalgenome.org/cgi-bin/QTLdb/BT/index).

Inbreeding coefficients calculation

The inbreeding coefficients (F_HOM and F_ROH) were calculated for each sample as follows:

• the F_HOM coefficient values were calculated by the dedicated function of the SVS v8.9 software.

• the F_ROH inbreeding coefficient has been calculated after ROH detection as:

$${\mathrm{F}}_{\mathrm{ROH}}\mathrm{ }=\mathrm{ }\frac{{{\sum }_i^n\mathrm{L}}_{\mathrm{ROH}}}{{\mathrm{L}}_{\mathrm{AUT}}}$$ where L_ROH is the length of the ith ROH segment, n is the number of detected ROH and L_AUT is the length of the autosomal genome covered by the SNPs (2,487,916,500 bp).

For the inbreeding values obtained, the minimum, maximum, mean, and standard deviation values for each breed were then obtained, to compare the results across them.

Genotypic frequencies for the traits of interest

Several SNP mapping causative loci are included in the Neogen GGP Bovine 100K SNP chip. Among them, the ones related to known cattle mendelian traits and disorders have been investigated (e. g. diseases, fertility haplotypes and milk proteins). Genotypic frequencies were obtained within breed for each locus using SVS v8.9.

Results and discussions

Population structure of the Aosta breeds

Figure 1A shows the dispersion plot of the first 2 components of the PCA representing 59% of the total variability. The ARP grouped in a separate cluster, suggesting being genetically differentiated from the ABC and ACH breeds that cluster in a unique group, supporting the hypothesis of their common origin. In Figure 1B it is possible to see the graphical representation of the PCA analysis performed including only genotypes of the subjects of the two breeds ABC and ACH, that cluster overlapping one another as in Figure 1A. This result is in line with Strillacci et al. (2020) where the Aosta Black Pied (ABP) and Chestnut (CAS) were classified as 2 different breeds. To compare their results with those of this study, it is to be accounted that the ABP in their study, is here part of the ABC, while the CAS here is classified as ABC and ACH according to the ascendants, as hereinbefore described. Additionally, the comparison with the results by Signer‐Hasler et al. (2023) confirmed the subdivision in two main clusters: one for the Aosta Black Pied, Chestnut, Herén (VPN, CAST and ER respectively in their study) and one for the Aosta Red Pied (VPR in their study). The sample size in this study is bigger and more evenly distributed across the Aosta cattle, reflecting the actual population-to-breed proportion. The cluster of the ARP is showing a small separated cluster. This is particularly interesting as the sub-cluster (lower left corner of Figure 1A) is composed of daughters of only one Aosta Red Pied bull. All sires in reproduction are registered in the Herd Book and must have, as this one, registered sire and dam and must be genomically tested for their paternity and maternity. Ancestors (i.e. paternal and maternal grandsires) are in common with other sires of the analysed females. The hypothesis we may speculate is that the difference could be due to a specific variation of its genome with respect to the rest of the population.

Figure 1. A) Graphical representation of Principal Component Analysis (PCA) result for the three breeds (ARP, ABC and ACH): Principal component 1 (PC_1) vs Principal Component 2 (PC_2) are plotted; B) Graphical representation of Principal Component Analysis (PCA) result limited to ABC and ACH: Principal component 1 (PC_1) vs Principal Component 2 (PC_2) are plotted; C) Aosta breeds Wright’s fixation index - F_ST.

The PCA results are supported by the F_ST index (Figure 1C): the difference between the ABC and ACH is 0.005, practically close to zero, meaning that the two breeds can be considered as a unique population. In fact, the value is much less than the F_ST obtained between ARP and the other two ones, 0.053 and 0.06 for the ABC and ACH respectively. Strillacci et al. (2020) identified similar F_ST values between the ARP and the other two breeds, i.e. 0.050 and 0.052 for the ABP and the CAS respectively.

The Admixture analysis identified two ancestors (K = 2). The results are shown in Figure S1, reporting the cross-validation error for K 1 to 4. At K = 2, which is the one with the lowest cross-validation error, the ABC and ACH were composed of 90% and 97% by the same population respectively, and the ARP breed consited of a different ancestor (94%).

Based on this evidence (low F_ST values and high genetic similarity), all the subsequent analyses have been developed considering the breeds grouped in ARP and ABCH, which is the union of ABC and ACH.

The Ne calculated for the ABCH and ARP is reported in Figure 2. In the most recent generation, i.e. 13^th, the Ne is 301 for the ARP and 443 for the ABCH depicting the good variability present in these breeds. In the Holstein breed, Makanjuola et al. (2020) reported a Ne varying from 43 to 66 according to the input information used to calculate it, i.e. SNP genotype vs. pedigree information.

Figure 2. Effective population size (Ne) per generation ago for the ARP and ABCH breeds.

Inbreeding

The inbreeding statistics are reported in Table 1. The average F_HOM coefficient is close to zero in both groups showing the lowest (–0.123) and highest (0.208) values for the ARP. The expected and observed heterozygosity are similar for both breeds. Values were for the ABCH and ARP respectively 0.3761 and 0.3721 for the observed heterozygosity and 0.3763 and 0.3710 for the expected one.

Table 1. Minimum (MIN), maximum (MAX), and average (AVG) values for inbreeding coefficients calculated for the Aosta population.

Inbreeding coefficients	ABCH			ARP
	MIN	MAX	AVG (sd)	MIN	MAX	AVG (sd)
FHOM	−0.078	0.169	0.001 (0.02)	−0.123	0.208	−0.003 (0.03)
FROH*	0.007	0.220	0.056 (0.02)	0.005	0.250	0.064 (0.02)
FROH class ≤2 Mbp	0.007	0.220	0,056 (0.02)	0.005	0.250	0,064 (0,02)
FROH class 2-4 Mbp	0.001	0.172	0,030 (0.02)	0.001	0.209	0,037 (0,02)
FROH class 4-8 Mbp	0.002	0.109	0,012 (0.01)	0.002	0.160	0,016 (0.01)
FROH class 8-16 Mbp	0.003	0.049	0,006 (0.01)	0.003	0.085	0,007 (0.01)
FROH class ≥16 Mbp	0.006	0.009	0,007 (0)	0.006	0.016	0,008 (0)

*F_ROH calculated considering five Runs of Homozygosity (ROH) length (in Mbp) classes.

F_ROH values calculated for each length class of ROH are very low especially when compared to values calculated with the same approach in a highly selected specialised breed, such as the Holstein (Makanjuola et al. 2020).

In their study Strillacci et al. (2020) reported very similar inbreeding coefficients to the ones obtained here: for the ARP (VPR in Strillacci et al. 2020) F_ROH was 0.067 compared to 0.064 in this study, while the F_HOM was −0.003 as in this study. More recently, Signer‐Hasler et al. (2023) calculated the F_ROH coefficient using the Illumina BovineHD bead chip for the Aosta Red Pied, Black Pied and Chestnut finding average values of 0.034, 0.029 and 0.025 respectively. The average F_ROH calculated in the Italian Holstein by Dadousis et al. (2022) is 0.15. It is interesting to note that, even if the ARP and the ABCH are two local breeds with a much smaller population size than the Holstein, the inbreeding coefficient is close to 1/3 with respect to the one found in the Holstein, where the intense selection is likely to be the cause of the loss of genetic variation.

The possibility of exploring at the genomic level the inbreeding in the Aosta cattle population in this study, allowed us to highlight the effectiveness of the selection programs applied for decades by farmers according to indications from the technical committee of their national breeder's association.

Runs of homozygosity

In all the cows of the two groups, ABCH and ARP, a total of 230,948 ROH were identified (Table 2). The largest total amount of ROH (156,711), the longest average length (2,125,789) and the highest average number (74) of ROH were found in the ARP. The obtained results showed that the Aosta population exhibits shorter ROH than the ones identified by Strillacci et al. (2020). This is particularly evident for the largest ROH that reaches a length up to 71,168,012 bp in the VRP (26,380,269 for ARP in this study) and 55,392,599 in the VBP (23,005,964 for ABCH in this study). Ferenčaković et al. (2013), concluded that the identification of ROH performed with the 50K SNP array is as reliable as the one obtained with the HD SNP chip: as such different SNP chip densities should not justify differences found in this study with respect to others as in Strillacci et al. (2020). We may speculate that the difference in ROH found in this study with respect to Strillacci et al. (2020) is likely due to the different composition of cohorts used in the two analyses.

Table 2. Minimum (MIN), maximum (MAX), and average number (AVG) and length (in bps) of Runs of Homozygosity (ROH) identified in the Aosta population.

	ROH statistics
	N				Length
Breed	Total	Min	Max	Avg (sd)	Min	Max	Avg
ABCH	74,237	11	179	68 (16)	1,000,004	23,005,964	2,023,899
ARP	156,711	9	183	74 (17)	1,000,000	26,380,269	2,125,789

*Average calculated per breed on total ROH length covered for each sample’s genome.

ROH_islands were assessed using SNP occurrences and are shown in Figure 3. A total of 6 and 4 ROH_islands were found in the ARP (threshold 51%) and in the ABCH (threshold 26%) respectively. In Table 3 the ROH_islands shared by more than 51% of cows, are reported together with the annotated genes and the associated QTL. The genes annotated in the ROH_islands identified in the ABCH breed are reported in Table 4.

Figure 3. Manhattan plots of the SNPs incidence. Red lines threshold represents the top 1% of SNP defining the ROH_islands.

Table 3. ROH_Islands identified in the ARP breed, with genes mapped in the ROH and associated QTL.

Breed	Chr	From	To	Genes	QTL*
ARP	2	61057756	61281214	CXCR4
ARP	2	61300454	61404353	DARS
ARP	5	18259786	18430163	KITLG	KITLG: Eye area pigmentation (QTL:21151)
ARP	6	69700321	69804132	PDGFRA
ARP	6	69826263	71305377	KIT, KDR, SRD5A3, TMEM165, CLOCK, PDCL2, NMU, EXOC1L, EXOC1	KIT: Eye area pigmentation (QTL:21154; QTL:21160), Somatic cell count (QTL:31634; QTL:31635); KDR: Bovine respiratory disease susceptibility (QTL:137210);
ARP	23	512811	1270657	KHDRBS2	KHDRBS2: Calving ease (maternal) (QTL:51867), Calving ease (QTL:51877), Dairy form (QTL:51868), Daughter pregnancy rate (QTL:51869), Foot angle (QTL:51870), Length of productive life (QTL:51874), Milk fat percentage (QTL:51871), Milk fat yield (QTL:51872), Milk protein percentage (QTL:51875), Net merit (QTL:51873), Rear leg placement - side view (QTL:51876), Somatic cell score (QTL:51878), Stillbirth (QTL:51879), Strength (QTL:51880), Teat length (QTL:51881)

*QTL (Quantitative Trait Loci) annotated in the Animal QTL Database for Cattle (https://www.animalgenome.org/cgi-bin/QTLdb/BT/index).

Table 4. ROH_Islands identified in the ABCH breed, with genes mapped in the ROH and associated QTL.

Breed	Chr	from	to	Genes	QTL*
ABCH	5	17,120,900	18,390,837	C5H12orf50, C5H12orf29, CEP290, TMTC3, KITLG	KITLG: Eye area pigmentation (QTL:21151)
ABCH	6	79,502,698	80,269,768	TECRL
ABCH	19	26,897,393	27,496,286	ASGR2, ASGR1, DLG4, ACADVL, DVL2, PHF23, GABARAP, CTDNEP1, ELP5, CLDN7, SLC2A4, YBX2, EIF5A, GPS2, NEURL4, ACAP1, KCTD11, TMEM95, TNK1, PLSCR3, TMEM256, NLGN2, SPEM1, SPEM2, TMEM102, FGF11, CHRNB1, ZBTB4, SLC35G6, POLR2A, TNFSF12, TNFSF13, SENP3, EIF4A1, CD68, MPDU1, SOX15, FXR2, SAT2, SHBG, ATP1B2, TP53, WRAP53, EFNB3, DNAH2	ACADVL: Bovine respiratory disease susceptibility (QTL:160277); ATP1B2: 305-day milk yield (QTL:14129; QTL:14128), Average daily milk yield (QTL:14134; QTL:14132; QTL:14133), Body temperature (QTL:14130), Milk fat percentage (QTL:14126), Milk protein percentage (QTL:14127), Red blood cell Na+/K + ATPase activity (QTL:14135; QTL:14136), Red blood cell potassium level (QTL:14131)
ABCH	21	497,917	2,006,266	MKRN3, MAGEL2, NDN, SNRPN, SNURF

*QTL (Quantitative Trait Loci) annotated in the Animal QTL Database for Cattle (https://www.animalgenome.org/cgi-bin/QTLdb/BT/index).

In accordance with the results of the ROH statistics, this study identified smaller ROH_islands respect to the results by Strillacci et al. (2020). The only exception is the new ROH_island on BTA21 found in the ABCH breed (Table 4). This region harbours five genes that are orthologs in many species, Homo sapiens included, where it has been identified as the causative region for the Prader-Willi syndrome (Cassidy 1997). A study, conducted on eight Italian beef cattle identified longer common ROH on the same chromosome shared between seven of them, starting at position 83,766 and ending at different positions (Fabbri et al. 2021). Furthermore, the genes annotated in this ROH_island, MKRN3, MAGEL2, NDN and SNRPN/SNURF, have a well-known role in the epigenetic regulation of precocious puberty onset, reproductive hormones synthesis, oocytes development, and pre- or post-implantation of embryos in cattle and humans (Alves et al. 2022). Lastly, in the study of Costilla et al. (2020), NDN, SNRPN and MAGEL2 genes have been associated with the temperament of cattle. Especially for the ABCH breed, farmers select cows for their fighting aptitude and the identified ROH_island may be related to a genetic component of the combativity of the Aosta cattle population.

On chromosome 5 the KITLG gene has been found in a ROH_island shared by both breeds; this gene has a QTL for the eye area pigmentation, and it has been associated with the roan phenotype in Belgian Blue and Shorthorn cattle (Seitz et al. 1999). Among the top 51% ROH_islands there are other genes that control coat colour such as the KIT gene that regulates the spotting (Fontanesi et al. 2010) and the PDGFRA gene that has been associated with white spotting (Fan et al. 2014), and also trypanotolerance in crossbred cattle (Yougbaré et al. 2021). Furthermore, the KDR gene has been associated with white pattern in Hereford cattle (Whitacre et al. 2013) and the TECRL gene associated with head pigmentation pattern (Mészáros et al. 2015), feed efficiency (de Almeida Santana et al. 2016) and fertility (Dias et al. 2017).

Among the other annotated genes in the identified ROH_islands, there are genes related to immunity, such as:

• CXCR4: that regulates the immune response (Revskij et al. 2022) and has been associated with BVDV persistently infected calves (Helal et al. 2013);

• CLDN7: that influences epithelial intestinal permeability. In intestine-specific claudin-7 knockout mice, this gene showed increased intestinal permeability and inflammation of mucosal structures (Meissner et al. 2017; Aschenbach et al. 2019);

• GPS2: that is involved in stress response and in immune system processes, particularly associated with bovine and ovine leukaemia virus response (Klener et al. 2006; Casas et al. 2020);

• TNFSF13 (also known as APRIL): is a proliferation-inducing ligand and is a member of the BAFF system molecules that play a vital role in mature B-cell survival and in the secretion of IgA antibodies. TNFSF13 can also have an adjuvant-like effect on the immune system to enhance antigen-specific humoral immunity (Zhang et al. 2010; Mallikarjunappa et al. 2019).

Furthermore, in the identified ROH_islands, in addition to the hereinbefore listed genes related to immunity, other interesting genes related to adaptability, feed efficiency and residual feed intake have been found and are reported in Table S1. Among those, very relevant for the Aosta cattle are: the ACADVL gene, associated with the mobilisation of fatty acids during periods of low energy intake and fatty acids’ beta-oxidation, resulting in a better adaptation to low feed periods (van Dorland and Bruckmaier 2013; Alaedin et al. 2021); the ACAP1 gene, involved in the cold stress response, it has been found to be an up-regulated gene in severe cold-exposed cattle (Xu et al. 2017); and the CLOCK gene that has a key role in positive modulation of the inflammatory response and regulates the circadian rhythm in mammals (Casey and Plaut 2012). All the genes found in the identified ROH_islands have been classified in Table S1 based on the available published studies.

The genes harboured in the ROH_islands of the Aosta population might be related to the ability of the Aosta breeds to adapt to harsh farming conditions. Particularly, the presence of adaptative genes in ROH_islands (Table S1) related to fatty acids mobilisation (e.g. ACADVL), cold stress response (e.g. ACAP1), and immune response (e.g. CXCR4, CLDN7, GPS2, TNFSF13) could be related to the summer pasture practice to which these breeds are undertaken.

Genotypic and gene frequencies for milk protein variant loci

Many studies investigated milk proteins and the effects of their variants on (i) the protein concentration and composition of milk (Heck et al. 2009; Huang et al. 2012), (ii) the milk coagulation properties (Jakob and Puhan 1992; Cipolat-Gotet et al. 2018); (iii) the effect on human nutrition (Tauzin et al. 2002; Caroli et al. 2004). Milk from Aosta cattle is used mainly to produce the PDO Fontina cheese. For this reason, the Aosta breeds’ selection plan pays attention to milk protein variants reporting for artificial insemination bulls their genotype for k-casein, β-casein and β-lactoglobulin.

Genotypic frequencies for the milk proteins loci, that are present on the 100K SNP chip, are reported in Table 5. For what concerns the casein fraction of the milk proteins: the k-casein registers a higher frequency of B allele with respect to A in the ARP (i.e. B allele frequency = 0.63) while in the ABCH the A allele is the one with the largest frequency, i.e. 0.60. The A variant of the k-casein has been associated in many studies with a lower k-casein content in milk, which results in a lower cheese yield (Heck et al. 2009; Huang et al. 2012).

Table 5. Genotypic frequencies for milk proteins and the lactoferrin loci.

Milk traits
Trait	Chr	Genotype	Frequency
			ABCH	ARP
α-S1-casein	6	AA	0.644	0.903
		AG	0.314	0.094
		GG	0.042	0.002
β-casein^#	6	A2A2	0,302	0,363
		A2B	0,270	0,023
		A2C	0,000	0,004
		A2E	0,000	0,000
		A2I	0,115	0,131
		A1A1	0,010	0,070
		A1A2	0,109	0,309
		A1B	0,061	0,013
		A1I	0,021	0,065
		BB	0,055	0,000
		BC	0,000	0,000
		BI	0,046	0,007
		II	0,011	0,014
β-lactoglobulin	11	AA	0.100	0.074
		AB	0.411	0.388
		BB	0.489	0.538
K-casein	6	AA	0.352	0.133
		AB	0.490	0.475
		BB	0.158	0.392
Lactoferrin*	22	GG	0.066	0.074
		AG	0.409	0.398
		AA	0.525	0.528

^#The β-casein gene variants have been identified considering the genotypes at the following SNPrs, as suggested by (Chessa et al. 2020): rs433954503, rs43703013, rs43703012, rs109299401, rs43703011, rs721259074 and rs3423226649.

*SNP markers on the GGP Bovine 100K.

Comparing the Aosta population with other autochthonous breeds, Sanchez et al. (2020) found that the most frequent variant in the Abondance and Vosgienne breeds was A variant, with frequencies of 64.8% and 65.5% respectively. On the contrary the B variant was found at higher frequencies in Brown Swiss, Montbéliarde, Tarentaise and Normande cattle, 74.4%, 61.6%, 56.7% and 85.5% respectively.

The β-casein variant has been recently brought to attention for its possible impact on the milk digestibility properties (Truswell 2005; Fernández-Rico et al. 2022). Even though the relationship with digestibility and gastrointestinal sensibility of A1 milk is still unclear, some milk packaging industries are nowadays marketing A2A2 milk to consumers because of the health claims that have been stated in the latest years. Recently a study on the Holstein Friesian cattle (Chessa et al. 2020) suggested that some selection is undergoing in this breed towards an increase of A2 allele in the last decades. In Anglo-Saxon countries, where the marketing of A2A2 milk was successful, sires with A2A2 genotype were in fact preferentially selected as reproducers. The use of these bulls as reproducers may partially explain the increased frequency of A2 variant that the authors reported, i.e. from 0.49 to 0.56 in the recent population. In our study the most frequent allele of the β-casein haplotype (CSN2 gene) was the A2 allele, being 0.60 in the ARP and 0.55 in the ABCH.

In the Pinzgauer alpine population, Caroli et al. (2010) found that the A1 and A2 alleles had a frequency 0.357 and 0.572 respectively, while the B variant is the least frequent one, f(B)=0.015. Some studies (Brooke-Taylor et al. 2017; Sebastiani et al. 2020) affirm that the A1 variant results in an improved curd consistency and milk coagulation as well as a better micelle size but results in a lower milk digestibility because of the bioactive peptide that forms during the milk digestion. On the other hand, Heck et al. (2009) found that Holstein cows with A1A1 genotype show a reduction in protein yield respect to those with A2A2 genotype, but this variant was also associated with a higher concentration of α-S1-casein and k-casein. They also identified that the A1A2 genotype resulted, respect to the homozygous A1A1 genotype, in a smaller content of α-lactalbumin, β-lactoglobulin and α-S1-casein and a greater content of β-casein.

For the α-S1-casein locus (CSN1S1 gene) the A mutation frequency was higher in both breeds being 0.95 in the ARP and 0.80 in the ABCH. The G mutation at the α-S1-CN SNP, is a very rare mutation also in the Holstein population (Heck et al. 2009). The G mutation was also found in a low frequency in the Italian Brown population by Mancini et al. (2013). Sanchez et al. (2020) reported a detailed scheme relating missense SNP mutation and protein variants. The mutation G for the SNP on chr 6 at position 85,427,427 (ARS_UCD1.2) is coded as C variant, while the A mutation in the B one. In the study of Caroli et al. (2010) the C variant in the Pinzgauer cattle had an allelic frequency of 0.27. The Pinzgauer is an alpine cattle population and has a frequency more similar to the ABCH (i.e. 0.2) respect to the Italian Holstein i.e. 0.003 (Chessa et al. 2020) and the Italian Brown, i.e. 0.08 (Mancini et al. 2013) which is more similar to the ARP having a frequency of 0.05.

Regarding the whey proteins, representing the 20% of milk proteins (Wood et al. 2021), the 100K SNP Chip releases the information for the β-lactoglobulin locus, i.e. BLG gene. This gene is particularly important in cheese making and the most frequent allele frequency in the Aosta population refers to allele B with 0.69 and 0.73 in ABCH and ARP breeds. The β-lactoglobulin B variant has been found as the most frequent variant also in the Abondance, Brown Swiss, Tarentaise, Normande and Vosgienne cattle with frequencies of 53.5, 68.8, 68.2, 57.0 and 68.8 respectively (Sanchez et al. 2020). A study of Heck et al. (2009) affirms that selecting for the B variant of the BLG gene will result in a higher concentration of caseins, but with the same milk composition and a different whey composition; for this reason, selection for the B variant of BLG gene can improve cheese production without negatively affect the cheese making properties.

However, the complexity of the relationship among the casein variants and the milk properties for cheese production make difficult to disclose an explicit selection for an overall combination of milk protein variants to exclusively enhance renneting properties.

In the Aosta breeds at the marker for the lactoferrin - LTF gene - the B allele has a frequency of 0.73 both for the ABCH and ARP (Table 5). The polymorphism at the LTF gene has been associated with milk Somatic Cell Count (SCC); in particular, Wojdak-Maksymiec (2016) found that the heterozygote genotype AG is associated with a higher SCC in milk while the AA has the lowest value.

Genotypic and gene frequencies for meat quality loci

In the last decades, as the consumers require a good quality product, the selection of beef cattle breeds has been also oriented to meat tenderness. The availability of the Calpain, Calpastatin and Leptin genes tests made possible to identify carriers of favourable meat quality characteristics before their slaughtering.

Since Aosta breeds are a double purpose cattle used to produce milk and meat, knowing the distribution of the different genotypes for the Calpain, Calpastatin and Leptin genes is important.

In Table 6 the genotypic frequencies of Calpain_316 and Calpain_4751 markers are reported. The G allele of Calpain_316 show frequencies of 0.92 and 0.93 in the ABCH and ARP breeds respectively. The most frequent genotype for Calpain_4751 is the heterozygote TC in the ABCH, while for the ARP the TT genotype outnumbers the others with a frequency of 0.546. Lisa and Di Stasio (2009) found allelic frequencies of 1 and 0.94 for the G allele at the CAPN_316 marker in the Aosta Black Pied and in the ARP respectively and for the CAPN_4751 frequencies of 0.62 and 0.61 for the T allele. The T allele frequency for the CAPN_4751 in our study is 0.56 for the ABCH and 0.74 for the ARP. The difference can be due to the effect of selection (ARP), to the analysis of Chestnut and Heren jointly, or to the smaller sample size of the study of Lisa and Di Stasio (2009).

Table 6. Genotypic frequencies at loci related to meat quality traits.

Meat traits
Trait	Chr	Genotype	Frequency
			ABCH	ARP
Calpain_316	29	GG	0.853	0.871
		GC	0.137	0.121
		CC	0.010	0.007
Calpain_4751	29	TT	0.307	0.546
		TC	0.502	0.385
		CC	0.190	0.069
CAST_2870*	7	GG	0.143	0.110
		AG	0.478	0.472
		AA	0.379	0.418
CAST_2959*	7	GG	0.005	0.037
		AG	0.105	0.296
		AA	0.890	0.667
Leptin 2FB	4	CC	0.340	0.588
		CT	0.489	0.361
		TT	0.171	0.051
UoGCAST1	7	GG	0.427	0.271
(CAST_282*)		GC	0.443	0.510
		CC	0.130	0.219

*SNP on the GGP Bovine 100K at the Calpastatin (CAST) gene.

For the Calpastatin (CAST) locus, allele frequencies varied across markers and breeds: the SNP CAST_2870 marker shows a frequency for the G allele of 0.38 (ABCH) and 0.35 (ARP); the SNP CAST_2959 marker shows a frequency for the G allele of 0.06 (ABCH) and 0.19 (ARP). Lastly, at the UoGCAST1 marker, the G allele have a frequency of 0.65 in the ABCH and 0.53 in the ARP.

The calpain gene, firstly identified by Smith et al. (2000), works in a complex system jointly with the calpastatin one, identified by Bishop et al. (1993). Their combined effect has an important role on meat tenderness. In fact, calpains are Ca²⁺ proteases that intervene during meat maturation processes, tenderising muscle fibres, while calpastatins are their inhibitor enzymes (Coria et al. 2018). It has been found that these proteins also influence flavour and juiciness (Casas et al. 2006) as well as colour, pH and water retention ability of meat (Reardon et al. 2010). Casas et al. (2006) found that the two homozygotes for the CAPN_4751 and 4959 markers have the most flavourful and juicy meat, with respect to the heterozygotes. Furthermore, Gill et al. (2009) reported that subjects with the CC genotype at the CAPN_316 marker show lower tenderometer values, higher hindquarter weights and tenderness scores by the taste panel, compared to the other genotypes. Many other studies identified the C allele as most favourable for meat tenderness at the markers: Calpain_316 (Page et al. 2002), Calpain_4751 (White et al. 2005) and UoGCAST1, also known as CAST_282 (Schenkel et al. 2006) (Casas et al. 2006; Gill et al. 2009; Chung et al. 2014). For the other two calpastatin gene makers (CAST_2870 and CAST_2959) identified by Cong et al. (1998) and Barendse (2002), the favourable allele is the A (Morris et al. 2006; Ribeca et al. 2009). Reardon et al. (2010) identified that at the UoGCAST1 marker, the GG genotype has a higher pH than the heterozygote and the CC. Making the Homozygote GG a good candidate for dark, firm, and dry beef (DFD) if the animals are not managed in the best way possible at the slaughterhouse.

Allele frequencies at the Leptin 2FB mutation of the leptin gene differ between the two populations: the most frequent allele is the C one found in the ARP (0.77). At the Leptin 2FB locus Schenkel et al. (2005) found that the C allele is associated with higher lean meat production with lower carcase fat. Brickell et al. (2010) found that primiparous heifers’ homozygote for the C allele showed a lower incidence of calf mortality at first calving.

The quality of meat (i.e. tenderness, juiciness, flavour), with respect to the genotypes at different markers with different combinations, has been tested by several authors in meat cattle breeds, but not specifically in the Aosta population and a targeted study could then be useful for the breeds valorisation and for their selection.

Genotypic and gene frequencies for fertility haplotypes/loci and mendelian inheritance diseases loci

In the Aosta population, the haplotypes affecting cows’ fertility and embryonic development found in cosmopolitan breeds (e.g. Brown Swiss Fertility Haplotypes - BHs, Holstein Fertility Haplotypes – HHs, Jerseys Fertility Haplotypes – JHs) are absent.

In Table 7 the genotypic frequencies for haplotypes and loci affecting fertility are reported. The variability identified in the Montbèliarde Haplotype 1 (MH1) could be the result of some ancestral cross between the Aosta breeds and the Montbèliarde.

Table 7. Genotypic frequencies for fertility traits.

Fertility
Trait^#	Chr	Genotype	Frequency
			ABCH	ARP
COQ9	18	AA	0.343	0.160
rs109301586*		AG	0.508	0.470
		GG	0.149	0.370
MH1	19	CC	0.958	0.997
		CT	0.042	0.003
		TT	0.000	0.000
PKP2-988*	5	DD	0.065	0.161
		DI	0.357	0.469
		II	0.578	0.370
STAT3_25402*	19	GG	0.034	0.018
		GT	0.283	0.258
		TT	0.683	0.724
STAT3_19069*	19	CC	0.027	0.019
		CT	0.247	0.254
		TT	0.726	0.727
STAT5_13244*	19	CC(GG)	0.059	0.020
(STAT5_13516*)		CT(GT)	0.339	0.262
		TT(TT)	0.602	0.718
STAT5_13319*	19	AA	0.001	0.000
		AG	0.080	0.004
		GG	0.919	0.996

*SNP markers on the GGP Bovine 100K.

^#Coenzyme Q9 (COQ9) (rs109301586); Montbéliarde Haplotype 1 (MH1); plakophilin 2 (PKP2_988). I → insertion; D → Deletion; Signal transducer and activator of transcription (STAT).

The COQ9-rs109301586 marker has been associated with improved fertility and higher oocyte mitochondrial content, with positive characteristics found for the A allele; this allele can improve oocyte competence needed for supporting the following embryonic development (Ortega et al. 2017).

Sugimoto et al. (2013) found that cows with the Del/Del genotype at the PKP2-988 marker had a higher conception rate than the TA/TA, for this reason the DD genotype is associated with more fertile cows than the II, in the Holstein population.

The Aosta breeds have also been tested for all the mendelian inheritance diseases markers on the 100K SNP Chip. In Table S2 only the markers with variability are reported and it is possible to see that the Aosta cattle are free from almost all the mendelian inherited diseases except for the presence of few heterozygote individuals.

For the Mulefoot_241 disease the recessive homozygous genotypes have been identified with a frequency of 0.004 and 0.001 for the ABCH and ARP, respectively.

The Infectious Bovine Keratoconjunctivitis (IBK) is a disease that has been found associated with the A/G mutation at the IBK marker on BTA8. The G mutation is the resistant allele with a substitution effect of 8.2% on reduction of pinkeye infection Kataria et al. (2011). In the Aosta population the G mutation has a frequency of 0.55 in the ABCH and 0.42 in the ARP.

Conclusions

The cost reduction in SNP genotyping favoured the genomic data production also in small local populations and in females. The Aosta breeds are part of the DUALBREEDING PSRN project envisaging a wide production of genotypes on the female population. The knowledge of the genomic characteristics of autochthonous populations is an important asset of information allowing a better understanding of their genetic background.

The autochthonous populations represent a unique source of biodiversity (Rognoni et al. 1983) and the genomic technologies made possible to accurately characterise them at the genomic level and to identify peculiarities linked to their extensive farming system and adaptation to the harsh mountain environment.

In this study the Aosta cattle female population has been analysed after the new Aosta breeds classification occurred in 2020. With this study and the one from Strillacci et al. (2020), it is possible to state that some ROH_islands are shared by bulls and cows. The genes annotated in the ROH_islands are related to efficiency, adaptability and resilience traits, including immune system. Indeed, we can suppose that the summer pasture practice favour individuals that are efficient in dealing with infections and have a prompt and fast response to stressors such as the harsh climatic conditions and the lower quantity of feed.

It is important to highlight that the Aosta breeds are free from the mendelian inheritance disorders identified in cosmopolitan breeds. The genotypic frequencies at the analysed loci well explain and highlight the selection that has occurred in these breeds, to improve milk and cheese production. Moreover, the ones related to genes influencing meat quality showed that there is possibility to improve the selection for this trait. Of particular interest are the calpain and calpastatin variants where the genotypic frequencies reflect the actual meat characteristics of the Aosta cattle. Since the favourable gene combinations are at low frequency leaving large room for improvement, it is possible to apply a gene assisted selection to improve the meat tenderness, a characteristic desired by consumers. Even though, good practices at the slaughterhouse and the use of a traditional long maturation periods can greatly improve the quality of meat; a gene assisted selection approach can be more extensively proposed and valorised since consumers are willing to pay more for good quality meat.

Author contributions

FB, MGS and AB: Conceived and designed the study; FB: data analysis; AD, RM and CP: participated to data analysis; MV and VB: data curation; MV and AB project administration; MV: funding acquisition; FB, MGS and AB: wrote the original draft. All authors reviewed and approved the final version of the manuscript.

Acknowledgments

The Associazione Nazionale Allevatori Razza Valdostana is kindly acknowledged for providing the genotypes and the data used for the analyses.

Disclosure statement

No potential conflict of interest was reported by the authors.

Data availability statement

The data supporting the conclusions of this manuscript are included in the Supplementary Materials.

Correction Statement

This article has been corrected with minor changes. These changes do not impact the academic content of the article.

Additional information

Funding

Project name: “Le razze bovine a duplice attitudine: un modello alternativo di zootecnia eco-sostenibile – DUAL BREEDING FASE 2” – PSRN Programma di Sviluppo Rurale Nazionale 2014/2020 – Sottomisura 10.2 – PERIODO 2020-2023 – funded by the European Agricultural Fund for Rural Development.