A 12 kb multi-allelic copy number variation encompassing a GC gene enhancer is associated with mastitis resistance in dairy cattle

Copy number variant scan in more than four thousand Holstein cows bred in Lombardy, Italy

Copy Number Variants (CNV) are modifications affecting the genome sequence of DNA, for instance, they can be duplications or deletions of a considerable number of base pairs (i.e., greater than 1000 bp and up to millions of bp). Their impact on the variation of the phenotypic traits has been widely demonstrated. In addition, CNVs are a class of markers useful to identify the genetic biodiversity among populations related to adaptation to the environment. The aim of this study was to detect CNVs in more than four thousand Holstein cows, using information derived by a genotyping done with the GGP (GeneSeek Genomic Profiler) bovine 100K SNP chip. To detect CNV the SVS 8.9 software was used, then CNV regions (CNVRs) were detected. A total of 123,814 CNVs (4,150 non redundant) were called and aggregated into 1,397 CNVRs. The PCA results obtained using the CNVs information, showed that there is some variability among animals. For many genes annotated within the CNVRs, the role in immune response is well known, as well as their association with important and economic traits object of selection in Holstein, such as milk production and quality, udder conformation and body morphology. Comparison with reference revealed unique CNVRs of the Holstein breed, and others in common with Jersey and Brown. The information regarding CNVs represents a valuable resource to understand how this class of markers may improve the accuracy in prediction of genomic value, nowadays solely based on SNPs markers.

Unraveling the genomic background of resilience based on variability in milk yield and milk production levels in North American Holstein cattle through genome-wide association study and Mendelian randomization analyses

Breeding more resilient animals will benefit the dairy cattle industry in the long term, especially as global climate changes become more severe. Previous studies have reported genetic parameters for various milk yield-based resilience indicators, but the underlying genomic background of these traits remain unknown. In this study, we conducted GWAS of 62,029 SNPs with 4 milk yield-based resilience indicators, including the weighted occurrence frequency (wfPert) and accumulated milk losses (dPert) of milk yield perturbations, and log-transformed variance (LnVar) and lag-1 autocorrelation (rauto) of daily yield residuals. These variables were previously derived from 5.6 million daily milk yield records from 21,350 lactations (parities 1-3) of 11,787 North American Holstein cows. The average daily milk yield (ADMY) throughout lactation was also included to compare the shared genetic background of resilience indicators with milk yield. The differential genetic background of these indicators was first revealed by the significant genomic regions identified and significantly enriched biological pathways of positional candidate genes, which confirmed the genetic difference among resilience indicators. Interestingly, the functional analyses of candidate genes suggested that the regulation of intestinal homeostasis is most likely affecting resilience derived based on variability in milk yield. Based on Mendelian randomization analyses of multiple instrumental SNPs, we further found an unfavorable causal association of ADMY with LnVar. In conclusion, the resilience indicators evaluated are genetically different traits, and there are causal associations of milk yield with some of the resilience indicators evaluated. In addition to providing biological insights into the molecular regulation mechanisms of resilience derived based on variability in milk yield, this study also indicates the need for developing selection indexes combining multiple indicator traits and taking into account their genetic relationship for breeding more resilient dairy cattle.

Identification of genetic variants and individual genes associated with postpartum hypocalcemia in Holstein cows

Periparturient hypocalcemia is a complex metabolic disorder that occurs at the onset of lactation because of a sudden irreversible loss of Ca incorporated into colostrum and milk. Some cows are unable to quickly adapt to this demand and succumb to clinical hypocalcemia, commonly known as milk fever, whereas a larger proportion of cows develop subclinical hypocalcemia. The main goal of this study was to identify causative mutations and candidate genes affecting postpartum blood calcium concentration in Holstein cows. Data consisted of blood calcium concentration measured in 2513 Holstein cows on the first three days after parturition. All cows had genotypic information for 79 k SNP markers. Two consecutive rounds of imputation were performed: first, the 2513 Holstein cows were imputed from 79 k to 312 k SNP markers. This imputation was performed using a reference set of 17,131 proven Holstein bulls with 312 k SNP markers. Then, the 2513 Holstein cows were imputed from 312 k markers to whole-genome sequence data. This second round of imputation used 179 Holstein animals from the 1000 Bulls Genome Project as a reference set. Three alternative phenotypes were evaluated: (1) total calcium concentration in the first 24 h postpartum, (2) total calcium concentration in the first 72 h postpartum calculated as the area under the curve; and (3) the recovery of total calcium concentration calculated as the difference in total calcium concentration between 72 and 24 h. The identification of genetic variants associated with these traits was performed using a two-step mixed model-based approach implemented in the R package MixABEL. The most significant variants were located within or near genes involved in calcium homeostasis and vitamin D transport (GC), calcium and potassium channels (JPH3 and KCNK13), energy and lipid metabolism (CA5A, PRORP, and SREBP1), and immune response (IL12RB2 and CXCL8), among other functions. This work provides the foundation for the development of novel breeding and management tools for reducing the incidence of periparturient hypocalcemia in dairy cattle.

Structural variants and short tandem repeats impact gene expression and splicing in bovine testis tissue

Structural variants (SVs) and short tandem repeats (STRs) are significant sources of genetic variation. However, the impacts of these variants on gene regulation have not been investigated in cattle. Here, we genotyped and characterized 19,408 SVs and 374,821 STRs in 183 bovine genomes and investigated their impact on molecular phenotypes derived from testis transcriptomes. We found that 71% STRs were multiallelic. The vast majority (95%) of STRs and SVs were in intergenic and intronic regions. Only 37% SVs and 40% STRs were in high linkage disequilibrium (LD) (R2 > 0.8) with surrounding SNPs/insertions and deletions (Indels), indicating that SNP-based association testing and genomic prediction are blind to a nonnegligible portion of genetic variation. We showed that both SVs and STRs were more than 2-fold enriched among expression and splicing QTL (e/sQTL) relative to SNPs/Indels and were often associated with differential expression and splicing of multiple genes. Deletions and duplications had larger impacts on splicing and expression than any other type of SV. Exonic duplications predominantly increased gene expression either through alternative splicing or other mechanisms, whereas expression- and splicing-associated STRs primarily resided in intronic regions and exhibited bimodal effects on the molecular phenotypes investigated. Most e/sQTL resided within 100 kb of the affected genes or splicing junctions. We pinpoint candidate causal STRs and SVs associated with the expression of SLC13A4 and TTC7B and alternative splicing of a lncRNA and CAPP1. We provide a catalog of STRs and SVs for taurine cattle and show that these variants contribute substantially to gene expression and splicing variation.

An organism-wide ATAC-seq peak catalog for the bovine and its use to identify regulatory variants

We report the generation of an organism-wide catalog of 976,813 cis-acting regulatory elements for the bovine detected by the assay for transposase accessible chromatin using sequencing (ATAC-seq). We regroup these regulatory elements in 16 components by nonnegative matrix factorization. Correlation between the genome-wide density of peaks and transcription start sites, correlation between peak accessibility and expression of neighboring genes, and enrichment in transcription factor binding motifs support their regulatory potential. Using a previously established catalog of 12,736,643 variants, we show that the proportion of single-nucleotide polymorphisms mapping to ATAC-seq peaks is higher than expected and that this is owing to an approximately 1.3-fold higher mutation rate within peaks. Their site frequency spectrum indicates that variants in ATAC-seq peaks are subject to purifying selection. We generate eQTL data sets for liver and blood and show that variants that drive eQTL fall into liver- and blood-specific ATAC-seq peaks more often than expected by chance. We combine ATAC-seq and eQTL data to estimate that the proportion of regulatory variants mapping to ATAC-seq peaks is approximately one in three and that the proportion of variants mapping to ATAC-seq peaks that are regulatory is approximately one in 25. We discuss the implication of these findings on the utility of ATAC-seq information to improve the accuracy of genomic selection.

Robust identification of regulatory variants (eQTLs) using a differential expression framework developed for RNA-sequencing

BACKGROUND

A gap currently exists between genetic variants and the underlying cell and tissue biology of a trait, and expression quantitative trait loci (eQTL) studies provide important information to help close that gap. However, two concerns that arise with eQTL analyses using RNA-sequencing data are normalization of data across samples and the data not following a normal distribution. Multiple pipelines have been suggested to address this. For instance, the most recent analysis of the human and farm Genotype-Tissue Expression (GTEx) project proposes using trimmed means of M-values (TMM) to normalize the data followed by an inverse normal transformation.

RESULTS

In this study, we reasoned that eQTL analysis could be carried out using the same framework used for differential gene expression (DGE), which uses a negative binomial model, a statistical test feasible for count data. Using the GTEx framework, we identified 35 significant eQTLs (P < 5 × 10^-8) following the ANOVA model and 39 significant eQTLs (P < 5 × 10^-8) following the additive model. Using a differential gene expression framework, we identified 930 and six significant eQTLs (P < 5 × 10^-8) following an analytical framework equivalent to the ANOVA and additive model, respectively. When we compared the two approaches, there was no overlap of significant eQTLs between the two frameworks. Because we defined specific contrasts, we identified trans eQTLs that more closely resembled what we expect from genetic variants showing complete dominance between alleles. Yet, these were not identified by the GTEx framework.

CONCLUSIONS

Our results show that transforming RNA-sequencing data to fit a normal distribution prior to eQTL analysis is not required when the DGE framework is employed. Our proposed approach detected biologically relevant variants that otherwise would not have been identified due to data transformation to fit a normal distribution.

High-resolution structural variants catalogue in a large-scale whole genome sequenced bovine family cohort data

BACKGROUND

Structural variants (SVs) are chromosomal segments that differ between genomes, such as deletions, duplications, insertions, inversions and translocations. The genomics revolution enabled the discovery of sub-microscopic SVs via array and whole-genome sequencing (WGS) data, paving the way to unravel the functional impact of SVs. Recent human expression QTL mapping studies demonstrated that SVs play a disproportionally large role in altering gene expression, underlining the importance of including SVs in genetic analyses. Therefore, this study aimed to generate and explore a high-quality bovine SV catalogue exploiting a unique cattle family cohort data (total 266 samples, forming 127 trios).

RESULTS

We curated 13,731 SVs segregating in the population, consisting of 12,201 deletions, 1,509 duplications, and 21 multi-allelic CNVs (> 50-bp). Of these, we validated a subset of copy number variants (CNVs) utilising a direct genotyping approach in an independent cohort, indicating that at least 62% of the CNVs are true variants, segregating in the population. Among gene-disrupting SVs, we prioritised two likely high impact duplications, encompassing ORM1 and POPDC3 genes, respectively. Liver expression QTL mapping results revealed that these duplications are likely causing altered gene expression, confirming the functional importance of SVs. Although most of the accurately genotyped CNVs are tagged by single nucleotide polymorphisms (SNPs) ascertained in WGS data, most CNVs were not captured by individual SNPs obtained from a 50K genotyping array.

CONCLUSION

We generated a high-quality SV catalogue exploiting unique whole genome sequenced bovine family cohort data. Two high impact duplications upregulating the ORM1 and POPDC3 are putative candidates for postpartum feed intake and hoof health traits, thus warranting further investigation. Generally, CNVs were in low LD with SNPs on the 50K array. Hence, it remains crucial to incorporate CNVs via means other than tagging SNPs, such as investigation of tagging haplotypes, direct imputation of CNVs, or direct genotyping as done in the current study. The SV catalogue and the custom genotyping array generated in the current study will serve as valuable resources accelerating utilisation of full spectrum of genetic variants in bovine genomes.

Resequencing of a Pekin duck breeding population provides insights into the genomic response to short-term artificial selection

BACKGROUND

Short-term, intense artificial selection drives fast phenotypic changes in domestic animals and leaves imprints on their genomes. However, the genetic basis of this selection response is poorly understood. To better address this, we employed the Pekin duck Z2 pure line, in which the breast muscle weight was increased nearly 3-fold after 10 generations of breeding. We denovo assembled a high-quality reference genome of a female Pekin duck of this line (GCA_003850225.1) and identified 8.60 million genetic variants in 119 individuals among 10 generations of the breeding population.

RESULTS

We identified 53 selected regions between the first and tenth generations, and 93.8% of the identified variations were enriched in regulatory and noncoding regions. Integrating the selection signatures and genome-wide association approach, we found that 2 regions covering 0.36 Mb containing UTP25 and FBRSL1 were most likely to contribute to breast muscle weight improvement. The major allele frequencies of these 2 loci increased gradually with each generation following the same trend. Additionally, we found that a copy number variation region containing the entire EXOC4 gene could explain 1.9% of the variance in breast muscle weight, indicating that the nervous system may play a role in economic trait improvement.

CONCLUSIONS

Our study not only provides insights into genomic dynamics under intense artificial selection but also provides resources for genomics-enabled improvements in duck breeding.

Single-step genome-wide association analyses of claw horn lesions in Holstein cattle using linear and threshold models

BACKGROUND

Lameness in dairy cattle is primarily caused by foot lesions including the claw horn lesions (CHL) of sole haemorrhage (SH), sole ulcers (SU), and white line disease (WL). This study investigated the genetic architecture of the three CHL based on detailed animal phenotypes of CHL susceptibility and severity. Estimation of genetic parameters and breeding values, single-step genome-wide association analyses, and functional enrichment analyses were performed.

RESULTS

The studied traits were under genetic control with a low to moderate heritability. Heritability estimates of SH and SU susceptibility on the liability scale were 0.29 and 0.35, respectively. Heritability of SH and SU severity were 0.12 and 0.07, respectively. Heritability of WL was relatively lower, indicating stronger environmental influence on the presence and development of WL than the other two CHL. Genetic correlations between SH and SU were high (0.98 for lesion susceptibility and 0.59 for lesion severity), whereas genetic correlations of SH and SU with WL also tended to be positive. Candidate quantitative trait loci (QTL) were identified for all CHL, including some on Bos taurus chromosome (BTA) 3 and 18 with potential pleiotropic effects associated with multiple foot lesion traits. A genomic window of 0.65 Mb on BTA3 explained 0.41, 0.50, 0.38, and 0.49% of the genetic variance for SH susceptibility, SH severity, WL susceptibility, and WL severity, respectively. Another window on BTA18 explained 0.66, 0.41, and 0.70% of the genetic variance for SH susceptibility, SU susceptibility, and SU severity, respectively. The candidate genomic regions associated with CHL harbour annotated genes that are linked to immune system function and inflammation responses, lipid metabolism, calcium ion activities, and neuronal excitability.

CONCLUSIONS

The studied CHL are complex traits with a polygenic mode of inheritance. Most traits exhibited genetic variation suggesting that animal resistance to CHL can be improved with breeding. The CHL traits were positively correlated, which will facilitate genetic improvement for resistance to CHL as a whole. Candidate genomic regions associated with lesion susceptibility and severity of SH, SU, and WL provide insights into a global profile of the genetic background underlying CHL and inform genetic improvement programmes aiming at enhancing foot health in dairy cattle.

An effect of large-scale deletions and duplications on transcript expression

Since copy number variants (CNVs) have been recognized as an important source of genetic and transcriptomic variation, we aimed to characterize the impact of CNVs located within coding, intergenic, upstream, and downstream gene regions on the expression of transcripts. Regions in which deletions occurred most often were introns, while duplications in coding regions. The transcript expression was lower for deleted coding (P = 0.008) and intronic regions (P = 1.355 × 10^-10), but it was not changed in the case of upstream and downstream gene regions (P = 0.085). Moreover, the expression was decreased if duplication occurred in the coding region (P = 8.318 × 10^-5). Furthermore, a negative correlation (r =  - 0.27) between transcript length and its expression was observed. The correlation between the percent of deleted/duplicated transcript and transcript expression level was not significant for all concerned genomic regions in five out of six animals. The exceptions were deletions in coding regions (P = 0.004) and duplications in introns (P = 0.01) in one individual. CNVs in coding (deletions, duplications) and intronic (deletions) regions are important modulators of transcripts by reducing their expression level. We hypothesize that deletions imply severe consequences by interrupting genes. The negative correlation between the size of the transcript and its expression level found in this study is consistent with the hypothesis that selection favours shorter introns and a moderate number of exons in highly expressed genes. This may explain the transcript expression reduction by duplications. We did not find the correlation between the size of deletions/duplications and transcript expression level suggesting that expression is modulated by CNVs regardless of their size.

Optical mapping compendium of structural variants across global cattle breeds

Structural variants (SV) have been linked to important bovine disease phenotypes, but due to the difficulty of their accurate detection with standard sequencing approaches, their role in shaping important traits across cattle breeds is largely unexplored. Optical mapping is an alternative approach for mapping SVs that has been shown to have higher sensitivity than DNA sequencing approaches. The aim of this project was to use optical mapping to develop a high-quality database of structural variation across cattle breeds from different geographical regions, to enable further study of SVs in cattle. To do this we generated 100X Bionano optical mapping data for 18 cattle of nine different ancestries, three continents and both cattle sub-species. In total we identified 13,457 SVs, of which 1,200 putatively overlap coding regions. This resource provides a high-quality set of optical mapping-based SV calls that can be used across studies, from validating DNA sequencing-based SV calls to prioritising candidate functional variants in genetic association studies and expanding our understanding of the role of SVs in cattle evolution.

Measurement(s)	Optical Mapping
Technology Type(s)	Optical Mapping
Factor Type(s)	Structural variants
Sample Characteristic - Organism	Bos taurus
Sample Characteristic - Location	United Kingdom • Kenya • Zambia • Uganda • Brazil • Nigeria

Assembly of a pangenome for global cattle reveals missing sequences and novel structural variations, providing new insights into their diversity and evolutionary history

A cattle pangenome representation was created based on the genome sequences of 898 cattle representing 57 breeds. The pangenome identified 83 Mb of sequence not found in the cattle reference genome, representing 3.1% novel sequence compared with the 2.71-Gb reference. A catalog of structural variants developed from this cattle population identified 3.3 million deletions, 0.12 million inversions, and 0.18 million duplications. Estimates of breed ancestry and hybridization between cattle breeds using insertion/deletions as markers were similar to those produced by single nucleotide polymorphism-based analysis. Hundreds of deletions were observed to have stratification based on subspecies and breed. For example, an insertion of a Bov-tA1 repeat element was identified in the first intron of the APPL2 gene and correlated with cattle breed geographic distribution. This insertion falls within a segment overlapping predicted enhancer and promoter regions of the gene, and could affect important traits such as immune response, olfactory functions, cell proliferation, and glucose metabolism in muscle. The results indicate that pangenomes are a valuable resource for studying diversity and evolutionary history, and help to delineate how domestication, trait-based breeding, and adaptive introgression have shaped the cattle genome.

Structural variant-based pangenome construction has low sensitivity to variability of haplotype-resolved bovine assemblies

Advantages of pangenomes over linear reference assemblies for genome research have recently been established. However, potential effects of sequence platform and assembly approach, or of combining assemblies created by different approaches, on pangenome construction have not been investigated. Here we generate haplotype-resolved assemblies from the offspring of three bovine trios representing increasing levels of heterozygosity that each demonstrate a substantial improvement in contiguity, completeness, and accuracy over the current Bos taurus reference genome. Diploid coverage as low as 20x for HiFi or 60x for ONT is sufficient to produce two haplotype-resolved assemblies meeting standards set by the Vertebrate Genomes Project. Structural variant-based pangenomes created from the haplotype-resolved assemblies demonstrate significant consensus regardless of sequence platform, assembler algorithm, or coverage. Inspecting pangenome topologies identifies 90 thousand structural variants including 931 overlapping with coding sequences; this approach reveals variants affecting QRICH2, PRDM9, HSPA1A, TAS2R46, and GC that have potential to affect phenotype.

Assessment of linkage disequilibrium patterns between structural variants and single nucleotide polymorphisms in three commercial chicken populations

BACKGROUND

Structural variants (SV) are causative for some prominent phenotypic traits of livestock as different comb types in chickens or color patterns in pigs. Their effects on production traits are also increasingly studied. Nevertheless, accurately calling SV remains challenging. It is therefore of interest, whether close-by single nucleotide polymorphisms (SNPs) are in strong linkage disequilibrium (LD) with SVs and can serve as markers. Literature comes to different conclusions on whether SVs are in LD to SNPs on the same level as SNPs to other SNPs. The present study aimed to generate a precise SV callset from whole-genome short-read sequencing (WGS) data for three commercial chicken populations and to evaluate LD patterns between the called SVs and surrounding SNPs. It is thereby the first study that assessed LD between SVs and SNPs in chickens.

RESULTS

The final callset consisted of 12,294,329 bivariate SNPs, 4,301 deletions (DEL), 224 duplications (DUP), 218 inversions (INV) and 117 translocation breakpoints (BND). While average LD between DELs and SNPs was at the same level as between SNPs and SNPs, LD between other SVs and SNPs was strongly reduced (DUP: 40%, INV: 27%, BND: 19% of between-SNP LD). A main factor for the reduced LD was the presence of local minor allele frequency differences, which accounted for 50% of the difference between SNP - SNP and DUP - SNP LD. This was potentially accompanied by lower genotyping accuracies for DUP, INV and BND compared with SNPs and DELs. An evaluation of the presence of tag SNPs (SNP in highest LD to the variant of interest) further revealed DELs to be slightly less tagged by WGS SNPs than WGS SNPs by other SNPs. This difference, however, was no longer present when reducing the pool of potential tag SNPs to SNPs located on four different chicken genotyping arrays.

CONCLUSIONS

The results implied that genomic variance due to DELs in the chicken populations studied can be captured by different SNP marker sets as good as variance from WGS SNPs, whereas separate SV calling might be advisable for DUP, INV, and BND effects.

Benchmarking phasing software with a whole-genome sequenced cattle pedigree

Background

Accurate haplotype reconstruction is required in many applications in quantitative and population genomics. Different phasing methods are available but their accuracy must be evaluated for samples with different properties (population structure, marker density, etc.). We herein took advantage of whole-genome sequence data available for a Holstein cattle pedigree containing 264 individuals, including 98 trios, to evaluate several population-based phasing methods. This data represents a typical example of a livestock population, with low effective population size, high levels of relatedness and long-range linkage disequilibrium.

Results

After stringent filtering of our sequence data, we evaluated several population-based phasing programs including one or more versions of AlphaPhase, ShapeIT, Beagle, Eagle and FImpute. To that end we used 98 individuals having both parents sequenced for validation. Their haplotypes reconstructed based on Mendelian segregation rules were considered the gold standard to assess the performance of population-based methods in two scenarios. In the first one, only these 98 individuals were phased, while in the second one, all the 264 sequenced individuals were phased simultaneously, ignoring the pedigree relationships. We assessed phasing accuracy based on switch error counts (SEC) and rates (SER), lengths of correctly phased haplotypes and the probability that there is no phasing error between a pair of SNPs as a function of their distance. For most evaluated metrics or scenarios, the best software was either ShapeIT4.1 or Beagle5.2, both methods resulting in particularly high phasing accuracies. For instance, ShapeIT4.1 achieved a median SEC of 50 per individual and a mean haplotype block length of 24.1 Mb (scenario 2). These statistics are remarkable since the methods were evaluated with a map of 8,400,000 SNPs, and this corresponds to only one switch error every 40,000 phased informative markers. When more relatives were included in the data (scenario 2), FImpute3.0 reconstructed extremely long segments without errors.

Conclusions

We report extremely high phasing accuracies in a typical livestock sample. ShapeIT4.1 and Beagle5.2 proved to be the most accurate, particularly for phasing long segments and in the first scenario. Nevertheless, most tools achieved high accuracy at short distances and would be suitable for applications requiring only local haplotypes.

Supplementary Information

The online version contains supplementary material available at 10.1186/s12864-022-08354-6.