Assessment of milk metabolites as biomarkers for predicting feed efficiency in dairy sheep

Estimating feed efficiency (FE) in dairy sheep is challenging due to the high cost of systems that measure individual feed intake. Identifying proxies that can serve as effective predictors of FE could make it possible to introduce FE into breeding programs. Here, 39 Assaf ewes in first lactation were evaluated regarding their FE by 2 metrics, residual feed intake (RFI) and feed conversion ratio (FCR). The ewes were classified into high, medium and low groups for each metric. Milk samples of the 39 ewes were subjected to untargeted metabolomics analysis. The complete milk metabolomic signature was used to discriminate the FE groups using partial least squares discriminant analysis. A total of 41 and 26 features were selected as the most relevant features for the discrimination of RFI and FCR groups, respectively. The predictive ability when utilizing the complete milk metabolomic signature and the reduced data sets were investigated using 4 machine-learning algorithms and a multivariate regression method. The Orthogonal Partial Least Square algorithm outperformed other ML algorithms for the FCR prediction in the scenarios using the complete milk metabolite signature (r2 = 0.62 ± 0.06) and the 26 selected features (0.62 ± 0.15). Regarding RFI predictions, the scenarios using the 41 selected features outperformed the scenario with the complete milk metabolite signature, where the Multilayer feedforward artificial neural network (r2 = 0.18 ± 0.14) and extreme gradient boosting (r2 = 0.17 ± 0.15) outperformed other algorithms. The functionality of the selected metabolites implied that the metabolism of glucose, galactose, fructose, sphingolipids, amino acids, insulin, and thyroid hormones was at play. Compared with the use of traditional methods, practical applications of these biomarkers might simplify and reduce costs in selecting feed-efficient ewes.

Epigenetic regulation of functional candidate genes for milk production traits in dairy sheep subjected to protein restriction in the prepubertal stage

BACKGROUND

As the prepubertal stage is a crucial point for the proper development of the mammary gland and milk production, this study aims to evaluate how protein restriction at this stage can affect methylation marks in milk somatic cells. Here, 28 Assaf ewes were subjected to 42.3% nutritional protein restriction (14 animals, NPR) or fed standard diets (14 animals, C) during the prepubertal stage. During the second lactation, the milk somatic cells of these ewes were sampled, and the extracted DNA was subjected to whole-genome bisulfite sequencing.

RESULTS

A total of 1154 differentially methylated regions (DMRs) were identified between the NPR and C groups. Indeed, the results of functional enrichment analyses of the genes harboring these DMRs suggested their relevant effects on the development of the mammary gland and lipid metabolism in sheep. The additional analysis of the correlations of the mean methylation levels within these DMRs with fat, protein, and dry extract percentages in the milk and milk somatic cell counts suggested associations between several DMRs and milk production traits. However, there were no phenotypic differences in these traits between the NPR and C groups.

CONCLUSION

In light of the above, the results obtained in the current study might suggest potential candidate genes for the regulation of milk production traits in the sheep mammary gland. Further studies focusing on elucidating the genetic mechanisms affected by the identified DMRs may help to better understand the biological mechanisms modified in the mammary gland of dairy sheep as a response to nutritional challenges and their potential effects on milk production.

Transcriptome analysis of perirenal fat from Spanish Assaf suckling lamb carcasses showing different levels of kidney knob and channel fat

Introduction

Suckling lamb meat is highly appreciated in European Mediterranean countries because of its mild flavor and soft texture. In suckling lamb carcasses, perirenal and pelvic fat depots account for a large fraction of carcass fat accumulation, and their proportions are used as an indicator of carcass quality.

Material and Methods

This study aimed to characterize the genetic mechanisms that regulate fat deposition in suckling lambs by evaluating the transcriptomic differences between Spanish Assaf lambs with significantly different proportions of kidney knob and channel fat (KKCF) depots in their carcasses (4 High-KKCF lambs vs. 4 Low-KKCF lambs).

Results

The analyzed fat tissue showed overall dominant expression of white adipose tissue gene markers, although due to the young age of the animals (17-36 days), the expression of some brown adipose tissue gene markers (e.g., UCP1, CIDEA) was still identified. The transcriptomic comparison between the High-KKCF and Low-KKCF groups revealed a total of 80 differentially expressed genes (DEGs). The enrichment analysis of the 49 DEGs with increased expression levels in the Low-KKCF lambs identified significant terms linked to the biosynthesis of lipids and thermogenesis, which may be related to the higher expression of the UCP1 gene in this group. In contrast, the enrichment analysis of the 31 DEGs with increased expression in the High-KKCF lambs highlighted angiogenesis as a key biological process supported by the higher expression of some genes, such as VEGF-A and THBS1, which encode a major angiogenic factor and a large adhesive extracellular matrix glycoprotein, respectively.

Discussion

The increased expression of sestrins, which are negative regulators of the mTOR complex, suggests that the preadipocyte differentiation stage is being inhibited in the High-KKCF group in favor of adipose tissue expansion, in which vasculogenesis is an essential process. All of these results suggest that the fat depots of the High-KKCF animals are in a later stage of development than those of the Low-KKCF lambs. Further genomic studies based on larger sample sizes and complementary analyses, such as the identification of polymorphisms in the DEGs, should be designed to confirm these results and achieve a deeper understanding of the genetic mechanisms underlying fat deposition in suckling lambs.

Influence of a temporary restriction of dietary protein in prepubertal ewe lambs on first lactation milk traits and response to a mammary gland inflammatory challenge

This study evaluated the influence of a temporary nutritional protein restriction (NPR) performed, under commercial conditions, in prepubertal female lambs on first lactation milk production traits and the inflammatory response triggered by an inflammatory challenge of the. From 40 Assaf female lambs, we defined a control group (Cn = 20), which received a standard diet for replacement lambs and the NPR group (n = 20), which received the same diet but without soybean meal between 3 and 5 months of age. About 150 days after lambing, 24 of these ewes (13 NPR, 11C) were subjected to an intramammary infusion of E. coli lipopolysaccharide (LPS). Our dynamic study identified indicator traits of local (SCC) and systemic (rectal T^a, IL-6, CXCL8, IL-10, IL-36RA, VEGF-A) response to the LPS challenge. The NPR did not show significant effects on milk production traits and did not affect the SCC and rectal T^a after the LPS challenge. However, the NPR had a significant influence on 8 of the 14 plasma biomarkers analysed, in all the cases with higher relative values in the C group. The effects observed on VEGF-A (involved in vasculogenesis during mammary gland development and vascular permeability) and IL-10 (a regulatory cytokine classically known by its anti-inflammatory action) are the most remarkable to explain the differences found between groups. Whereas further studies should be undertaken to confirm these results, our findings are of interest considering the current concern about the future world's demand for protein and the need for animal production systems to evolve toward sustainability.

Feed efficiency in dairy sheep: An insight from the milk transcriptome

Introduction

As higher feed efficiency in dairy ruminants means a higher capability to transform feed nutrients into milk and milk components, differences in feed efficiency are expected to be partly linked to changes in the physiology of the mammary glands. Therefore, this study aimed to determine the biological functions and key regulatory genes associated with feed efficiency in dairy sheep using the milk somatic cell transcriptome.

Material and methods

RNA-Seq data from high (H-FE, n = 8) and low (L-FE, n = 8) feed efficiency ewes were compared through differential expression analysis (DEA) and sparse Partial Least Square-Discriminant analysis (sPLS-DA).

Results

In the DEA, 79 genes were identified as differentially expressed between both conditions, while the sPLS-DA identified 261 predictive genes [variable importance in projection (VIP) > 2] that discriminated H-FE and L-FE sheep.

Discussion

The DEA between sheep with divergent feed efficiency allowed the identification of genes associated with the immune system and stress in L-FE animals. In addition, the sPLS-DA approach revealed the importance of genes involved in cell division (e.g., KIF4A and PRC1) and cellular lipid metabolic process (e.g., LPL, SCD, GPAM, and ACOX3) for the H-FE sheep in the lactating mammary gland transcriptome. A set of discriminant genes, commonly identified by the two statistical approaches, was also detected, including some involved in cell proliferation (e.g., SESN2, KIF20A, or TOP2A) or encoding heat-shock proteins (HSPB1). These results provide novel insights into the biological basis of feed efficiency in dairy sheep, highlighting the informative potential of the mammary gland transcriptome as a target tissue and revealing the usefulness of combining univariate and multivariate analysis approaches to elucidate the molecular mechanisms controlling complex traits.

Intergenerational impact of dietary protein restriction in dairy ewes on epigenetic marks in the perirenal fat of their suckling lambs

In sheep, nutrition during the prepubertal stage is essential for growth performance and mammary gland development. However, the potential effects of nutrient restriction in a prepuberal stage over the progeny still need to be better understood. Here, the intergenerational effect of maternal protein restriction at prepubertal age (2 months of age) on methylation patterns was evaluated in the perirenal fat of Assaf suckling lambs. In total, 17 lambs from ewes subjected to dietary protein restriction (NPR group, 44% less protein) and 17 lambs from control ewes (C group) were analyzed. These lambs were ranked based on their carcass proportion of perirenal and cavitary fat and classified into HighPCF and LowPCF groups. The perirenal tissue from 4 NPR-LowPCF, 4 NPR-HighPCF, 4 C-LowPCF, and 4 C-HighPCF lambs was subjected to whole-genome bisulfite sequencing and differentially methylated regions (DMRs) were identified. Among other relevant processes, these DMRs were mapped in genes responsible for regulating the transition of brown to white adipose tissue and nonshivering thermoregulation, which might be associated with better adaptation/survival of lambs in the perinatal stage. The current study provides important biological insights about the intergenerational effect on the methylation pattern of an NPR in replacement ewes.

Elucidating genes and gene networks linked to individual susceptibility to milk fat depression in dairy goats

Dietary supplementation with marine lipids modulates ruminant milk composition toward a healthier fatty acid profile for consumers, but it also causes milk fat depression (MFD). Because the dairy goat industry is mainly oriented toward cheese manufacturing, MFD can elicit economic losses. There is large individual variation in animal susceptibility with goats more (RESPO+) or less (RESPO-) responsive to diet-induced MFD. Thus, we used RNA-Seq to examine gene expression profiles in mammary cells to elucidate mechanisms underlying MFD in goats and individual variation in the extent of diet-induced MFD. Differentially expression analyses (DEA) and weighted gene co-expression network analysis (WGCNA) of RNA-Seq data were used to study milk somatic cell transcriptome changes in goats consuming a diet supplemented with marine lipids. There were 45 differentially expressed genes (DEGs) between control (no-MFD, before diet-induced MFD) and MFD, and 18 between RESPO+ and RESPO-. Biological processes and pathways such as "RNA transcription" and "Chromatin modifying enzymes" were downregulated in MFD compared with controls. Regarding susceptibility to diet-induced MFD, we identified the "Triglyceride Biosynthesis" pathway upregulated in RESPO- goats. The WGCNA approach identified 9 significant functional modules related to milk fat production and one module to the fat yield decrease in diet-induced MFD. The onset of MFD in dairy goats is influenced by the downregulation of SREBF1, other transcription factors and chromatin-modifying enzymes. A list of DEGs between RESPO+ and RESPO- goats (e.g., DBI and GPD1), and a co-related gene network linked to the decrease in milk fat (ABCD3, FABP3, and PLIN2) was uncovered. Results suggest that alterations in fatty acid transport may play an important role in determining individual variation. These candidate genes should be further investigated.

Integrated analyses of the methylome and transcriptome to unravel sex differences in the perirenal fat from suckling lambs

In sheep, differences were observed regarding fat accumulation and fatty acid (FA) composition between males and females, which may impact the quality and organoleptic characteristics of the meat. The integration of different omics technologies is a relevant approach for investigating biological and genetic mechanisms associated with complex traits. Here, the perirenal tissue of six male and six female Assaf suckling lambs was evaluated using RNA sequencing and whole-genome bisulfite sequencing (WGBS). A multiomic discriminant analysis using multiblock (s)PLS-DA allowed the identification of 314 genes and 627 differentially methylated regions (within these genes), which perfectly discriminate between males and females. These candidate genes overlapped with previously reported QTLs for carcass fat volume and percentage of different FAs in milk and meat from sheep. Additionally, differentially coexpressed (DcoExp) modules of genes between males (nine) and females (three) were identified that harbour 22 of these selected genes. Interestingly, these DcoExp were significantly correlated with fat percentage in different deposits (renal, pelvic, subcutaneous and intramuscular) and were associated with relevant biological processes for adipogenesis, adipocyte differentiation, fat volume and FA composition. Consequently, these genes may potentially impact adiposity and meat quality traits in a sex-specific manner, such as juiciness, tenderness and flavour.

Microbial community in resistant and susceptible Churra sheep infected by Teladorsagia circumcincta

Gastrointestinal nematodes (GIN) are a major threat to health and welfare in small ruminants worldwide. Teladorsagia circumcincta is a nematode that inhabits the abomasum of sheep, especially in temperate regions, causing important economic losses. Given that T. circumcincta and microbiome share the same niche, interactions between them and the host are expected. Although it is known that within a sheep breed there are animals that are more resistant than others to infection by GIN, it is not known if the microbiome influences the phenotype of these animals. Under this condition, 12 sheep were classified according to their cumulative faecal egg count (cFEC) at the end of a first experimental infection, 6 as resistant group (RG) and 6 as susceptible group (SG) to T. circumcincta infection. Then, all sheep were experimentally infected with 70,000 L3 of T. circumcincta and at day 7 days post-infection were euthanized. At necropsy, gastric mucosa and gastric content from abomasum were collected to extract bacterial DNA and sequence V3-V4 region from 16S rRNA gene using Ilumina technology. After bioanalysis performed, results showed that α-diversity and β-diversity remained similar in both groups. However, resistant phenotype sheep showed a higher number of bacteria butyrate-fermenting species as Clostridium sensu stricto 1 (abundance in RG: 1.29% and in SG: 0.069%; p = 0.05), and Turicibacter (abundance in RG: 0.31% and in SG: 0.027%; p = 0.07) in gastric content but also Serratia spp in gastric mucosa (abundance in RG: 0.12% and in SG: 0.041%; p = 0.07). A trend towards a significant negative correlation between cFEC and Clostridium sensu stricto 1 abundance in gastric content was detected (r = - 0.537; p = 0.08). These data suggest that microbiome composition could be another factor associated with the development of the resistant phenotype modifying the interaction with the host and the in last instance affecting the individual risk of infection.

Low-density SNP panel for efficient imputation and genomic selection of milk production and technological traits in dairy sheep

The present study aimed to ascertain how different strategies for leveraging genomic information enhance the accuracy of estimated breeding values for milk and cheese-making traits and to evaluate the implementation of a low-density (LowD) SNP chip designed explicitly for that aim. Thus, milk samples from a total of 2,020 dairy ewes from 2 breeds (1,039 Spanish Assaf and 981 Churra) were collected and analyzed to determine 3 milk production and composition traits and 2 traits related to milk coagulation properties and cheese yield. The 2 studied populations were genotyped with a customized 50K Affymetrix SNP chip (Affymetrix Inc.) containing 55,627 SNP markers. The prediction accuracies were obtained using different multitrait methodologies, such as the BLUP model based on pedigree information, the genomic BLUP (GBLUP), and the BLUP at the SNP level (SNP-BLUP), which are based on genotypic data, and the single-step GBLUP (ssGBLUP), which combines both sources of information. All of these methods were analyzed by cross-validation, comparing predictions of the whole population with the test population sets. Additionally, we describe the design of a LowD SNP chip (3K) and its prediction accuracies through the different methods mentioned previously. Furthermore, the results obtained using the LowD SNP chip were compared with those based on the 50K SNP chip data sets. Finally, we conclude that implementing genomic selection through the ssGBLUP model in the current breeding programs would increase the accuracy of the estimated breeding values compared with the BLUP methodology in the Assaf (from 0.19 to 0.39) and Churra (from 0.27 to 0.44) dairy sheep populations. The LowD SNP chip is cost-effective and has proven to be an accurate tool for estimating genomic breeding values for milk and cheese-making traits, microsatellite imputation, and parentage verification. The results presented here suggest that the routine use of this LowD SNP chip could potentially increase the genetic gains of the breeding selection programs of the 2 Spanish dairy sheep breeds considered here.

Effect of genotyping strategies on the sustained benefit of single-step genomic BLUP over multiple generations

BACKGROUND

Single-step genomic best linear unbiased prediction (ssGBLUP) allows the inclusion of information from genotyped and ungenotyped individuals in a single analysis. This avoids the need to genotype all candidates with the potential benefit of reducing overall costs. The aim of this study was to assess the effect of genotyping strategies, the proportion of genotyped candidates and the genotyping criterion to rank candidates to be genotyped, when using ssGBLUP evaluation. A simulation study was carried out assuming selection over several discrete generations where a proportion of the candidates were genotyped and evaluation was done using ssGBLUP. The scenarios compared were: (i) three genotyping strategies defined by their protocol for choosing candidates to be genotyped (RANDOM: candidates were chosen at random; TOP: candidates with the best genotyping criterion were genotyped; and EXTREME: candidates with the best and worse criterion were genotyped); (ii) eight proportions of genotyped candidates (p); and (iii) two genotyping criteria to rank candidates to be genotyped (candidates' own phenotype or estimated breeding values). The criteria of the comparison were the cumulated gain and reliability of the genomic estimated breeding values (GEBV).

RESULTS

The genotyping strategy with the greatest cumulated gain was TOP followed by RANDOM, with EXTREME behaving as RANDOM at low p and as TOP with high p. However, the reliability of GEBV was higher with RANDOM than with TOP. This disparity between the trend of the gain and the reliability is due to the TOP scheme genotyping the candidates with the greater chances of being selected. The extra gain obtained with TOP increases when the accuracy of the selection criterion to rank candidates to be genotyped increases.

CONCLUSIONS

The best strategy to maximise genetic gain when only a proportion of the candidates are to be genotyped is TOP, since it prioritises the genotyping of candidates which are more likely to be selected. However, the strategy with the greatest GEBV reliability does not achieve the largest gain, thus reliability cannot be considered as an absolute and sufficient criterion for determining the scheme which maximises genetic gain.

Study on the concordance between different SNP-genotyping platforms in sheep

Different SNP genotyping technologies are commonly used in multiple studies to perform QTL detection, genotype imputation, and genomic predictions. Therefore, genotyping errors cannot be ignored, as they can reduce the accuracy of different procedures applied in genomic selection, such as genomic imputation, genomic predictions, and false-positive results in genome-wide association studies. Currently, whole-genome resequencing (WGR) also offers the potential for variant calling analysis and high-throughput genotyping. WGR might overshadow array-based genotyping technologies due to the larger amount and precision of the genomic information provided; however, its comparatively higher price per individual still limits its use in larger populations. Thus, the objective of this work was to evaluate the accuracy of the two most popular SNP-chip technologies, namely, Affymetrix and Illumina, for high-throughput genotyping in sheep considering high-coverage WGR datasets as references. Analyses were performed using two reference sheep genome assemblies, the popular Oar_v3.1 reference genome and the latest available version Oar_rambouillet_v1.0. Our results demonstrate that the genotypes from both platforms are suggested to have high concordance rates with the genotypes determined from reference WGR datasets (96.59% and 99.51% for Affymetrix and Illumina technologies, respectively). The concordance results provided in the current study can pinpoint low reproducible markers across multiple platforms used for sheep genotyping data. Comparing results using two reference genome assemblies also informs how genome assembly quality can influence genotype concordance rates among different genotyping platforms. Moreover, we describe an efficient pipeline to test the reliability of markers included in sheep SNP-chip panels against WGR datasets available on public databases. This pipeline may be helpful for discarding low-reliability markers before exploiting genomic information for gene mapping analyses or genomic prediction.

Genome-wide association studies (GWAS) and post-GWAS analyses for technological traits in Assaf and Churra dairy breeds

This study aimed to perform a GWAS to identify genomic regions associated with milk and cheese-making traits in Assaf and Churra dairy sheep breeds; second, it aimed to identify possible positional and functional candidate genes and their interactions through post-GWAS studies. For 2,020 dairy ewes from 2 breeds (1,039 Spanish Assaf and 981 Churra), milk samples were collected and analyzed to determine 6 milk production and composition traits and 6 traits related to milk coagulation properties and cheese yield. The genetic profiles of the ewes were obtained using a genotyping chip array that included 50,934 SNP markers. For both milk and cheese-making traits, separate single-breed GWAS were performed using GCTA software. The set of positional candidate genes identified via GWAS was subjected to guilt-by-association-based prioritization analysis with ToppGene software. Totals of 84 and 139 chromosome-wise significant associations for the 6 milk traits and the 6 cheese-making traits were identified in this study. No significant SNPs were found in common between the 2 studied breeds, possibly due to their genetic heterogeneity of the phenotypes under study. Additionally, 63 and 176 positional candidate genes were located in the genomic intervals defined as confidence regions in relation to the significant SNPs identified for the analyzed traits for Assaf and Churra breeds. After the functional prioritization analysis, 71 genes were identified as promising positional and functional candidate genes and proposed as targets of future research to identify putative causative variants in relation to the traits under examination. In addition, this multitrait study allowed us to identify variants that have a pleiotropic effect on both milk production and cheese-related traits. The incorporation of variants among the proposed functional and positional candidate genes into genomic selection strategies represent an interesting approach for achieving rapid genetic gains, specifically for those traits difficult to measure, such as cheese-making traits.

Identification of Regulatory Functions of LncRNAs Associated With T. circumcincta Infection in Adult Sheep

Several recent studies have demonstrated the role of long non-coding RNAs (lncRNAs) in regulating the defense mechanism against parasite infections, but no studies are available that investigated their relevance for immune response to nematode infection in sheep. Thus, the aim of the current study was to (i) detect putative lncRNAs that are expressed in the abomasal lymph node of adult sheep after an experimental infection with the gastrointestinal nematode (GIN) Teladorsagia circumcincta and (ii) to elucidate their potential functional role associated with the differential host immune response. We hypothesized that putative lncRNAs differentially expressed (DE) between samples from animals that differ in resistance to infection may play a significant regulatory role in response to nematode infection in adult sheep. To obtain further support for our hypothesis, we performed co-expression and functional gene enrichment analyses with the differentially expressed lncRNAs (DE lncRNAs). In a conservative approach, we included for this predictive analysis only those lncRNAs that are confirmed and supported by documentation of expression in gastrointestinal tissues in the current sheep gene atlas. We identified 9,105 putative lncRNA transcripts corresponding to 7,124 gene loci. Of these, 457 were differentially expressed lncRNA loci (DELs) with 683 lncRNA transcripts. Based on a gene co-expression analysis via weighted gene co-expression network analysis, 12 gene network modules (GNMs) were found significantly correlated with at least one of 10 selected target DE lncRNAs. Based on the principle of “guilt-by-association,” the DE genes from each of the three most significantly correlated GNMs were subjected to a gene enrichment analysis. The significant pathways associated with DE lncRNAs included ERK5 Signaling, SAPK/JNK Signaling, RhoGDI Signaling, EIF2 Signaling, Regulation of eIF4 and p70S6K Signaling and Oxidative Phosphorylation pathways. They belong to signaling pathway categories like Cellular Growth, Proliferation and Development, Cellular Stress and Injury, Intracellular and Second Messenger Signaling and Apoptosis. Overall, this lncRNA study conducted in adult sheep after GIN infection provided first insights into the potential functional role of lncRNAs in the differential host response to nematode infection.

A multiple-phenotype imputation procedure as a method for prediction of cheese-making efficiency in Spanish Assaf sheep

Sheep milk is mainly intended to manufacture a wide variety of high-quality cheeses. The ovine cheese industry would benefit from an improvement, through genetic selection, of traits related to the milk coagulation properties (MCPs) and cheese yield-related traits, broadly denoted as "cheese-making traits." Considering that routine measurements of these traits needed for genetic selection are expensive and time-consuming, this study aimed to evaluate the accuracy of a cheese-making phenotype imputation method based on the information from official milk control records combined with the pH of the milk. For this study, we analyzed records of milk production traits, milk composition traits, and measurements of cheese-making traits available from a total of 1,145 dairy ewes of the Spanish Assaf sheep breed. Cheese-making traits included five related to the MCPs and two cheese yield-related traits. The milk and cheese-making phenotypes were adjusted for significant effects based on a general linear model. The adjusted phenotypes were used to define a multiple-phenotype imputation procedure for the cheese-making traits based on multivariate normality and Markov chain Monte Carlo sampling. Five of the seven cheese-making traits considered in this study achieved a prediction accuracy of 0.60 computed as the correlation between the adjusted phenotypes and the imputed phenotypes. Particularly the logarithm of curd-firming time since rennet addition (logK20) (0.68), which has been previously suggested as a potential candidate trait to improve the cheese ability in this breed, and the logarithm of the ratio between the rennet clotting time and the curd firmness at 60 min (logRCT/A60) (0.65), which has been defined by other studies as an indicator trait of milk coagulation efficiency. This study represents a first step toward the possible use of the phenotype imputation of cheese-making traits to develop a practical methodology for the dairy sheep industry to impute cheese-making traits only based on the analysis of a milk sample without the need of pedigree information. This information could be also used in future planning of specific breeding programs considering the importance of the cheese-making efficiency in dairy sheep and highlights the potential of phenotype imputation to leverage sample size on expensive, hard-to-measure phenotypes.

Estimation of genetic parameters for cheese-making traits in Spanish Churra sheep

The global production of sheep milk is growing, and the main industrial use of sheep milk is cheese making. The Spanish Churra sheep breed is one of the most important native dairy breeds in Spain. The present study aimed to estimate genetic parameters for a wide range of traits influencing the cheese-making ability of Churra sheep milk. Using a total of 1,049 Churra ewes, we studied the following cheese-making traits: 4 traits related to milk coagulation properties (rennet coagulation time, curd-firming time, and curd firmness at 30 and 60 min after addition of rennet), 2 traits related to cheese yield (individual laboratory cheese yield and individual laboratory dried curd yield), and 3 traits measuring curd firmness over time (maximum curd firmness, time to attain maximum curd firmness, and syneresis). In addition, a list of milk traits, including the native pH of the milk and several milk production and composition traits (milk yield; the fat, protein, and dried extract percentages; and the somatic cell count), were also analyzed for the studied animals. After discarding the noncoagulating samples (only 3.7%), data of 1,010 ewes were analyzed with multiple-trait animal models by using the restricted maximum likelihood method to estimate (co)variance components, heritabilities, and genetic correlations. In general, the heritability estimates were low to moderate, ranging from 0.08 (for the individual laboratory dried curd yield trait) to 0.42 (for the fat percentage trait). High genetic correlations were found within pairs of related traits (i.e., 0.93 between fat and dried extract percentages, -0.93 between the log of the curd-firming time and curd firmness at 30 min, 0.70 between individual laboratory cheese yield and individual laboratory dried curd yield, and -0.94 between time to attain maximum curd firmness and syneresis). Considering all the information provided here, we suggest that in addition to the current consideration of the protein percentage trait for improving cheese yield traits, the inclusion of the pH of milk as a measured trait in the Churra dairy breeding program would represent an efficient strategy for improving the cheese-making ability of milk from this breed.

Differences within Churra breed sheep in the early immune response to the infection by Teladorsagia circumcincta

This study describes early immunological mechanisms that underlie resistance to Teladorsagia circumcincta infection in adult Churra sheep. After a first experimental infection, 6 animals were classified as resistant (RG) and 6 as susceptible (SG) to T. circumcincta infection based on their cumulative faecal egg count (cFEC) at the end of the infection. RG showed higher IgA levels against somatic antigen of T. circumcincta fourth-larvae stage (L4) in serum at day 3 post-infection (pi) (p < 0.05) and close to significance at day 21 pi (p = 0.06). Moreover, a strong negative correlation between cFEC and specific IgA was only significant in RG at day 3 pi (r = - 0.870; p < 0.05), but absent in SG. At the end of this infection, sheep were treated with moxidectin and infected again 3 weeks later to be slaughtered at day 7 pi. At necropsy, the specific IgA levels in gastric mucosa were similar between groups; the absence differences at day 7 pi could be due to a previous increase in the IgA response, probably around day 3 pi, as described during the first infection. L4 burden, 68% lower in RG than in SG, was influenced by the specific IgA in gastric mucus and the number of γδ T cells. RG group showed a positive correlation between γδ T cells and eosinophils (r = 0.900; p = 0.037); however, this correlation was not found in SG. These results show that these two phenotypes show different early immune response pattern to T. circumcincta infection in Churra sheep.

The Milk Microbiota of the Spanish Churra Sheep Breed: New Insights into the Complexity of the Milk Microbiome of Dairy Species

Simple Summary

In the last decade, the number of studies focused on the study of the microbiota of different tissues, organs, or physiological fluids has considerably increased. The milk of dairy species is an important and continuous source of commensal, mutualistic and potentially probiotic bacteria. Second-generation sequencing technologies have been applied to characterise the milk microbiota of dairy cows, whereas the study of the sheep milk microbiota is scarce. In the present study, we aimed to explore the bacterial diversity and composition of milk samples from the Churra sheep breed, a rustic autochthonous breed from the region of Castilla y León (Spain). Moreover, this study tries to clarify the complex bacterial composition of sheep milk comparing the results presented here with previous research on the milk microbiota of the Assaf sheep breed. This assessment has shown that the milk microbiota of ewes from one flock of the Assaf breed is more diverse than the milk microbiota reported here for two different flocks of Churra sheep. The study also provides a step into a better understanding of the link between the bacterial milk composition in these two sheep breeds and somatic cell count, an indicator trait of subclinical mastitis resistance in dairy sheep.

Abstract

Milk from healthy animals has classically been considered a sterile fluid. With the development of massively parallel sequencing and its application to the study of the microbiome of different body fluids, milk microbiota has been documented in several animal species. In this study, the main objective of this work was to access bacterial profiles of healthy milk samples using the next-generation sequencing of amplicons from the 16S rRNA gene to characterise the milk microbiome of the Churra breed. A total of 212 samples were collected from two Churra dairy farms with a different management system. The core milk microbiota in Churra ewes includes lesser genera (only two taxa: Staphylococcus and Escherichia/Shigella) than studies reported in other dairy species or even in a previous study in Assaf sheep milk. We found that diversity values in the two flocks of Churra breed were lower than the diversity of the milk microbiota in Assaf. The non-metric multidimensional scaling (NMDS) ordination using Bray-Curtis distance separates samples based on their microbiota composition. The information reported here might be used to understand the complex issue of milk microbiota composition.

Gene Networks Driving Genetic Variation in Milk and Cheese-Making Traits of Spanish Assaf Sheep

Most of the milk produced by sheep is used for the production of high-quality cheese. Consequently, traits related to milk coagulation properties and cheese yield are economically important to the Spanish dairy industry. The present study aims to identify candidate genes and their regulators related to 14 milk and cheese-making traits and to develop a low-density panel of markers that could be used to predict an individual's genetic potential for cheese-making efficiency. In this study, we performed a combination of the classical genome-wide association study (GWAS) with a stepwise regression method and a pleiotropy analysis to determine the best combination of the variants located within the confidence intervals of the potential candidate genes that may explain the greatest genetic variance for milk and cheese-making traits. Two gene networks related to milk and cheese-making traits were created using the genomic relationship matrices built through a stepwise multiple regression approach. Several co-associated genes in these networks are involved in biological processes previously found to be associated with milk synthesis and cheese-making efficiency. The methodology applied in this study enabled the selection of a co-association network comprised of 374 variants located in the surrounding of genes showing a potential influence on milk synthesis and cheese-making efficiency.

Identification of potential functional variants underlying ovine resistance to gastrointestinal nematode infection by using RNA-Seq

In dairy sheep flocks from Mediterranean countries, replacement and adult ewes are the animals most affected by gastrointestinal nematode (GIN) infections. In this study, we have exploited the information derived from an RNA-Seq experiment with the aim of identifying potential causal mutations related to GIN resistance in sheep. Considering the RNA-Seq samples from 12 ewes previously classified as six resistant and six susceptible animals to experimental infection by Teladorsagia circumcincta, we performed a variant calling analysis pipeline using two different types of software, gatk version 3.7 and Samtools version 1.4. The variants commonly identified by the two packages (high-quality variants) within two types of target regions - (i) QTL regions previously reported in sheep for parasite resistance based on SNP-chip or sequencing technology studies and (ii) functional candidate genes selected from gene expression studies related to GIN resistance in sheep - were further characterised to identify mutations with a potential functional impact. Among the genes harbouring these potential functional variants (930 and 553 respectively for the two types of regions), we identified 111 immune-related genes in the QTL regions and 132 immune-related genes from the initially selected candidate genes. For these immune-related genes harbouring potential functional variants, the enrichment analyses performed highlighted significant GO terms related to apoptosis, adhesion and inflammatory response, in relation to the QTL related variants, and significant disease-related terms such as inflammation, adhesion and necrosis, in relation to the initial candidate gene list. Overall, the study provides a valuable list of potential causal mutations that could be considered as candidate causal mutations in relation to GIN resistance in sheep. Future studies should assess the role of these suggested mutations with the aim of identifying genetic markers that could be directly implemented in sheep breeding programmes considering not only production traits, but also functional traits such as resistance to GIN infections.

Conjugated linoleic acid (CLA)-induced milk fat depression: application of RNA-Seq technology to elucidate mammary gene regulation in dairy ewes

Milk fat depression (MFD) is characterized by a reduction in the content of milk fat, presumably caused by the anti-lipogenic effects of rumen biohydrogenation intermediates, such as trans-10 cis-12 conjugated linoleic acid (CLA). In this study, RNA-Seq technology was used to help elucidate the mammary responses involved in CLA-induced MFD in lactating ewes. To this end, we compared the milk somatic cell transcriptome of ewes suffering from CLA-induced MFD with control ewes (i.e., those without MFD), as well as with ewes fed a diet supplemented with fish oil (FO-MFD) that we previously reported affects the mammary transcriptome. In the differential expression analysis between CLA-MFD and controls, we identified 1,524 differentially expressed genes (DEGs), whereas 653 were detected between CLA- and FO-MFD groups. Although this article focuses on lipid metabolism, CLA affected the expression of many genes related to other biological processes, especially immunity. Among the 55 genes shared by both MFD conditions, some genes linked to fatty acid synthesis, such as ACACA, AACS, ACSS2, or ACSS3, were downregulated. In addition, this study provides a list of candidate genes that are not usually considered in the nutrigenomics of MFD but that may act as key regulators of this syndrome in dairy ewes.

Early adipose deposits in sheep: comparative analysis of the perirenal fat transcriptome of Assaf and Churra suckling lambs

Adipose deposits influence the quality of ruminant carcasses, and in suckling lambs, internal types of adipose deposits represent a notable proportion of total fat. The aim of this study was to perform a comparative analysis of the perirenal fat transcriptomes of suckling lambs from two breeds with different growth and carcass characteristics. The perirenal fat tissue from 14 suckling lambs (Assaf, n = 8; Churra, n = 6) was used for the RNA-seq analysis. The functional enrichment analysis of the 670 highly expressed genes (>150 fragments per kilobase of exon per million fragments mapped) in the perirenal fat transcriptome of both breeds revealed that the majority of these genes were involved in energy processes. The expression of the UCP1 gene, a classical biomarker of brown fat, and the presence of multilocular adipocytes in the two breeds supported the presence of brown fat at the transition stage towards white fat tissue. The differential expression analysis performed identified 373 differentially expressed genes (DEGs) between the two compared breeds. Brown/white fat gene biomarkers were not included in the list of DEGs. In Assaf lambs, DEGs were enriched in Gene Ontology (GO) biological processes related to fatty-acid oxidation, whereas in Churra lambs, the majority of the significantly enriched GO terms were related to cholesterol synthesis, which suggests that upregulated DEGs in Assaf lambs are implicated in fat burning, whereas the Churra upregulated DEGs are linked to fat accumulation. These results can help to increase knowledge of the genes controlling early fat deposition in ruminants and shed light on fundamental aspects of adipose tissue growth.

Detection of quantitative trait loci and putative causal variants affecting somatic cell score in dairy sheep by using a 50K SNP chip and whole-genome sequencing

This study presents a scan of the ovine genome to identify quantitative trait loci (QTL) influencing the somatic cell score (SCS), a classical indicator of subclinical mastitis in sheep, and a subsequent high-resolution analysis of one of the identified QTL regions based on the analysis of whole-genome sequence data sets. A half-sib commercial population of Churra sheep genotyped with a 50K SNP chip was analyzed using linkage analysis (LA) and combined linkage and linkage disequilibrium analysis (LDLA). By LA, 2 5% chromosome-wide significant QTL on OAR5 and OAR25 and one 5% genome-wide significant QTL on ovine chromosome 20 (OAR20) were detected, whereas 22 significant associations were identified by LDLA. Two of the associations detected by LDLA replicated LA-detected effects (OAR20, OAR25). We compared the detected associations with previously reported QTL in sheep and cattle, and functional candidate genes were identified within the estimated confidence intervals. We then performed a high-resolution analysis of the OAR20 QTL region, the most significant QTL region identified by LA that replicated a QTL previously described in Churra sheep for SCS using microsatellite markers. For that, 2 segregating trios of 2 segregating families for the OAR20 QTL (each including the Qq sire and 2 daughters, QQ and qq) were selected for whole-genome sequencing. The bioinformatic analysis of the 6 sequenced samples performed across the genomic interval considered (14.2-41.7 Mb) identified a total of 227,030 variants commonly identified by 2 independent software packages. For the 3 different concordance tests considered, due to discrepancies regarding the QTL peak in the segregating families, the list of mutations concordant with the QTL segregating pattern was processed to identify the variants identified in immune-related genes that show a moderate/high impact on the encoded protein function. Among a list of 85 missense variants concordant with the QTL segregation pattern that were within candidate immune-related genes, 13 variants distributed across 7 genes [PKHD1, NOTCH4, AGER, ENSOARG00000009395 (HLA-C, Homo sapiens), ENSOARG00000015002 (HLA-B, H. sapiens), MOG, and ENSOARG00000018075 (BoLA, Bos taurus, orthologous to human HLA-A] were predicted to cause deleterious effects on protein function. Future studies should assess the possible associations of the candidate variants identified herein in commercial populations with indicator traits of udder inflammation (SCS, clinical mastitis).

Exploring the mechanisms of resistance to Teladorsagia circumcincta infection in sheep through transcriptome analysis of abomasal mucosa and abomasal lymph nodes

The present study exploited the RNA-seq technology to analyze the transcriptome of target tissues affected by the Teladorsagia circumcincta infection in two groups of adult ewes showing different statuses against gastrointestinal nematode (GIN) infection with the aim of identifying genes linked to GIN infection resistance in sheep. For this, based on the accumulated faecal egg count of 18 adult Churra ewes subjected to a first experimental infection with T. circumcincta, six ewes were classified as resistant and six others as susceptible to the infection. These 12 animals were dewormed and infected again. After humanitarian sacrifice of these 12 animals at day 7 post-infection, RNA samples were obtained from abomasal mucosa and lymph node tissues and RNA-Seq datasets were generated using an Illumina HiSeq 2000 sequencer. The distribution of the genes based on their expression level were very similar among the two different tissues and conditions. The differential expression analysis performed with two software (DESeq and EdgeR) only identified common differentially expressed genes (DEGs), a total of 106, in the lymph node samples which were considered as GIN-activated. The enrichment analysis performed for these GIN-activated genes identified some pathways related to cytokine-mediated immune response and the PPARG signaling pathway as well as disease terms related to inflammation and gastro-intestinal diseases as enriched. A systematic comparison with the results of previous studies confirmed the involvement of genes such as ITLN2, CLAC1 and galectins, in the immune mechanism activated against T. circumcincta in resistant sheep.

High-resolution analysis of selection sweeps identified between fine-wool Merino and coarse-wool Churra sheep breeds

Background

With the aim of identifying selection signals in three Merino sheep lines that are highly specialized for fine wool production (Australian Industry Merino, Australian Merino and Australian Poll Merino) and considering that these lines have been subjected to selection not only for wool traits but also for growth and carcass traits and parasite resistance, we contrasted the OvineSNP50 BeadChip (50 K-chip) pooled genotypes of these Merino lines with the genotypes of a coarse-wool breed, phylogenetically related breed, Spanish Churra dairy sheep. Genome re-sequencing datasets of the two breeds were analyzed to further explore the genetic variation of the regions initially identified as putative selection signals.

Results

Based on the 50 K-chip genotypes, we used the overlapping selection signals (SS) identified by four selection sweep mapping analyses (that detect genetic differentiation, reduced heterozygosity and patterns of haplotype diversity) to define 18 convergence candidate regions (CCR), five associated with positive selection in Australian Merino and the remainder indicating positive selection in Churra. Subsequent analysis of whole-genome sequences from 15 Churra and 13 Merino samples identified 142,400 genetic variants (139,745 bi-allelic SNPs and 2655 indels) within the 18 defined CCR. Annotation of 1291 variants that were significantly associated with breed identity between Churra and Merino samples identified 257 intragenic variants that caused 296 functional annotation variants, 275 of which were located across 31 coding genes. Among these, four synonymous and four missense variants (NPR2_His847Arg, NCAPG_Ser585Phe, LCORL_Asp1214Glu and LCORL_Ile1441Leu) were included.

Conclusions

Here, we report the mapping and genetic variation of 18 selection signatures that were identified between Australian Merino and Spanish Churra sheep breeds, which were validated by an additional contrast between Spanish Merino and Churra genotypes. Analysis of whole-genome sequencing datasets allowed us to identify divergent variants that may be viewed as candidates involved in the phenotypic differences for wool, growth and meat production/quality traits between the breeds analyzed. The four missense variants located in the NPR2, NCAPG and LCORL genes may be related to selection sweep regions previously identified and various QTL reported in sheep in relation to growth traits and carcass composition.

Electronic supplementary material

The online version of this article (doi:10.1186/s12711-017-0354-x) contains supplementary material, which is available to authorized users.

Elucidating fish oil-induced milk fat depression in dairy sheep: Milk somatic cell transcriptome analysis

In this study, RNA sequencing was used to obtain a comprehensive profile of the transcriptomic changes occurring in the mammary gland of lactating sheep suffering from fish oil-induced milk fat depression (FO-MFD). The milk somatic cell transcriptome analysis of four control and four FO-MFD ewes generated an average of 42 million paired-end reads per sample. In both conditions, less than 220 genes constitute approximately 89% of the total counts. These genes, which are considered as core genes, were mainly involved in cytoplasmic ribosomal proteins and electron transport chain pathways. In total, 117 genes were upregulated, and 96 genes were downregulated in FO-MFD samples. Functional analysis of the latter indicated a downregulation of genes involved in the SREBP signaling pathway (e.g., ACACA, ACSL, and ACSS) and Gene Ontology terms related to lipid metabolism and lipid biosynthetic processes. Integrated interpretation of upregulated genes indicated enrichment in genes encoding plasma membrane proteins and proteins regulating protein kinase activity. Overall, our results indicate that FO-MFD is associated with the downregulation of key genes involved in the mammary lipogenesis process. In addition, the results also suggest that this syndrome may be related to upregulation of other genes implicated in signal transduction and codification of transcription factors.

Variant discovery in the sheep milk transcriptome using RNA sequencing

BACKGROUND

The identification of genetic variation underlying desired phenotypes is one of the main challenges of current livestock genetic research. High-throughput transcriptome sequencing (RNA-Seq) offers new opportunities for the detection of transcriptome variants (SNPs and short indels) in different tissues and species. In this study, we used RNA-Seq on Milk Sheep Somatic Cells (MSCs) with the goal of characterizing the genetic variation within the coding regions of the milk transcriptome in Churra and Assaf sheep, two common dairy sheep breeds farmed in Spain.

RESULTS

A total of 216,637 variants were detected in the MSCs transcriptome of the eight ewes analyzed. Among them, a total of 57,795 variants were detected in the regions harboring Quantitative Trait Loci (QTL) for milk yield, protein percentage and fat percentage, of which 21.44% were novel variants. Among the total variants detected, 561 (2.52%) and 1,649 (7.42%) were predicted to produce high or moderate impact changes in the corresponding transcriptional unit, respectively. In the functional enrichment analysis of the genes positioned within selected QTL regions harboring novel relevant functional variants (high and moderate impact), the KEGG pathway with the highest enrichment was "protein processing in endoplasmic reticulum". Additionally, a total of 504 and 1,063 variants were identified in the genes encoding principal milk proteins and molecules involved in the lipid metabolism, respectively. Of these variants, 20 mutations were found to have putative relevant effects on the encoded proteins.

CONCLUSIONS

We present herein the first transcriptomic approach aimed at identifying genetic variants of the genes expressed in the lactating mammary gland of sheep. Through the transcriptome analysis of variability within regions harboring QTL for milk yield, protein percentage and fat percentage, we have found several pathways and genes that harbor mutations that could affect dairy production traits. Moreover, remarkable variants were also found in candidate genes coding for major milk proteins and proteins related to milk fat metabolism. Several of the SNPs found in this study could be included as suitable markers in genotyping platforms or custom SNP arrays to perform association analyses in commercial populations and apply genomic selection protocols in the dairy production industry.

Implementation of an extended ZINB model in the study of low levels of natural gastrointestinal nematode infections in adult sheep

Background

In this study, two traits related with resistance to gastrointestinal nematodes (GIN) were measured in 529 adult sheep: faecal egg count (FEC) and activity of immunoglobulin A in plasma (IgA). In dry years, FEC can be very low in semi-extensive systems, such as the one studied here, which makes identifying animals that are resistant or susceptible to infection a difficult task. A zero inflated negative binomial model (ZINB) model was used to calculate the extent of zero inflation for FEC; the model was extended to include information from the IgA responses.

Results

In this dataset, 64 % of animals had zero FEC while the ZINB model suggested that 38 % of sheep had not been recently infected with GIN. Therefore 26 % of sheep were predicted to be infected animals with egg counts that were zero or below the detection limit and likely to be relatively resistant to nematode infection. IgA activities of all animals were then used to decide which of the sheep with zero egg counts had been exposed and which sheep had not been recently exposed. Animals with zero FEC and high IgA activity were considered resistant while animals with zero FEC and low IgA activity were considered as not recently infected. For the animals considered as exposed to the infection, the correlations among the studied traits were estimated, and the influence of these traits on the discrimination between unexposed and infected animals was assessed.

Conclusions

The model presented here improved the detection of infected animals with zero FEC. The correlations calculated here will be useful in the development of a reliable index of GIN resistance that could be of assistance for the study of host resistance in studies based on natural infection, especially in adult sheep, and also the design of breeding programs aimed at increasing resistance to parasites.

Electronic supplementary material

The online version of this article (doi:10.1186/s12917-016-0723-7) contains supplementary material, which is available to authorized users.

Transcriptome expression analysis of candidate milk genes affecting cheese-related traits in 2 sheep breeds

Because ewe milk is principally used for cheese making, its quality is related to its content of total solids and the way in which milk constituents influence cheese yield and determine the technological and organoleptic characteristics of dairy products. Therefore, an in-depth knowledge of the expression levels of milk genes influencing cheese-related traits is essential. In the present study, the milk transcriptome data set of 2 dairy sheep breeds, Assaf and Spanish Churra, was used to evaluate the expression levels of 77 transcripts related to cheese yield and quality traits. For the comparison between both breeds, we selected the RNA sequencing (RNA-Seq) data at d 10 of lactation because this is the time point at which within and between breed differences due to lactation length are minimal. The evaluated genes encode major milk proteins (caseins and whey proteins), endogenous proteases, and enzymes related to fatty acid metabolism and citrate content. Through this analysis, we identified the genes predominantly expressed in each of the analyzed pathways that appear to be key genes for traits related to sheep milk cheese. Among the highly expressed genes in both breeds were the genes encoding caseins and whey proteins (CSN2, CSN3, CSN1S1, ENSOARG00000005099/PAEP, CSN1S2, LALBA), genes related to lipid metabolism (BTN1A1, XDH, FASN, ADFP, SCD, H-FABP, ACSS2), and one endogenous protease (CTSB). Moreover, a differential expression analysis between Churra and Assaf sheep allowed us to identify 7 genes that are significantly differentially expressed between the 2 breeds. These genes were mainly linked to endogenous protease activity (CTSL, CTSK, KLK10, KLK6, SERPINE2). Additionally, there were 2 differentially expressed genes coding for an intracellular fatty acid transporter (FABP4), an intermediate molecule of the citric acid cycle (SUCNR1), and 2 heat shock proteins (HSP70, HSPB8) that could be related to high protein production. The differential expression of these genes could have a direct influence on the different phenotypes observed between the 2 analyzed breeds.

Detection and replication of QTL underlying resistance to gastrointestinal nematodes in adult sheep using the ovine 50K SNP array

BACKGROUND

Persistence of gastrointestinal nematode (GIN) infection and the related control methods have major impacts on the sheep industry worldwide. Based on the information generated with the Illumina OvineSNP50 BeadChip (50 K chip), this study aims at confirming quantitative trait loci (QTL) that were previously identified by microsatellite-based genome scans and identifying new QTL and allelic variants that are associated with indicator traits of parasite resistance in adult sheep. We used a commercial half-sib population of 518 Spanish Churra ewes with available data for fecal egg counts (FEC) and serum levels of immunoglobulin A (IgA) to perform different genome scan QTL mapping analyses based on classical linkage analysis (LA), a combined linkage disequilibrium and linkage analysis (LDLA) and a genome-wide association study (GWAS).

RESULTS

For the FEC and IgA traits, we detected a total of three 5 % chromosome-wise significant QTL by LA and 63 significant regions by LDLA, of which 13 reached the 5 % genome-wise significance level. The GWAS also revealed 10 significant SNPs associated with IgAt, although no significant associations were found for LFEC. Some of the significant QTL for LFEC that were detected by LA and LDLA on OAR6 overlapped with a highly significant QTL that was previously detected in a different half-sib population of Churra sheep. In addition, several new QTL and SNP associations were identified, some of which show correspondence with effects that were reported for different populations of young sheep. Other significant associations that did not coincide with previously reported associations could be related to the specific immune response of adult animals.

DISCUSSION

Our results replicate a FEC-related QTL located on OAR6 that was previously reported in Churra sheep and provide support for future research on the identification of the allelic variant that underlies this QTL. The small proportion of genetic variance explained by the detected QTL and the large number of functional candidate genes identified here are consistent with the hypothesis that GIN resistance/susceptibility is a complex trait that is not determined by individual genes acting alone but rather by complex multi-gene interactions. Future studies that combine genomic variation analysis and functional genomic information may help elucidate the biology of GIN disease resistance in sheep.

An interpretive review of selective sweep studies in Bos taurus cattle populations: identification of unique and shared selection signals across breeds

This review compiles the results of 21 genomic studies of European Bos taurus breeds and thus provides a general picture of the selection signatures in taurine cattle identified by genome-wide selection-mapping scans. By performing a comprehensive summary of the results reported in the literature, we compiled a list of 1049 selection sweeps described across 37 cattle breeds (17 beef breeds, 14 dairy breeds, and 6 dual-purpose breeds), and four different beef-vs.-dairy comparisons, which we subsequently grouped into core selective sweep (CSS) regions, defined as consecutive signals within 1 Mb of each other. We defined a total of 409 CSSs across the 29 bovine autosomes, 232 (57%) of which were associated with a single-breed (Single-breed CSSs), 134 CSSs (33%) were associated with a limited number of breeds (Two-to-Four-breed CSSs) and 39 CSSs (9%) were associated with five or more breeds (Multi-breed CSSs). For each CSS, we performed a candidate gene survey that identified 291 genes within the CSS intervals (from the total list of 5183 BioMart-extracted genes) linked to dairy and meat production, stature, and coat color traits. A complementary functional enrichment analysis of the CSS positional candidates highlighted other genes related to pathways underlying behavior, immune response, and reproductive traits. The Single-breed CSSs revealed an over-representation of genes related to dairy and beef production, this was further supported by over-representation of production-related pathway terms in these regions based on a functional enrichment analysis. Overall, this review provides a comparative map of the selection sweeps reported in European cattle breeds and presents for the first time a characterization of the selection sweeps that are found in individual breeds. Based on their uniqueness, these breed-specific signals could be considered as “divergence signals,” which may be useful in characterizing and protecting livestock genetic diversity.

Combining GWAS and RNA-Seq Approaches for Detection of the Causal Mutation for Hereditary Junctional Epidermolysis Bullosa in Sheep

In this study, we demonstrate the use of a genome-wide association mapping together with RNA-seq in a reduced number of samples, as an efficient approach to detect the causal mutation for a Mendelian disease. Junctional epidermolysis bullosa is a recessive genodermatosis that manifests with neonatal mechanical fragility of the skin, blistering confined to the lamina lucida of the basement membrane and severe alteration of the hemidesmosomal junctions. In Spanish Churra sheep, junctional epidermolysis bullosa (JEB) has been detected in two commercial flocks. The JEB locus was mapped to Ovis aries chromosome 11 by GWAS and subsequently fine-mapped to an 868-kb homozygous segment using the identical-by-descent method. The ITGB4, which is located within this region, was identified as the best positional and functional candidate gene. The RNA-seq variant analysis enabled us to discover a 4-bp deletion within exon 33 of the ITGB4 gene (c.4412_4415del). The c.4412_4415del mutation causes a frameshift resulting in a premature stop codon at position 1472 of the integrin β4 protein. A functional analysis of this deletion revealed decreased levels of mRNA in JEB skin samples and the absence of integrin β4 labeling in immunohistochemical assays. Genotyping of c.4412_4415del showed perfect concordance with the recessive mode of the disease phenotype. Selection against this causal mutation will now be used to solve the problem of JEB in flocks of Churra sheep. Furthermore, the identification of the ITGB4 mutation means that affected sheep can be used as a large mammal animal model for the human form of epidermolysis bullosa with aplasia cutis. Our approach evidences that RNA-seq offers cost-effective alternative to identify variants in the species in which high resolution exome-sequencing is not straightforward.

Application of selection mapping to identify genomic regions associated with dairy production in sheep

In Europe, especially in Mediterranean areas, the sheep has been traditionally exploited as a dual purpose species, with income from both meat and milk. Modernization of husbandry methods and the establishment of breeding schemes focused on milk production have led to the development of "dairy breeds." This study investigated selective sweeps specifically related to dairy production in sheep by searching for regions commonly identified in different European dairy breeds. With this aim, genotypes from 44,545 SNP markers covering the sheep autosomes were analysed in both European dairy and non-dairy sheep breeds using two approaches: (i) identification of genomic regions showing extreme genetic differentiation between each dairy breed and a closely related non-dairy breed, and (ii) identification of regions with reduced variation (heterozygosity) in the dairy breeds using two methods. Regions detected in at least two breeds (breed pairs) by the two approaches (genetic differentiation and at least one of the heterozygosity-based analyses) were labeled as core candidate convergence regions and further investigated for candidate genes. Following this approach six regions were detected. For some of them, strong candidate genes have been proposed (e.g. ABCG2, SPP1), whereas some other genes designated as candidates based on their association with sheep and cattle dairy traits (e.g. LALBA, DGAT1A) were not associated with a detectable sweep signal. Few of the identified regions were coincident with QTL previously reported in sheep, although many of them corresponded to orthologous regions in cattle where QTL for dairy traits have been identified. Due to the limited number of QTL studies reported in sheep compared with cattle, the results illustrate the potential value of selection mapping to identify genomic regions associated with dairy traits in sheep.

Identification of a 31-bp deletion in the RELN gene causing lissencephaly with cerebellar hypoplasia in sheep

Lissencephaly is an inherited developmental disorder in which neuronal migration is impaired. A type of lissencephaly associated with cerebellar hypoplasia (LCH) was diagnosed in a commercial flock of Spanish Churra sheep. The genotyping of 7 affected animals and 33 controls with the OvineSNP50 BeadChip enabled the localization of the causative mutation for ovine LCH to a 4.8-Mb interval on sheep chromosome 4 using genome-wide association and homozygosity mapping. The RELN gene, which is located within this interval, was considered a strong positional and functional candidate because it plays critical roles in neuronal migration and layer formation. By performing a sequencing analysis of this gene’s specific mRNA in a control lamb, we obtained the complete CDS of the ovine RELN gene. The cDNA sequence from an LCH-affected lamb revealed a deletion of 31 bp (c.5410_5440del) in predicted exon 36 of RELN, resulting in a premature termination codon. A functional analysis of this mutation revealed decreased levels of RELN mRNA and a lack of reelin protein in the brain cortex and blood of affected lambs. This mutation showed a complete concordance with the Mendelian recessive pattern of inheritance observed for the disease. The identification of the causal mutation of LCH in Churra sheep will facilitate the implementation of gene-assisted selection to detect heterozygous mutants, which will help breeders avoid at-risk matings in their flocks. Moreover, the identification of this naturally occurring RELN mutation provides an opportunity to use Churra sheep as a genetically characterized large animal model for the study of reelin functions in the developing and mature brain.

Identification of quantitative trait loci underlying milk traits in Spanish dairy sheep using linkage plus combined linkage disequilibrium and linkage analysis approaches

In this study, 2 procedures were used to analyze a data set from a whole-genome scan, one based on linkage analysis information and the other combing linkage disequilibrium and linkage analysis (LDLA), to determine the quantitative trait loci (QTL) influencing milk production traits in sheep. A total of 1,696 animals from 16 half-sib families were genotyped using the OvineSNP50 BeadChip (Illumina Inc., San Diego, CA) and analysis was performed using a daughter design. Moreover, the same data set has been previously investigated through a genome-wide association (GWA) analysis and a comparison of results from the 3 methods has been possible. The linkage analysis and LDLA methodologies yielded different results, although some significantly associated regions were common to both procedures. The linkage analysis detected 3 overlapping genome-wise significant QTL on sheep chromosome (OAR) 2 influencing milk yield, protein yield, and fat yield, whereas 34 genome-wise significant QTL regions were detected using the LDLA approach. The most significant QTL for protein and fat percentages was detected on OAR3, which was reported in a previous GWA analysis. Both the linkage analysis and LDLA identified many other chromosome-wise significant associations across different sheep autosomes. Additional analyses were performed on OAR2 and OAR3 to determine the possible causality of the most significant polymorphisms identified for these genetic effects by the previously reported GWA analysis. For OAR3, the analyses demonstrated additional genetic proof of the causality previously suggested by our group for a single nucleotide polymorphism located in the α-lactalbumin gene (LALBA). In summary, although the results shown here suggest that in commercial dairy populations, the LDLA method exhibits a higher efficiency to map QTL than the simple linkage analysis or linkage disequilibrium methods, we believe that comparing the 3 analysis methods is the best approach to obtain a global picture of all identifiable QTL segregating in the population at both family-based and population-based levels.

GWA analysis for milk production traits in dairy sheep and genetic support for a QTN influencing milk protein percentage in the LALBA gene

In this study, we used the Illumina OvineSNP50 BeadChip to conduct a genome-wide association (GWA) analysis for milk production traits in dairy sheep by analyzing a commercial population of Spanish Churra sheep. The studied population consisted of a total of 1,681 Churra ewes belonging to 16 half-sib families with available records for milk yield (MY), milk protein and fat yields (PY and FY) and milk protein and fat contents (PP and FP). The most significant association identified reached experiment-wise significance for PP and FP and was located on chromosome 3 (OAR3). These results confirm the population-level segregation of a previously reported QTL affecting PP and suggest that this QTL has a significant pleiotropic effect on FP. Further associations were detected at the chromosome-wise significance level on 14 other chromosomal regions. The marker on OAR3 showing the highest significant association was located at the third intron of the alpha-lactalbumin (LALBA) gene, which is a functional and positional candidate underlying this association. Sequencing this gene in the 16 Churra rams of the studied resource population identified additional polymorphisms. One out of the 31 polymorphisms identified was located within the coding gene sequence (LALBA_g.242T>C) and was predicted to cause an amino acid change in the protein (Val27Ala). Different approaches, including GWA analysis, a combined linkage and linkage disequilibrium study and a concordance test with the QTL segregating status of the sires, were utilized to assess the role of this mutation as a putative QTN for the genetic effects detected on OAR3. Our results strongly support the polymorphism LALBA_g.242T>C as the most likely causal mutation of the studied OAR3 QTL affecting PP and FP, although we cannot rule out the possibility that this SNP is in perfect linkage disequilibrium with the true causal polymorphism.

Analysis of genetic diversity in Bolivian llama populations using microsatellites

South American camelids (SACs) have a major role in the maintenance and potential future of rural Andean human populations. More than 60% of the 3.7 million llamas living worldwide are found in Bolivia. Due to the lack of studies focusing on genetic diversity in Bolivian llamas, this analysis investigates both the genetic diversity and structure of 12 regional groups of llamas that span the greater part of the range of distribution for this species in Bolivia. The analysis of 42 microsatellite markers in the considered regional groups showed that, in general, there were high levels of polymorphism (a total of 506 detected alleles; average PIC across per marker: 0.66), which are comparable with those reported for other populations of domestic SACs. The estimated diversity parameters indicated that there was high intrapopulational genetic variation (average number of alleles and average expected heterozygosity per marker: 12.04 and 0.68, respectively) and weak genetic differentiation among populations (FST range: 0.003-0.052). In agreement with these estimates, Bolivian llamas showed a weak genetic structure and an intense gene flow between all the studied regional groups, which is due to the exchange of reproductive males between the different flocks. Interestingly, the groups for which the largest pairwise FST estimates were observed, Sud Lípez and Nor Lípez, showed a certain level of genetic differentiation that is probably due to the pattern of geographic isolation and limited communication infrastructures of these southern localities. Overall, the population parameters reported here may serve as a reference when establishing conservation policies that address Bolivian llama populations.

Linkage disequilibrium and inbreeding estimation in Spanish Churra sheep

Background

Genomic technologies, such as high-throughput genotyping based on SNP arrays, have great potential to decipher the genetic architecture of complex traits and provide background information concerning genome structure in domestic animals, including the extent of linkage disequilibrium (LD) and haplotype blocks. The objective of this study was to estimate LD, the population evolution (past effective population size) and the level of inbreeding in Spanish Churra sheep.

Results

A total of 43,784 SNPs distributed in the ovine autosomal genome was analyzed in 1,681 Churra ewes. LD was assessed by measuring r² between all pairs of loci. For SNPs up to 10 kb apart, the average r² was 0.329; for SNPs separated by 200–500 kb the average r² was 0.061. When SNPs are separated by more than 50 Mbp, the average r² is the same as between non-syntenic SNP pairs (0.003). The effective population size has decreased through time, faster from 1,000 to 100 years ago and slower since the selection scheme started (15–25 generations ago). In the last generation, four years ago, the effective population size was estimated to be 128 animals. Inbreeding coefficients, although differed depending on the estimation approaches, were generally low and showed the same trend, which indicates that since 2003, inbreeding has been slightly increasing in the studied resource population.

Conclusions

The extent of LD in Churra sheep persists over much more limited distances than reported in dairy cattle and seems to be similar to other ovine populations. Churra sheep show a wide genetic base, with a long-term viable effective population size that has been slightly decreasing since selection scheme began in 1986. The genomic dataset analyzed provided useful information for identifying low-level inbreeding in the sample, whereas based on the parameters reported here, a higher marker density than that analyzed here will be needed to successfully conduct accurate mapping of genes underlying production traits and genomic selection prediction in this sheep breed. Although the Ovine Assembly development is still in a draft stage and future refinements will provide a more accurate physical map that will improve LD estimations, this work is a first step towards the understanding of the genetic architecture in sheep.

Analysis of mitochondrial DNA in Bolivian llama, alpaca and vicuna populations: a contribution to the phylogeny of the South American camelids

The objectives of this work were to assess the mtDNA diversity of Bolivian South American camelid (SAC) populations and to shed light on the evolutionary relationships between the Bolivian camelids and other populations of SACs. We have analysed two different mtDNA regions: the complete coding region of the MT-CYB gene and 513 bp of the D-loop region. The populations sampled included Bolivian llamas, alpacas and vicunas, and Chilean guanacos. High levels of genetic diversity were observed in the studied populations. In general, MT-CYB was more variable than D-loop. On a species level, the vicunas showed the lowest genetic variability, followed by the guanacos, alpacas and llamas. Phylogenetic analyses performed by including additional available mtDNA sequences from the studied species confirmed the existence of the two monophyletic clades previously described by other authors for guanacos (G) and vicunas (V). Significant levels of mtDNA hybridization were found in the domestic species. Our sequence analyses revealed significant sequence divergence within clade G, and some of the Bolivian llamas grouped with the majority of the southern guanacos. This finding supports the existence of more than the one llama domestication centre in South America previously suggested on the basis of archaeozoological evidence. Additionally, analysis of D-loop sequences revealed two new matrilineal lineages that are distinct from the previously reported G and V clades. The results presented here represent the first report on the population structure and genetic variability of Bolivian camelids and may help to elucidate the complex and dynamic domestication process of SAC populations.

Investigation of the genetic architecture of a bone carcass weight QTL on BTA6

A previous analysis of an F(2) /Backcross Charolais × Holstein cross population identified the presence of a highly significant QTL on chromosome 6 (BTA6) affecting the proportion of bone in the carcass. Two closely linked QTL affected birth weight (BW) and body length at birth (BBL). In this report, the marker density around the QTL on BTA6 was increased, adding four additional microsatellite markers across the chromosome and 46 SNPs within the target QTL confidence interval. Of the SNPs, 26 were in positional candidate genes and the remaining 20 provided an even distribution of markers in the target QTL region. As a bone-related trait, the sum of the bone weight for all the left fore- and hindquarter joints of the carcass was analysed. We also studied the BW and BBL. Analyses of the data substantially reduced the QTL confidence interval. No strong evidence was found that the QTL for the three traits studied are different, and we conclude that the results are consistent with a single pleiotropic QTL influencing the three traits, with the largest effects on the proportion of bone in the carcass. The analyses also suggest that none of the SNPs tested is the sole causative variant of the QTL effects. Specifically, the SNP in the NCAPG gene previously reported as a causal mutation for foetal growth and carcass traits in other cattle populations was excluded as the causal mutation for the QTL reported here. Polymorphisms located in other previously identified candidate genes including SPP1, ABCG2, IBSP, MEPE and PPARGC1A were also excluded. The results suggest that SNP51_BTA-119876 is the polymorphism in strongest linkage disequilibrium with the causal mutation(s). Further research is required to identify the causal variant(s) associated with this bone-related QTL.

Replication and refinement of a quantitative trait locus influencing milk protein percentage on ovine chromosome 3

A previous genome scan that was conducted in Spanish Churra sheep identified a significant quantitative trait locus (QTL) for milk protein percentage (PP) on chromosome 3 (OAR3), between markers KD103 and OARVH34. The aim of this study was to replicate these results and to refine the mapped position of this QTL. To accomplish this goal, we analysed 14 new half-sib families of Spanish Churra sheep including 1661 ewes from 29 different flocks. These animals were genotyped for 21 microsatellite markers mapping to OAR3. In addition to a classical linkage analysis (LA), a combined linkage disequilibrium and linkage analysis (LDLA) was performed with the aim of enhancing the resolution of the QTL mapping. The LA that was performed in this sheep population identified the presence of a highly significant QTL for PP near marker KD103 (P(c) < 0.001; P(exp) < 0.001). The phenotypic variance that was owing to the QTL was 2.74%. Two segregating families for the target QTL were identified in this population with QTL effect estimates of 0.47 and 0.95 SD. The LDLA identified the same QTL as the previous analyses with a high level of statistical significance (P = 9.184 E-11) and narrowed the confidence interval (CI) to a 13 cM region. These results confirm the segregation of the previously identified OAR3 QTL that influences PP in Spanish Churra sheep. Future research will aim to increase the marker density across the refined CI and to analyse the corresponding candidate genes to identify the allelic variant or variants that underlie this genetic effect.

The role of bovine causal genes underlying dairy traits in Spanish Churra sheep

In dairy cattle, quantitative trait nucleotides (QTNs) underlying quantitative trait loci (QTL) for milk production traits have been identified in bovine DGAT1, GHR and ABCG2 genes. The SPP1 gene has also been proposed to be a regulator of lactation. In sheep, QTL underlying milk production traits have been reported only recently, and no proven QTN has been identified. Taking into account the close phylogenetic relationship between sheep and cattle, this study examined the possible effects of the aforementioned genes on sheep milk production traits. We first studied the genetic variability of the DGAT1, GHR, ABCG2 and SPP1 genes in 15 rams of the Spanish Churra dairy sheep breed. Second, we performed an association analysis between SNPs identified in these genes and three milk production traits recorded in a commercial population of Churra sheep. This analysis revealed only three significant associations at the nominal level (P-value <0.05) involving allelic variants of the ABCG2 gene, whereas no significant association was found for the DGAT1, GHR and SPP1 genes. When the Bonferroni correction was applied to take into account the multiple tests performed, none of the associations identified at the nominal level remained significant. Nevertheless, taking into account the high level of false-negative findings that can arise when applying the stringent Bonferroni correction, we think that our results provide a valuable primary assessment of strong candidate genes for milk traits in sheep.

Information content in genome-wide scans: concordance between patterns of genetic differentiation and linkage mapping associations

BACKGROUND

Scanning the genome with high density SNP markers has become a standard approach for identifying regions of the genome showing substantial between-population genetic differentiation, and thus evidence of diversifying selection. Such regions may contain genes of large phenotypic effect. However, few studies have attempted to address the power or efficacy of such an approach.

RESULTS

In this study, the patterns of allele frequency differences between two cattle breeds based on the Bovine HapMap study were compared with statistical evidence for QTL based on a linkage mapping study of an experimental population formed by a cross between the same breeds. Concordance between the two datasets was seen for chromosomes carrying QTL with strong statistical support, such as BTA5 and BTA18, which carry genes associated with coat color. For these chromosomes, there was a correspondence between the strength of the QTL signal along the chromosome and the degree of genetic differentiation between breeds. However, such an association was not seen in a broader comparison that also included chromosomes carrying QTL with lower significance levels. In addition, other chromosomal regions with substantial QTL effects did not include markers showing extreme between-breed genetic differentiation. Furthermore, the overall consistency between the two studies was weak, with low genome-wide correlation between the statistical values obtained in the linkage mapping study and between-breed genetic differentiation from the HapMap study.

CONCLUSIONS

These results suggest that genomic diversity scans are capable of detecting regions associated with qualitative traits but may be limited in their power to detect regions associated with quantitative phenotypic differences between populations, which may depend on the marker resolution of the study and the level of LD in the populations under investigation.

Detection of quantitative trait loci affecting the milk fatty acid profile on sheep chromosome 22: role of the stearoyl-CoA desaturase gene in Spanish Churra sheep

Sheep milk fat contains several components that may provide human health benefits, such as monounsaturated fatty acids and conjugated linoleic acid (CLA). Most of the CLA in ruminant milk is synthesized in the mammary gland by the action of the enzyme stearoyl-CoA desaturase (SCD) on circulating vaccenic acid (trans-11 C18:2; VA). Previous studies have found significant associations between polymorphisms in the SCD gene and the fatty acid composition of ruminant products, including sheep milk. Based on this, we performed a quantitative trait loci (QTL) analysis of an ovine chromosome (22) that harbors the SCD gene for effects on milk fatty acid composition traits and classical milk production traits. We identified a suggestive QTL influencing the CLA/VA ratio with the maximum statistic at position 26 cM of the studied chromosome, whereas the SCD gene has been mapped to position 41.6 cM. The individual introduction of 4 SCD single nucleotide polymorphisms in the QTL model did not cause a reduction of the variance explained by the QTL, which suggests that the SCD gene is not directly responsible for the detected effect in the Churra population studied herein. This conclusion was supported by the lack of any significant association identified between the 4 SCD single nucleotide polymorphisms and the CLA/VA ratio. This association analysis suggested a possible effect of the SCD gene on milk fat percentage in Churra sheep. An independent confirmation of these primary results will be required before attempting its practical implementation in selection programs.

The identification of QTL that affect the fatty acid composition of milk on sheep chromosome 11

In this work, we analysed 11 genetic markers localized on OAR11 in a commercial population of Spanish Churra sheep to detect QTL that underlie milk fatty acid (FA) composition traits. Following a daughter design, we analysed 799 ewes distributed in 15 half-sib families. Eight microsatellite markers and three novel SNPs identified in two genes related to fatty acid metabolism, acetyl-CoA carboxylase alpha (ACACA) and fatty acid synthase (FASN), were genotyped in the whole population under study. The phenotypic traits considered in the study included 22 measurements related to the FA composition of the milk and three other milk production traits (milk protein percentage, milk fat percentage and milk yield). Across-family regression analysis revealed four significant QTL at the 5% chromosome-wise level influencing contents of capric acid (C10:0), lauric acid (C12:0), linoleic conjugated acid (CLA) and polyunsaturated fatty acids (PUFA) respectively. The peaks of the QTL affecting C10:0 and PUFA contents in milk map close to the FASN gene, which has been evaluated as a putative positional candidate for these QTL. The QTL influencing C12:0 content reaches its maximum significance at 58 cM, close to the gene coding for the glucose-dependent insulinotropic polypeptide. We were not able to find any candidate genes related to fat metabolism at the QTL influencing CLA content, which is located at the proximal end of the chromosome. Further research efforts will be needed to confirm and refine the QTL locations reported here.