A multi-tissue atlas of allelic-specific expression reveals the characteristics, mechanisms, and relationship with dominant effects in cattle

A multi-tissue atlas of allelic-specific expression reveals the characteristics, mechanisms, and relationship with dominant effects in cattle | Research Square window.SnipcartSettings = { analytics: { enabled: false } }; (function() { var accessVector = localStorage.getItem('access_vector') || ''; window.dataLayer = window.dataLayer || []; if (accessVector) { window.dataLayer.push({ user: { profile: { profileInfo: { snid: accessVector } } } }); } })(); (function(w,d,s,l,i){w[l]=w[l]||[];w[l].push({'gtm.start':new Date().getTime(),event:'gtm.js'});var f=d.getElementsByTagName(s)[0],j=d.createElement(s),dl=l!='dataLayer'?'&l='+l:'';j.async=true;j.src='https://www.googletagmanager.com/gtm.js?id='+i+dl;f.parentNode.insertBefore(j,f);})(window,document,'script','dataLayer','GTM-K279D39R'); Browse Preprints In Review Journals COVID-19 Preprints AJE Video Bytes Research Tools Research Promotion AJE Professional Editing AJE Rubriq About Preprint Platform In Review Editorial Policies Our Team Advisory Board Help Center Sign In Submit a Preprint Cite Share Download PDF Research Article A multi-tissue atlas of allelic-specific expression reveals the characteristics, mechanisms, and relationship with dominant effects in cattle Jiaqi Li, Lei Xu, Xiaoyun Liang, Letian Li, Xixia Huang, Qiuming Chen This is a preprint; it has not been peer reviewed by a journal. https://doi.org/ 10.21203/rs.3.rs-5530951/v1 This work is licensed under a CC BY 4.0 License Status: Under Review Version 1 posted 11 You are reading this latest preprint version Abstract Background Allele-specific expression (ASE) analysis is a crucial tool for validating expression quantitative trait loci (eQTLs), identifying causal variants associated with complex traits, and investigating the genetic mechanisms underlying heterosis. In this study, we characterized ASE variants across 35 tissues using 7,532 publicly available RNA-seq datasets. Additionally, we explored the mechanisms driving ASE through integration with epigenomic data and examined the relationship between ASE and dominance effects on gene expression and milk-related traits in Holstein cattle. Results ASE variants exhibited stronger tissue specificity and lower reproducibility compared to eQTLs. Interestingly, variants with opposite directional effects demonstrated greater resilience across diverse environments. Functional annotation revealed that ASE variants were predominantly located in enhancer regions during transcription, rather than promoter regions. Furthermore, ASE variants were implicated in post-transcriptional and translational processes, including mutations affecting mRNA splicing and triggering nonsense-mediated decay. Analysis of eQTLs, splicing QTLs (sQTLs), and validated QTLs associated with milk-related traits in Holstein cattle, coupled with enrichment analysis in QTL databases and effect size evaluation, indicated that ASE variants were more closely aligned with dominant effects than additive effects, particularly in reproductive and immune-related tissues/traits, which exhibited higher levels of heterosis. Conclusions Our findings not only enhance our understanding of the genetic mechanisms underlying heterosis and ASE formation but also provide a valuable resource of regulatory variants that can be leveraged to improve economic traits through molecular breeding or the strategic exploitation of heterosis. Figures Figure 1 Figure 2 Figure 3 Figure 4 Figure 5 Introduction Over the past decade, genome-wide association studies (GWAS) have consistently demonstrated that most variants associated with diseases and other complex traits are located in non-coding regions of the genome [ 1 , 2 ]. Understanding how these genetic variants influence gene expression is crucial for uncovering the genetic basis of complex traits. To address this, several Genotype-Tissue Expression (GTEx) projects have been initiated across multiple species [ 3 – 5 ]. These projects aim to identify expression quantitative trait loci (eQTLs), which are genetic variants that affect gene expression levels. Depending on their proximity to the target genes, eQTLs can be classified into cis -eQTLs, which act on nearby genes, and trans -eQTLs, which influence distant genes, sometimes on different chromosomes. In these GTEx initiatives, allelic-specific expression (ASE) has been employed as a robust tool for validating cis -eQTL [ 4 , 6 , 7 ]. ASE measures deviations from the expected equal expression of two alleles at a heterozygous locus within an individual. The internal consistency in environmental and technical factors between alleles within the same individual makes ASE a powerful approach. By integrating ASE with other approaches, such as eQTL mapping [ 8 ], GWAS [ 9 ], and selective sweep analysis [ 10 ], researchers have identified a smaller yet more plausible set of causal variants that affect complex traits [ 11 ]. However, comprehensive ASE analyses across multiple tissues have largely been limited to studies in humans [ 3 , 12 ]. While ASE can identify cis-regulatory elements [ 13 ], epigenomics remains the primary approach for uncovering these elements. This method has been extensively utilized in projects, such as the Encyclopedia of DNA Elements (ENCODE), which aims to identify causative genetic variants related to human diseases [ 14 ], and the Functional Annotation of Animal Genomes (FAANG), which seeks to improve the quantity and quality of agricultural products through molecular selection [ 15 ]. Some studies have started integrating ASE and epigenomic data to mutually validate findings [ 16 ] or to identify causal variants [ 17 ], but the relationship between ASE and epigenomic regulatory elements remains poorly understood. ASE also holds promise in unraveling the genetic basis of heterosis (hybrid vigor) [ 18 , 19 ]. The mechanisms behind heterosis have long been debated, primarily through two genetic models: dominance and overdominance, which are based on the degree of dominance of genetic variants, defined by the ratio of dominant to additive effect [ 20 ]. Although the relationship between dominant effects and ASE events has been discussed [ 18 , 19 ], comprehensive and quantitative analyses are still needed. In this study, we present a novel pipeline that uniformly integrates 7,532 public RNA-seq datasets to identify ASE variants in cattle. Using these variants, we performed three key analyses: (1) characterization of ASE events, including tissue specificity and directionality; (2) exploration of the relationship between ASE and dominant cis-eQTLs linked to molecular phenotypes of gene expression and dominant QTLs linked to physiological phenotypes of milk-related traits in Holstein cattle; and (3) functional annotation of ASE variants using epigenomic data, QTL databases, and conservation scores of sequences. Our findings not only advance the understanding of the genetic mechanisms underlying ASE events and heterosis but also provide valuable genetic resources for genomics and breeding programs in cattle. Results ASE discovery We downloaded 8,518 public RNA-seq samples from the NCBI SRA database, generating approximately 363 billion clean reads. After data cleaning (Materials and Methods), we retained 7,532 samples across 157 subsets, which spanned 35 distinct tissues and 126 NCBI BioProjects (Fig. 1 A and Additional file 1: Table S1 ). Among these subsets, 57 were well-defined as single-breed samples (Additional file 2: Fig. S1 ). The median sample size across these tissues was 80, with a range from 20 in ovary to 1,688 in whole blood. Based on the transcripts per million (TPM) values for 25,365 expressed genes (mean TPM > 0.1), hierarchical cluster analysis recapitulated tissue types. Unlike previous human GTEx studies [ 3 ], which identified blood-related samples as the primary outgroup, our research revealed the inclusion of mammary gland and milk somatic cell within the outgroup, underscoring the significance and similarity between these two tissues (Additional file 2: Fig. S2 ). A neighbor-joining (NJ) tree of imputed genotype primarily reflected two ancestries ( Bos taurus and Bos indicus ) (Additional file 2: Fig. S3). The observed separations of tissue type in gene expression and the distinctions of ancestry in genotype are consistent with findings from the Cattle Genotype-Tissue Expression (CattleGTEx) atlas [ 4 ], suggesting the high quality and utility of our data processing for follow-up analyses. To mitigate heterogeneity across different BioProjects, our ASE analysis focused on subsets with at least 20 individuals per tissue within each BioProject. We used the GATK ASEReadCounter module to calculate allele counts at heterozygous sites for each individual. Subsequently, we applied stringent filtering based on read counts and allele ratios, conducted a binomial test with false discovery rate (FDR) adjustment at the individual level, and further filtered by the number of heterozygous sites and the ratio of significant allelic imbalance at the population level. Through these processes, we identified a total of 452,052 ASE variants. After removing redundancies across subsets, we identified 161,059 unique ASE variants (Fig. 1 ). These variants were annotated to 13,136 genes, accounting for 66.1% of all protein-coding genes, which is lower than the 94.7% reported in human eQTL studies [ 21 ]. Genes lacking ASE variants were enriched for those lacking expression in our analyzed tissues, such as those involved in keratin filament and acrosomal vesicle (Additional file 2: Table S2 ). ASE variants were primarily distributed across exons, introns, 3' untranslated regions (UTRs), and intergenic regions, collectively accounting for 92.33% of the total, with 45,835 in exons, 39,590 in introns, 35,991 in 3' UTRs, and 27,302 in intergenic regions (Additional file 2: Table S3). The distribution pattern of ASE variants shows a similar trend to that of eQTLs in humans [ 21 ] and pigs [ 5 ], with protein-coding regions (e.g., stop gain, stop loss, nonsynonymous and synonymous mutation) showing the highest enrichment (Fig. 2 A). This pattern is further supported by the observation that more than 90% of the ASE variants were located downstream of the transcriptional start site (TSS) of the nearest gene (Fig. 2 B), reflecting the transcriptomic nature of ASE detection. Furthermore, we found a higher enrichment of variants affecting mRNA splicing, consistent with previous eQTL studies in humans [ 21 ] and pigs [ 5 ]. Additionally, we observed a more significant enrichment in the downstream compared to upstream regions and in the 3′ UTR compared to the 5′ UTR of protein-coding genes (Fig. 2 A). This contrasts with findings from humans [ 3 ] and cattle [ 4 ] eQTL studies, underscoring the distinct characteristics of ASE events. Among the identified ASE variants, 72.8% were unique to a single tissue, 14.2% were found in two tissues, and 5.3% appeared in three tissues (Fig. 2 C). To reduce heterogeneity across different BioProjects, we compared ASE variants within each BioProject, revealing that on average, 68.4% of ASE variants (ranging from 8.15–98.7%) were tissue-specific (Additional file 2: Fig. S4 and S5A). Regarding the shared expression of heterozygosity across tissues, an average of 40.0% of ASE variants (ranging from 6.3–71.3%) are exclusive to one tissue (Additional file 2: Fig. S5B), while an average of 53.3% of ASE variants (ranging from 22.6–87.9%) were replicated across different tissues (Additional file 2: Fig. S6). The distribution curve of shared ASE variants among tissues exhibited an inverse S-shaped pattern (Fig. 2 C), contrasting with the U-shaped pattern typically observed in eQTL studies in humans [ 3 ], cattle [ 4 ], and pigs [ 5 ]. These results suggest that ASE events exhibit a higher degree of tissue specificity. Additionally, we calculated the Pearson correlation coefficient of effect size for shared ASE variants across different tissues, yielding a median value of 0.7549 (ranging from − 0.6731 to 0.9941) for the 7,585 significant correlations out of 8,342 pairwise comparisons (Additional file 2: Fig. S7). This indicates that most shared ASE variants across tissues are likely driven by the same underlying factors. The discovery of ASE variants showed a significant correlation with sample size (Pearson r = 0.55; two-sided Student’s t-test: P = 6.61 × 10 − 14 ) (Fig. 2 D), a correlation lower than those reported for eGenes (0.85) and sGenes (0.63) in cattle [ 4 ], underscoring the unique nature of ASE events. A similar correlation was observed between the number of ASE variants and heterozygosity (Pearson r = 0.54; two-sided Student’s t-test: P = 2.09 × 10 − 13 ) (Fig. 2 E). However, substantial variability in the number of ASE variants was observed among immune-related tissues, further emphasizing the tissue specificity of ASE events. The average replication rate of ASE variants between different BioProjects of the same tissue was 49.1%, ranging from 1.9–100%. Notably, the replication rate was lower in reproductive tissues, suggesting a higher complexity in these tissues (Fig. 2 F). Despite the lower replication rate, the median Pearson correlation coefficient of effect size for shared ASE variants was relatively high (0.85, ranging from 0.21 to 0.99) in the 707 significant correlations out of 740 pairwise comparisons of the same tissue between different BioProjects (Additional file 2: Fig. S8), indicating that the majority of detected ASE variants were likely genuine. In our study, 80% of ASE variants exhibited a more than twofold effect on gene expression (Fig. 2 G), which is higher than the 22% observed in human cis -eQTL studies [ 21 ]. This underscores the importance and directness of ASE events in identifying regulatory elements. Interestingly, ASE variants in tissues related to reproduction and immune systems exhibited higher effect sizes (Additional file 2: Fig. S9), suggesting a potential heterozygous advantage in reproductive and immune traits. Furthermore, we identified 23,061 ASE variants with opposite directionality, accounting for 14.3% of the total. Among the 35 tissues examined, immune-related tissues, including milk somatic cells, whole blood, and white cells, exhibited a higher prevalence of ASE variants with opposite directionality (Fig. 2 H). Notably, the variants rs208685250 of MHC Class I JSP.1 and rs136860823 of MHC Class I BOLA demonstrated the highest number of subsets with opposite ASE directionality (Additional file 2: Fig. S10). This phenomenon of ASE variants with opposite directionality within the MHC family has also been documented in humans [ 22 ]. Gene ontology (GO) enrichment analysis of genes with ASE variants exhibiting opposite directionality highlighted innate immune response as the most significantly enriched biological process (Additional file 2: Table S4). Additionally, we observed a Pearson correlation coefficient of 0.86 between the number of ASE variants with opposite directionality and the variance explained by the first PEER factor (two-sided Student’s t-test: P = 3.03 × 10 − 6 ) (Fig. 2 I), indicating a strong influence of data heterogeneity on ASE directionality. Moreover, the absolute effect size of ASE variants with opposite directionality was significantly lower than that of variants with consistent directionality (two-sided Student’s t-test: P < 2.2 × 10 − 16 ) (Fig. 2 J). These results suggest that ASE variants with opposite directionality, particularly in immune-related genes, are more resilient. We also calculated linkage disequilibrium to examine allelic heterogeneity in gene expression associated with ASE events. Our analysis revealed that 10% of ASE genes contained more than one independent ASE variant (r 2 > 0.6) in subsets with sample sizes exceeding 80 (Fig. 2 K). This proportion is lower than the 46% reported for eGenes [ 4 ], suggesting that ASE variants may involve less complex interactions compared to eQTLs. Relationship between ASE variants and eQTL To better understand the relationship between ASE and eQTL, we performed an extensive analysis of cis -eQTL and cis -sQTL using datasets comprising over 80 individuals. This analysis, based on imputed genotype data from RNA-seq, identified an average of 79,814 additive cis -eQTL variants (ranging from 164 to 315,589), 29,958 dominant cis -eQTL variants (ranging from 342 to 116,913), 112,686 additive cis -sQTL variants (ranging from 974 to 375,496), and 50,195 dominant cis -sQTL variants (ranging from 378 to 220,710) (Additional file 2: Table S5). Among the identified ASE variants, we found that, on average, 282 were shared with additive cis -sQTL variants, 156 with dominant cis -eQTL variants, 410 with additive cis -sQTL variants, and 265 with dominant cis -sQTL variants. These correspond to average fold enrichments of 0.4136 (ranging from 0 to 1.0321), 0.6306 (ranging from 0 to 2.8491), 0.4357 (ranging from 0 to 1.2773), and 0.6089 (ranging from 0 to 1.6658). Although most of these enrichments were not statistically significant, the consistently higher enrichment observed in dominant cis -eQTLs and cis -sQTLs compared to their additive counterparts remains noteworthy. Specifically, in 13 out of the 17 subsets where ASE variants overlapped with dominant cis -eQTL variants, the enrichment of ASE variants was higher within dominant cis -eQTLs compared to their additive counterparts. Similarly, in 15 out of the 18 subsets where ASE variants overlapped with dominant cis -sQTL variants, the enrichment was greater in dominant cis -sQTLs relative to additive cis -sQTLs. These findings underscore the critical role of dominance effects in ASE regulation and highlight the importance of considering these effects when investigating the genetic mechanisms underlying ASE (Fig. 3 A and Additional file 2: Table S5). To further investigate the link between ASE and dominance, we calculated the Pearson correlation coefficient between the effect sizes of shared ASE and eQTL variants. The correlation between ASE and dominant cis -eQTL variants (Pearson r = 0.5121; two-sided Student’s t-test: P < 2.2 × 10 − 16 ) was notably higher than that between ASE and additive cis -eQTL variants (Pearson r = 0.3752; two-sided Student’s t-test: P < 2.2 × 10 − 16 ) (Fig. 3 B and 3 C). We also assessed the degree of dominance by computing the ratio of T-statistics ( \(\:{t}_{Dom}/{t}_{Add}\) ) for each eQTL, where \(\:{t}_{Dom}\) and \(\:{t}_{Add}\) represent the dominant and additive effects, respectively. Interestingly, the correlation between the effect sizes of shared ASE variants and the \(\:{t}_{Dom}/{t}_{Add}\) ratio (Pearson r = 0.5210; two-sided Student’s t-test: P < 2.2 × 10 − 16 ) (Fig. 3 D) was even higher than the correlation between shared ASE and either dominant or additive cis -eQTL variants individually. Additionally, in five out of six subsets with significant correlations of ASE effect sizes, including two subsets with a single breed, the \(\:{t}_{Dom}/{t}_{Add}\) ratio was higher than the correlations with either additive or dominant cis -eQTL variants alone (Additional file 2: Fig. S11). A similar pattern emerged when examining the relationship between shared ASE and cis -sQTL variants. The correlation between the effect sizes of shared ASE variants and the \(\:{t}_{Dom}/{t}_{Add}\) ratio (Pearson r = 0.9291; two-sided Student’s t-test: P < 2.2 × 10 − 16 ) was higher than that between ASE and dominant cis -sQTL variants (Pearson r = 0.7227; two-sided Student’s t-test: P < 2.2 × 10 − 16 ). Moreover, the correlation between ASE and dominant cis -sQTL variants was higher than that between ASE and additive cis -sQTL variants (Pearson r = 0.7092; two-sided Student’s t-test: P < 2.2 × 10 − 16 ) (Fig. 3 E, 3 F, and 3 G). Notably, in all 11 subsets with significant correlations of ASE effect sizes, including four subsets with a single breed, the \(\:{t}_{Dom}/{t}_{Add}\) ratio was higher than the correlations with either additive or dominant cis-sQTL variants alone (Additional file 2: Fig. S11). Functional annotation of ASE variants The neutral theory of molecular evolution posits that most mutations at the molecular level are neutral, with minimal or no phenotypic impact. However, mutations occurring in conserved genomic regions can have significant phenotypic effects [ 23 ]. To explore the relationship between ASE variants and sequence conservation, we employed PhastCons and PhyloP scores from the UCSC Genome Browser [ 24 ]. Our analysis identified 16,746 ASE variants with phyloP scores greater than 2.0 and 14,307 ASE variants with phastCons scores above 0.8, representing approximately a twofold enrichment in conserved bases (Fig. 4 A). Additionally, we observed that 93.57% of derived alleles among ASE variants exhibited negative effects (Additional file 2: Fig. S12), suggesting that most of these derived alleles are likely deleterious. ASE analysis is a powerful approach for identifying cis -regulatory elements involved in the regulation of diseases and other complex traits [ 13 ]. These regulatory elements are predominantly located within accessible chromatin regions, which can be identified using Assay for Transposase Accessible Chromatin Sequencing (ATAC-seq) [ 25 ]. To further investigate the relationship between ASE variants and chromatin accessibility, we analyzed ATAC-seq peak signals across 241 samples from 20 tissues obtained from the NCBI SRA database. The median sample size across these tissues was four, with a range from three samples in the hypothalamus, embryonic stem cells, spleen, and testes to 111 samples in the embryo (Additional file 1: Table S6). This analysis revealed that 50,645 ASE variants overlapped with ATAC-seq peaks, accounting for 31.44% of the total ASE variants. The proportion of ASE variants within ATAC-seq peaks ranged from 2.86–13.39% across the 11 tissues with both RNA-seq and ATAC-seq data (Fig. 4 B). In addition to ATAC-seq, Chromatin Immunoprecipitation Sequencing (ChIP-seq) is another method for globally identifying regulatory elements. ChIP-seq involves immunoprecipitation to selectively enrich DNA fragments associated with specific proteins, such as transcription factors or histone modifications. We examined the relationship between ASE variants and ChIP-seq peak signals across 158 samples from 20 tissues using five different antibodies. The median sample size across these tissues was three, with a range from two samples in the placenta and six other tissues to 103 samples in the mammary gland. The number of samples corresponding to the antibodies CTCF, H3K27ac, H3K27me3, H3K4me1, and H3K4me3 were 49, 90, 50, 153, and 157, respectively (Additional file 1: Table S7). Our analysis identified 90,556 ASE variants located within ChIP-seq peaks, accounting for 56.2% of all ASE variants. Specifically, the number of ASE variants within ChIP-seq peaks for the antibodies CTCF, H3K27ac, H3K27Me3, H3K4Me1, and H3K4Me3 were 45,828, 81,190, 57,010, 79,365, and 55,046, respectively. The proportion of ASE variants in ChIP-seq peaks varied from 0.7–58.7% across the eight tissues with both RNA-seq and ChIP-seq data. Notably, among the five antibodies, the proportions of ASE variants located within peaks associated with H3K4me1 and H3K27ac, both markers of active enhancers [ 26 ], were substantially higher (Fig. 4 C), indicating that ASE events are more frequently observed in enhancer regions. Furthermore, we explored the relationship between detected ASE variants and QTL data from the cattle QTL Database [ 27 ] to investigate the potential impact of ASE variants on phenotypes. This analysis identified 7,196 ASE variants within QTL regions. When evaluating the enrichment of ASE variants relative to the QTL database, we found that traits with medium or low heritability, such as meat and carcass traits, as well as reproduction and health traits, exhibited higher levels of enrichment. In contrast, traits with medium or high heritability, including milk and exterior traits, showed lower levels of enrichment (Fig. 4 D). Relationship between ASE variants and complex traits The primary goal of the ASE atlas is to serve as a resource for elucidating the genetic mechanisms underlying complex traits. By focusing on shared variants between ASE and QTLs for milk traits (milk yield, milk fat percentage, milk fat yield, milk protein percentage, milk protein yield, and somatic cell score) in the cattle QTL Database [ 27 ], as well as SNPs associated with milk traits reported in a previous study [ 28 ], we designed multiplex PCR experiments to validate the relationship between milk traits and ASE variants in Holstein cattle (Fig. 5 A and Materials and Methods). Of the 161 SNPs that were designed and successfully detected, 155 were retained after filtering for a minor allele frequency threshold of < 0.05 in a cohort of 1,052 cows. Using a genome-wide significant threshold ( P < 5 × 10 − 8 ) [ 29 ] and a mixed linear model, we identified 28 SNPs with additive effects for six traits and 13 SNPs with dominant effects for five traits, including 12 SNPs with both additive and dominant effects (Fig. 5 B and 5 C). Interestingly, the most significant locus was located on chromosome 14, a well-known locus responsible for milk traits, with the genes CPSF1, GPAA1 , PLEC , FAM83H , and LYNX1 mapped to this locus. Other candidate genes identified in our analysis include CSN2 , GHR , NUCB2 , ELFN2 , TBC1D22A , LGALS12 , and EP300 (Additional file 2: Table S8). Similar to eQTLs, we also investigated the relationship between the effect size of ASE variants and dominant QTLs, along with the degree of dominance. We found a strong, statistically significant correlation between the effect size of ASE variants and the \(\:{t}_{Dom}/{t}_{Add}\) ratio (Pearson r = 0.7459; two-sided Student’s t-test: P = 0.005) (Fig. 5 D). Even with a suggestive threshold ( P < 1 × 10 − 6 ) (Fig. 5 E), the correlation remained significant (Pearson r = 0.6843; two-sided Student’s t-test: P = 0.002). Unfortunately, owing to the limited number of QTLs identified for the same traits, the correlation between ASE effect sizes and dominant QTLs was not statistically significant. Discussion In this study, we have constructed a comprehensive ASE Atlas, representing a valuable resource of regulatory variants across 35 distinct cattle tissues. Our approach involves an in silico protocol for ASE analysis, encompassing the generation, characterization, and functional annotation using epigenomic data, as well as the analysis of genetic effects (additive and dominant) based on publicly available datasets. While our methodology is both time-efficient and cost-effective, comparable to the CattleGTEx atlas protocol [ 4 ], there are potential improvements. For instance, collecting RNA-seq samples from multiple tissues within the same individuals could help mitigate spatiotemporal heterogeneity, thereby enhancing comparability. Furthermore, integrating whole genome sequencing data with RNA-seq would allow for the generation of personalized genomes, which could reduce alignment bias [ 16 ]. It would also enable the phasing of haplotypes, thereby increasing the power of ASE analysis [ 12 ], and facilitate the identification of intergenic or non-expressed variants, thereby providing a more precise understanding of the relationship between ASE and QTLs related to both physiological and molecular phenotypes. Additionally, matching epigenomic data with RNA-seq would offer deeper insights into the mechanisms driving ASE events. Our findings highlight several distinctive features of ASE compared to eQTLs [ 3 – 5 , 21 ], particularly in terms of localization and tissue specificity. We observed a higher proportion of ASE variants in protein-coding genes and a higher degree of tissue specificity. This distinction likely arises because ASE variants are detected only in expressed sequences, whereas common variants in promotor regions are underrepresented. Specifically, only 1% of ASE variants were located in the promotor regions of protein-coding genes (within 2 kb upstream of the TSS), with 93% found downstream of the TSS. This contrasts with the ~ 60% of eQTLs found upstream of the TSS [ 3 ]. Additionally, the reproducibility of ASE variants was lower than that of pig eQTLs [ 5 ] in independent datasets (49% vs. 77%) and across different tissues within the same BioProject (53% vs. 92%). These observations suggest that ASE may exhibit greater specificity than eQTLs. Supporting this, we found that ASE variants are less frequent in promotor regions of protein-coding genes and more commonly located within ChIP-seq peaks marked by active enhancers (H3K4me1 and H3K27ac) rather than active promoters (H3K4me3) [ 26 ]. Moreover, ASE variants were less commonly found in ATAC-seq peaks, which are more typically associated with promoters [ 30 , 31 ], compared to ChIP-seq peaks marked by active enhancers. These results suggest that enhancer mutations may be a primary source of ASE events. The future availability of matched RNA-seq and H3K27ac histone ChIP-seq data will likely facilitate the identification of causative mutations in enhancer RNA [ 32 ], potentially driving ASE. When examining the functional enrichment of ASE variants, we observed a significant presence of variants affecting splicing, consistent with our finding that ASE variants are more closely related to sQTLs than eQTLs. The highest enrichment was observed in variants leading to an in-frame stop codon. These findings strongly suggest that mutations altering mRNA splicing are key mechanisms underlying ASE [ 11 ]. Similar to eQTLs in humans [ 21 ] and pigs [ 5 ], ASE variants were also enriched in nonsynonymous, synonymous, stop-loss, and UTR variants, indicating the involvement of post-transcriptional and translational mechanisms, such as altered RNA-binding proteins, RNA stability, RNA editing, and the creation or disruption of upstream initiation codons [ 11 , 33 ]. Our study also showed that ASE variants are more enriched in dominant cis -eQTLs and cis -sQTLs compared to additive ones, with a higher correlation coefficient between ASE variant effect sizes and those of dominant cis -eQTLs/ cis -sQTLs. Previous studies in maize reported that dominant expression patterns were more prevalent than additive patterns in ASE genes [ 19 ], and in poplar, 53% of ASE variants exhibited dominant effects on physiological and photosynthetic traits [ 34 ]. While research in rice suggests that biased expression of the favorable allele leads to dominant effects [ 18 ], our findings indicate a more complex relationship, where the degree of dominance in physiological and molecular phenotypes correlates with ASE effect sizes. Among the identified ASE variants, 14.3% exhibited opposite directional effects, particularly in immune-related tissues, genes, and GO terms, such as the MHC family. ASE variants with opposite directional effects in the MHC gene family have also been reported in human T cells [ 22 ]. Indeed, immune response eQTLs with opposite directional effects have also been observed in different contexts in humans [ 35 ]. Our study also found the number of ASE variants with opposite directional effects was closely related to dataset heterogeneity. Interestingly, ASE variants with opposite directional effects had significantly lower absolute effect sizes than those with consistent directional effects, suggesting that smaller effect sizes may facilitate directional shifts, thereby enhancing resilience in varying environments. This hypothesis is further supported by the relatively high variability in the number of ASE variants in immune-related tissues. Similar findings in maize suggest that ASE variants with opposite directional effects may contribute to adaptation in diverse environments [ 19 ], challenging the hypothesis that direction-shifting ASE causes overdominance in rice [ 18 ]. Although ASE has been widely used to explore the mechanism of heterosis [ 18 , 19 ], to our knowledge, this is the first study to report a linear correlation between ASE and the degree of dominance implicated in heterosis models [ 36 , 37 ]. We also found that tissues related to immune and reproductive functions, which are associated with high levels of heterosis [ 38 ], exhibited higher ASE effect sizes. Consistently, ASE variants were more enriched in functional traits with high levels of heterosis and low heritability [ 39 ] in the cattle QTL database. The predominance of derived alleles with negative effects supports the hypothesis that the dominance of advantageous ancestral alleles complements the deleterious effects of derived alleles [ 37 ]. Based on our discovery of a linear relationship between ASE and dominance, along with the hypothesis that overdominance results from allele interactions at a single locus [ 40 ], we propose that overdominance in physiological phenotypes stems from the interaction of two alleles at the molecular level, with certain loci manifesting overdominant effects. To validate our hypothesis inferred from the molecular phenotype of gene expression, we also observed the linear correlation between ASE effect size and the degree of dominance in QTLs affecting milk-related traits. Notably, the well-known QTL harboring CPSF1 on chromosome 14, affecting milk yield [ 41 , 42 ], exhibited both additive and dominant effects in our study. In fact, this dominant QTL was also observed in a previous study [ 43 ]. Other well-established candidate genes associated with milk traits, including CSN2 [ 41 , 44 , 45 ], GHR [ 43 , 46 ], and NUCB2 [ 47 ], also exhibited dominant effects. Additionally, genes with limited or no previous reports in association with milk traits, such as EP300, TBC1D22A , ELFN2 , and LGALS12 , also displayed dominant effects. These results supported the hypothesis that dominance is pervasive in mammals [ 36 ]. Compared with previous studies, and thanks to the homogeneity of ASE data in RNA-seq, the variants we identified in association with milk traits are closer to causal. Conclusions Our comprehensive ASE atlas offers a valuable resource of regulatory variants for uncovering the genetic mechanisms underlying complex traits and enhancing economic traits through molecular selection. The characterization and functional annotation of ASE variants deepen our understanding of ASE formation and aid in identifying causal variants associated with complex traits. The observed correlation between ASE and dominance effects provides new insights into the genetic basis of heterosis, with promising applications in cattle production. Materials and Methods Data origin, alignment, and clustering A total of 33,999 RNA-seq entries were retrieved from the NCBI SRA database (4 July 2023) by searching for “cattle” and selecting “RNA” as the source. An in-house Perl script was used to refine this dataset by filtering for biological samples with the assay type of “RNA-seq” and the organisms of “ Bos taurus , Bos indicus , or Bos taurus × Bos indicus .” Further filtering was applied to include only samples sequenced on either the ILLUMINA or BGISEQ platforms. We also required a minimum of 20 samples per well-defined tissue type within each BioProject. This curation process yielded a final set of 8,518 transcriptomic datasets. The raw reads were processed using Trimmomatic v0.39 [ 48 ] with the parameters “LEADING:3 TRAILING:3 SLIDINGWINDOW:4:15 MINLEN:36” to ensure high-quality reads. The maximum read lengths were determined using seqkit v2.5.1 [ 49 ]. The filtered reads were aligned to the Bos taurus reference genome assembly (ARS-UCD1.2_Btau5.0.1Y) using STAR v2.7.10b [ 50 ] with the parameters “--outFilterMultimapNmax 1 --outFilterIntronMotifs RemoveNoncanonical Unannotated --outFilterMismatchNmax 10 --outSAMstrandField intronMotif --outSJfilterReads Unique.” Unaligned reads were further aligned using HISAT2 v2.2.1 [ 51 ]. The alignment files from STAR and HISAT2 were merged using Picard v3.0.0, and alignment statistics were computed using Samtools v1.17 (stat command) [ 52 ]. Samples were retained for downstream analysis if they met the criteria of a read mapping ratio ≥ 0.6 and a minimum of 50,000,000 mapped bases. We further filtered for sample size, requiring a minimum sample size of ≥ 20 per clearly defined tissue within each BioProject, resulting in a final dataset of 7,532 samples. Gene expression levels were quantified as TPM using Stringtie v2.2.1 [ 53 ], with the GFF annotation file of the reference assembly. Hierarchical clustering of autosomal gene expression levels with an average TPM greater than 0.1 was performed using the R v4.2.3 package dendextend [ 54 ], employing Euclidean distance. SNP calling and imputation To ensure high-quality variant calling, reads with a mapping quality score below 1 and unmapped reads were excluded using Samtools v1.17 (view command) [ 52 ]. The filtered alignment files were then processed with Picard tools, utilizing the ReorderSam, SortSam, and AddOrReplaceReadGroups commands. The SplitNCigarReads function in GATK v4.4.0.0 [ 55 ] was used to split reads at exon junctions and hard-clip overhanging sequences within intronic regions. Base quality recalibration was performed using GATK’s BaseRecalibrator and ApplyBQSR modules. Variant calling was conducted using GATK’s HaplotypeCaller, followed by joint genotyping using the GenomicsDBImport and GenotypeGVCFs modules. SNPs were filtered using the VariantFiltration tool with the expression “QD 30.0.” The SelectVariants module of GATK was used to extract biallelic SNPs with a minor allele frequency of 0.01. Clean SNPs across all autosomes were merged using the MergeVcfs module of GATK and subsequently annotated using ANNOVAR v2016-02-01 [ 56 ]. To further ensure data quality, SNPs were additionally filtered using VCFtools v0.1.17 [ 57 ] with the parameters “--maf 0.05 --max-missing 0.7 --minDP 10 --minGQ 20.” This filtering was applied both to the entire dataset of 7,532 samples and to 19 subsets, each containing more than 80 samples. Genotype imputation was then performed using Beagle v27Jan18.7e1 [ 58 ], leveraging a reference SNP panel [ 59 ]. After imputation, SNPs were filtered using BCFtools v1.17 [ 60 ] with thresholds of MAF > 0.05 and DR2 ≥ 0.8 to maintain high data quality. If imputed genotypes from the entire dataset of 7,532 samples were absent in the 19 subsets, they were incorporated into the respective subsets for further eQTL and related analyses. Hamming distances were calculated on the imputed genotypes of the entire dataset using PLINK v1.90b7.1 [ 61 ]. A neighbor-joining (NJ) tree was constructed based on these distances using MEGA11 [ 62 ] and then visualized using iTOL v6 [ 63 ]. Linkage disequilibrium (LD) analysis for each dataset was also conducted using PLINK v1.90b7.1. ASE detection The procedure for detecting ASE was adapted from a previous study [ 64 ] with slight modifications. Based on SNP data from all 7,532 samples, the ASEReadCounter module of GATK v4.4.0.0 was employed to count the reads corresponding to each allele at heterozygous sites for each individual. To be included in the analysis, both the reference and alternative allele counts were required to exceed three, and the minor allele read count ratio had to be greater than 0.01. A binomial test was then applied to these sites, followed by Benjamini-Hochberg multiple testing correction to control the FDR. Significant allelic imbalance at the individual level was determined by an FDR threshold of 0.05. At the population level, we focused on sites where at least six individuals were heterozygous. Sites showing significant allelic imbalance (FDR < 0.05) in over 90% of these individuals were retained for further analysis. The effect size of allelic imbalance was quantified using the log allelic fold change (aFC), calculated according to a modified model from a previous study [ 65 ]. The calculations were performed as follows: (1) when the direction of allelic imbalance was consistent across individuals: $$\:{\delta\:}_{\text{1,0}}=\begin{array}{c}median\\\:n=1...N\end{array}\frac{{c}_{1,n}}{{c}_{0,n}}$$ The reported effect size: $$\:{s}_{\text{1,0}}={{log}}_{2}{\delta\:}_{\text{1,0}}$$ (2) when the directions of allelic imbalance varied among individuals: $$\:{\delta\:}_{\text{1,0}}=\left|{log}_{2}\left(\frac{{c}_{1,n}}{{c}_{0,n}}\right)\right|$$ The reported effect size: $$\:{s}_{\text{1,0}}=\begin{array}{c}median\\\:n=1...N\end{array}{\delta\:}_{\text{1,0}}$$ In these equations, \(\:{c}_{1,n}\) represents the read count for the alternative allele, \(\:{c}_{0,n}\) represents the read count of the reference allele, \(\:n\) is the index of the heterozygous site, and \(\:N\) is the number of heterozygous sites showing significant allelic imbalance in the population. The distance from ASE variants to the TSS was defined based on annotations from ANNOVAR and the GFF annotation file of the reference assembly, measuring the distance from each ASE variant to the TSS of the nearest genes. GO enrichment analysis was conducted using the DAVID server [ 66 ]. Covariate analysis for eQTL discovery For the eQTL analysis, we employed a covariate analysis pipeline based on the methodology outlined in the CattleGTEx project [ 4 ]. First, gene expression levels for 25,365 expressed genes, each with a mean TPM value greater than 0.1, were normalized using the quantile-quantile normalization method [ 67 ]. To identify hidden factors contributing to transcriptome-wide variation in gene expression, Bayesian methods were employed to estimate latent covariates using PEER v1.0 [ 68 ]. The first ten PEER factors, where the cumulative posterior variances reached or nearly reached a plateau [ 4 ], were included as covariates. Additionally, to account for population structure in eQTL analysis, principal component (PC) analysis was performed based on the imputed genotypes using PLINK v1.90b7.1 [ 61 ]. Following the recommendations of CattleGTEx [ 4 ], the number of PCs included as covariates was determined by sample size: the first three PCs for sample sizes less than 150, the first five PCs for sample sizes between 150 and 249, and the first ten PCs for sample sizes of 250 or greater. Cis -eQTL mapping For cis -eQTL mapping, normalized gene expression levels were used as the molecular phenotype, focusing on genes with a mean TPM value greater than 0.1. The analysis was conducted across 19 subsets, each containing more than 80 samples, using QTLtools v1.2 [ 69 ] with the parameters “--nominal 0.01.” To account for confounding factors, PEER factors and principal components were included as covariates in the cis-eQTL detection. A significance threshold of P < 1 × 10 − 6 was applied for SNP associations. When evaluating additive effects, genotypes were coded as follows: homozygous reference (0), heterozygous (1), and homozygous alternative (2). For dominant effects, the genotypes were coded as homozygous reference (0), heterozygous (1), and homozygous alternative (0). The effect size of cis -eQTL was described using log allelic fold change (aFC), calculated according to a modified linear regression model [ 65 ]: $$\:y={\beta\:}_{0}+{\beta\:}_{1}{X}_{1}+{\beta\:}_{2}{X}_{2}+\epsilon\:$$ Where \(\:y\) denotes the normalized molecular phenotype, \(\:{\beta\:}_{0}\) represents the intercept, \(\:{\beta\:}_{1}\) captures the effect of SNP markers, \(\:{X}_{1}\) represents the maker genotypes, \(\:{\beta\:}_{2}\) accounts for the effect of covariates, \(\:{X}_{2}\) represents the covariates, \(\:\epsilon\:\) represents the random residuals. The effect size ( \(\:{s}_{\text{1,0}}\) ) was reported as: $$\:{\delta\:}_{\text{1,0}}=\frac{\:2{\beta\:}_{1}}{{\beta\:}_{0}}+1$$ $$\:{s}_{\text{1,0}}={{log}}_{2}{\delta\:}_{\text{1,0}}$$ The degree of dominance was assessed by calculating the ratios \(\:{t}_{Dom}/{t}_{Add}\) according to established methods [ 20 ]. The T-statistics for additive and dominant effects, \(\:{t}_{Add}\) and \(\:{t}_{Dom}\) , were computed as follows: $$\:{t}_{Add}=\frac{{\beta\:}_{Add}}{se\left({\beta\:}_{Add}\right)}\:\:\:\:\:\:\:\:\:\:\:{t}_{Dom}=\frac{{\beta\:}_{Dom}}{se\left({\beta\:}_{Dom}\right)}\:$$ Here, \(\:{\beta\:}_{Add}\) and \(\:{\beta\:}_{Dom}\) represent the additive and dominant effects, respectively, while \(\:se\left({\beta\:}_{Add}\right)\) and \(\:se\left({\beta\:}_{Dom}\right)\) denote their corresponding standard errors. If the effect of the minor allele is negative, the degree of dominance is defined as \(\:-{t}_{Dom}/{t}_{Add}\) ; otherwise, it is defined as \(\:{t}_{Dom}/{t}_{Add}\) . Cis -sQTL mapping For cis -sQTL mapping, splicing junctions were identified from alignment files using the Leafcutter v0.2.7 tool [ 70 ]. Initially, the bam2junc.sh script was used to convert each individual’s BAM file into a junction file. The identified junction files were then clustered across the population using the leafcutter_cluster.py script, with a minimum of 50 reads per junction and a maximum intron length of 500,000 base pairs. Phenotype tables required for downstream splicing analysis were prepared using the prepare_phenotype_table.py script from the Leafcutter tool and subsequently normalized using quantile-quantile normalization. The covariates and analytical methods for cis -sQTL mapping were identical to those employed for cis -eQTL mapping. Overlap analysis of ASE variants with conservation sites Based on the multiple alignment format (maf) file from UCSC [ 24 ] ( https://hgdownload.soe.ucsc.edu/goldenPath/hg38/multiz470way/maf/ ), we established the correspondence between human autosomal physical positions and the cattle ARS-UCD1.2 assembly using an in-house Perl script. Conservation scores for humans were retrieved from UCSC ( https://hgdownload.soe.ucsc.edu/goldenPath/hg38/phyloP470way/hg38.470way.phyloP/ and https://hgdownload.soe.ucsc.edu/goldenPath/hg38/phastCons470way/hg38.470way.phastCons/ ) were subsequently translated to corresponding cattle positions using the established mappings. We then extracted the relevant conservation scores for ASE variants in cattle through another custom Perl script. To infer ancestral alleles, we obtained the whole-genome sequences of five outgroup species (American Bison, Banteng, Yak, European Bison, and Gaur) [ 71 ] (Additional file 1: Table S9). After performing quality control similar to the RNA-seq data processing, we aligned the trimmed reads to the Bos taurus reference assembly using bwa-mem v0.7.17-r1188 [ 46 ] with default parameters. The alignment files were sorted by coordinate using Picard's SortSam tool, and duplicate reads were marked and removed using Picard's MarkDuplicates tool. Variant calling was conducted using BCFtools v1.17 mpileup with parameters “-q 30 -C 50 -Q 20 -B.” The ancestral allele was inferred using a consensus sequence from at least two of the outgroup species through an in-house Perl script. ATAC-seq analysis We retrieved 528 entries of run information from the NCBI SRA database (14 December 2023) by searching for "cattle ATAC-seq,” specifying "illumina" as the platform, “DNA” as the source, “EpiGenomics” as the strategy, and “ Bos taurus ” as the organism. An in-house Perl script was used to further refine the dataset, applying criteria including a minimum sample size of three and clearly defined tissue information within each BioProject, resulting in 241 ATAC-seq data. The raw reads were quality-trimmed, and adapters were removed using Trim Galore v0.6.10 [ 72 ] with the parameters “--q 25 --phred33 --stringency 3 --length 35 -e 0.1.” The trimmed reads were aligned to the cattle reference genome (ARS-UCD1.2) using Bowtie2 v2.5.3 [ 73 ] with the parameters “--very-sensitive -X 2000.” Reads mapped to the sex chromosomes (X and Y) were filtered out. Duplicate reads were marked and removed using Sambamba v1.0.1 [ 74 ], followed by sorting the BAM files by coordinate and read name using Samtools v1.17 [ 52 ]. Peaks representing open chromatin regions were called using MACS3 v3.0.0 [ 75 ] with the parameters “--shift 75 --extsize 150 --nomodel --call-summits --nolambda --keep-dup all -q 0.01”. We then used the intersect command of bedtools v2.26.0 [ 76 ] to obtain overlapping peak regions from the same tissue type, which were subsequently utilized to identify ASE variants located within ATAC-seq peaks. ChIP-seq analysis A total of 5,251 entries of run information were retrieved from the NCBI SRA database (6 May 2024) by searching for "cattle ChIP-seq,” specifying "illumina" as the platform, “DNA” as the source, “EpiGenomics” as the strategy, and “ Bos taurus ” as the organism. Using an in-house Perl script, we selected samples from four BioProjects (PRJEB41939 [ 77 ], PRJEB52456 [ 78 ], PRJEB53044 [ 79 ] and PRJEB6906 [ 80 ]) that had at least two samples per tissue per antibody, resulting in 499 ChIP-seq datasets. Both input and antibody-treated samples were initially processed using Trim Galore v0.6.10 with the parameters “--q 25 --phred33 --length 25 -e 0.1 --stringency 4 --paired” to remove low-quality bases and adapters. The trimmed reads were aligned to the reference genome (ARS-UCD1.2) using Bowtie2 v2.5.3, and the resulting SAM files were converted to BAM format using Samtools v1.17. Duplicate reads were marked and removed using Sambamba v1.0.1. The deduplicated BAM files were then sorted using Samtools v1.17. To identify significant DNA-protein interaction regions, peaks were called using MACS3 v3.0.0 [ 75 ] with a stringent cutoff (-q 0.01) to ensure high-confidence detection. Antibody-treated samples were compared against input controls to filter out background signals and enhance the specificity of the detected peaks. As with ATAC-seq, we used the intersect command of bedtools v2.26.0 [ 76 ] to find overlapping peak regions from the same tissue and antibody type and then identified ASE variants located within these ChIP-seq peaks. Validation of shared variants between ASE and QTL for milk traits Among the 8,733 SNPs associated with milk traits reported by a previous study [ 28 ], 156 overlapped with ASE variants. Additionally, of the 15,958 SNPs linked to milk traits in the cattle QTL database, 104 overlapped with ASE variants. After merging these two sets, we retained 248 SNPs. Filtering out those with a minor allele frequency < 0.05 in Holstein cattle based on a reference panel from a previous study [ 59 ] reduced the number to 208 SNPs. Further pruning using PLINK with the parameter “--indep-pairwise 5 1 0.8” in the reference panel resulted in a final set of 171 SNPs. An additional 10 SNPs were discarded due to being flanked by repeat sequences as determined by Adsen Biotechnology Co., Ltd. (Urumchi, China) during genotyping. A total of 1,052 Chinese Holstein cows from Xinjiang were used in this study. Nine milk-related traits were considered: milk yield (kg), milk fat percentage (%), milk fat yield (kg), milk protein percentage (%), milk protein yield (kg), somatic cell count (10 4 /ml), milk lactose percentage (%), total solids percentage (%), and urea nitrogen concentration (mg/dl). The total number of test-day records was 8,839. Descriptive statistics of the test-day traits are provided in Additional file 2: Table S10. Whole blood was collected from the jugular vein of each individual using a blood collection needle by an experienced veterinarian and then placed in a blood collection tube containing EDTA. The genomic DNA was extracted using the standard phenol-chloroform procedure, and then its quantity and quality were assessed using agarose gel electrophoresis and a Qubit fluorometer (Invitrogen, Carlsbad, USA), respectively. The qualified DNA was transported to Adsen Biotechnology Co., Ltd. for multiplex PCR experiments, following sequencing and variant detection (Additional file 2: Supplementary Materials and Methods). We used the following mixed linear model, implemented with the lmerTest package in R [ 6 ]. $$\:y={\beta\:}_{0}+{\beta\:}_{1}{X}_{1}+{\beta\:}_{2}{X}_{2}+{\beta\:}_{3}{X}_{3}+\epsilon\:$$ Where \(\:y\) denotes the phenotype, \(\:{\beta\:}_{0}\) denotes the intercept of the linear regression, \(\:{\beta\:}_{1}\) represents the effects of SNP markers, and \(\:{X}_{1}\) denotes the maker genotypes. The term \(\:{\beta\:}_{2}\) represents the fixed effects of herd, parity, and milk month, while \(\:{X}_{2}\) denotes these factors. The term \(\:{\beta\:}_{3}\) represents the random effect of test day, with \(\:{X}_{3}\) denoting the test day, and \(\:\epsilon\:\) represents random residuals. Descriptive statistics for the fixed and random effects are provided in Additional file 2: Table S11. Similar to cis -eQTL analysis, for examining additive effects, genotypes were coded as follows: homozygous reference (0), heterozygous (1), and homozygous alternative (2). For dominant effects, genotypes were coded as homozygous reference (0), heterozygous (1), and homozygous alternative (0). The ratios of \(\:{t}_{Dom}/{t}_{Add}\) , used to measure the degree of dominance, were calculated as described previously [ 20 ] and were consistent with the cis -eQTL analysis. Abbreviations ASE Allelic-specific expression QTL quantitative trait loci eQTL expression QTL sQTL splicing QTL GWAS genome-wide association study GTEx Genotype-Tissue Expression ENCODE Encyclopedia of DNA Elements FAANG Functional Annotation of Animal Genomes TPM transcripts per million CattleGTEx Cattle Genotype-Tissue Expression UTR untranslated region TSS transcriptional start site ATAC-seq Assay for Transposase Accessible Chromatin sequencing ChIP-seq Chromatin Immunoprecipitation Sequencing LD Linkage disequilibrium MAF minor allele frequency NJ neighbor-joining FDR false discovery rate aFC allelic fold change GO Gene Ontology PC principal component. Declarations Acknowledgements We would like to express our sincere gratitude to Professor Tom Druet from University of Liège, Professor Yu Wang from Northwest A&F University and Associate Professor Han Xu from Anhui Agricultural University for their valuable suggestions and insightful contributions to this study. Authors’ contributions Q.C. conceptualized and designed the study. Q.C. and J.L. conducted the allele-specific expression analysis. X.L. carried out the ATAC-seq and ChIP-seq analyses. Q.C. and L.X. performed mixed linear model analysis on milk related traits. L.L. prepared the schematic diagram of cattle tissues. L.X. coordinated with the Holstein farm for phenotype data collection and blood sampling. X.H. secured funding and supervised the entire study. Q.C. wrote the manuscript. Funding This study was financially supported by the National Key R&D Program of China (2021YFD1200903). Data Availability All raw data analyzed in this study are publicly available for download without restrictions from the NCBI SRA database (https://www.ncbi.nlm.nih.gov/sra/). Details of RNA-seq, ATAC-seq, ChIP-seq and whole genome sequence can be found in Additional file 1: Table S1, S6, S7 and S9, respectively. All the computational scripts and codes for RNA-seq, ATAC-seq. and ChIP-seq data quality control, gene expression normalization, ASE identification, SNP detection, genotype imputation, cis -eQTL and cis -sQTL mapping, and functional annotation are available at the github website (https://github.com/Qiuming1986/ASE-in-cattle). All the original results of ASE, eQTL, sQTL, ATAC-seq, and Chip-seq have been deposited in the Figshare database with the following digital object identifier: 10.6084/m9.figshare.26808418. Ethics approval and consent to participate All animal procedures were conducted in accordance with the Regulations for the Administration of Affairs Concerning Experimental Animals of China and were approved by the Animal Care Committee of Xinjiang Agricultural University, which oversees the ethical use of animals in research at the university. Competing interests The authors declare that they have no competing interests. References Ward LD, Kellis M. Interpreting noncoding genetic variation in complex traits and human disease. Nature biotechnology. 2012;30(11):1095–106. Maurano MT, Humbert R, Rynes E, Thurman RE, Haugen E, Wang H, et al. Systematic localization of common disease-associated variation in regulatory DNA. Science. 2012;337(6099):1190–5. Consortium TG, Ardlie KG, DeLuca DS, Segrè AV, Sullivan TJ, Young TR, et al. The Genotype-Tissue Expression (GTEx) pilot analysis: multitissue gene regulation in humans. Science. 2015;348(6235):648–60. Liu S, Gao Y, Canela-Xandri O, Wang S, Yu Y, Cai W, et al. A multi-tissue atlas of regulatory variants in cattle. Nat Genet. 2022;54(9):1438–47. Epub 2022/08/12. doi: 10.1038/s41588-022-01153-5 . PubMed PMID: 35953587; PubMed Central PMCID: PMCPMC7613894. Teng J, Gao Y, Yin H, Bai Z, Liu S, Zeng H, et al. A compendium of genetic regulatory effects across pig tissues. Nature Genetics. 2024:1–12. Kuznetsova A, Brockhoff PB, Christensen RHB. lmerTest package: tests in linear mixed effects models. Journal of statistical software. 2017;82(13). Mapel XM, Kadri NK, Leonard AS, He Q, Lloret-Villas A, Bhati M, et al. Molecular quantitative trait loci in reproductive tissues impact male fertility in cattle. nature communications. 2024;15(1):674. Zou J, Hormozdiari F, Jew B, Castel SE, Lappalainen T, Ernst J, et al. Leveraging allelic imbalance to refine fine-mapping for eQTL studies. PLoS genetics. 2019;15(12):e1008481. Cavalli M, Pan G, Nord H, Arzt EW, Wallerman O, Wadelius C. Allele-specific transcription factor binding in liver and cervix cells unveils many likely drivers of GWAS signals. Genomics. 2016;107(6):248–54. Magris G, Jurman I, Fornasiero A, Paparelli E, Schwope R, Marroni F, et al. The genomes of 204 Vitis vinifera accessions reveal the origin of European wine grapes. Nature Communications. 2021;12(1):7240. Cleary S, Seoighe C. Perspectives on allele-specific expression. Annual Review of Biomedical Data Science. 2021;4(1):101–22. Castel SE, Aguet F, Mohammadi P, Ardlie KG, Lappalainen T. A vast resource of allelic expression data spanning human tissues. Genome biology. 2020;21:1–12. Delbare SY, Clark AG. Allele-specific expression elucidates cis-regulatory logic. PLoS Genetics. 2018;14(11):e1007690. Moore JE, Purcaro MJ, Pratt HE, Epstein CB, Shoresh N, Adrian J, et al. Expanded encyclopaedias of DNA elements in the human and mouse genomes. Nature. 2020;583(7818):699–710. Clark EL, Archibald AL, Daetwyler HD, Groenen MA, Harrison PW, Houston RD, et al. From FAANG to fork: application of highly annotated genomes to improve farmed animal production. Genome Biology. 2020;21:1–9. Quan J, Yang M, Wang X, Cai G, Ding R, Zhuang Z, et al. Multi-omic characterization of allele-specific regulatory variation in hybrid pigs. Nature Communications. 2024;15(1):5587. Kundu K, Tardaguila M, Mann AL, Watt S, Ponstingl H, Vasquez L, et al. Genetic associations at regulatory phenotypes improve fine-mapping of causal variants for 12 immune-mediated diseases. Nature genetics. 2022;54(3):251–62. Shao L, Xing F, Xu C, Zhang Q, Che J, Wang X, et al. Patterns of genome-wide allele-specific expression in hybrid rice and the implications on the genetic basis of heterosis. Proceedings of the National Academy of Sciences. 2019;116(12):5653-8. Zhan W, Cui L, Yang S, Zhang K, Zhang Y, Yang J. Natural variations of heterosis-related allele-specific expression genes in promoter regions lead to allele-specific expression in maize. BMC genomics. 2024;25(1):476. Cui L, Yang B, Pontikos N, Mott R, Huang L. ADDO: a comprehensive toolkit to detect, classify and visualize additive and non-additive quantitative trait loci. Bioinformatics. 2020;36(5):1517–21. Consortium TG, Aguet F, Anand S, Ardlie KG, Gabriel S, Getz GA, et al. The GTEx Consortium atlas of genetic regulatory effects across human tissues. Science. 2020;369(6509):1318–30. Gutierrez-Arcelus M, Baglaenko Y, Arora J, Hannes S, Luo Y, Amariuta T, et al. Allele-specific expression changes dynamically during T cell activation in HLA and other autoimmune loci. Nature genetics. 2020;52(3):247–53. Eyre-Walker A, Keightley PD. The distribution of fitness effects of new mutations. Nature Reviews Genetics. 2007;8(8):610–8. Raney BJ, Barber GP, Benet-Pagès A, Casper J, Clawson H, Cline MS, et al. The UCSC Genome Browser database: 2024 update. Nucleic Acids Research. 2024;52(D1):D1082-D8. Lu Z, Hofmeister BT, Vollmers C, DuBois RM, Schmitz RJ. Combining ATAC-seq with nuclei sorting for discovery of cis-regulatory regions in plant genomes. Nucleic acids research. 2017;45(6):e41-e. Shlyueva D, Stampfel G, Stark A. Transcriptional enhancers: from properties to genome-wide predictions. Nature Reviews Genetics. 2014;15(4):272–86. Hu Z-L, Park CA, Reecy JM. Bringing the animal QTLdb and CorrDB into the future: meeting new challenges and providing updated services. Nucleic acids research. 2022;50(D1):D956-D61. Jiang J, Cole JB, Freebern E, Da Y, VanRaden PM, Ma L. Functional annotation and Bayesian fine-mapping reveals candidate genes for important agronomic traits in Holstein bulls. Communications biology. 2019;2(1):212. Schaid DJ, Chen W, Larson NB. From genome-wide associations to candidate causal variants by statistical fine-mapping. Nature Reviews Genetics. 2018;19(8):491–504. Yuan C, Tang L, Lopdell T, Petrov VA, Oget-Ebrad C, Moreira GCM, et al. An organism-wide ATAC-seq peak catalog for the bovine and its use to identify regulatory variants. Genome Research. 2023;33(10):1848–64. Alexandre PA, Naval-Sánchez M, Menzies M, Nguyen LT, Porto-Neto LR, Fortes MR, et al. Chromatin accessibility and regulatory vocabulary across indicine cattle tissues. Genome biology. 2021;22:1–20. Wang C, Chen C, Lei B, Qin S, Zhang Y, Li K, et al. Constructing eRNA-mediated gene regulatory networks to explore the genetic basis of muscle and fat-relevant traits in pigs. Genetics Selection Evolution. 2024;56(1):1–21. Flynn ED, Lappalainen T. Functional characterization of genetic variant effects on expression. Annual Review of Biomedical Data Science. 2022;5(1):119–39. Xuan A, Song Y, Bu C, Chen P, El-Kassaby YA, Zhang D. Changes in DNA methylation in response to 6-benzylaminopurine affect allele-specific gene expression in Populus tomentosa. International Journal of Molecular Sciences. 2020;21(6):2117. Kim-Hellmuth S, Bechheim M, Pütz B, Mohammadi P, Nédélec Y, Giangreco N, et al. Genetic regulatory effects modified by immune activation contribute to autoimmune disease associations. Nature communications. 2017;8(1):1–10. Cui L, Yang B, Xiao S, Gao J, Baud A, Graham D, et al. Dominance is common in mammals and is associated with trans-acting gene expression and alternative splicing. Genome Biology. 2023;24(1):215. Chen ZJ. Genomic and epigenetic insights into the molecular bases of heterosis. Nature Reviews Genetics. 2013;14(7):471–82. Wakchaure R, Ganguly S, Praveen PK, Sharma S, Kumar A, Mahajan T, et al. Importance of heterosis in animals: a review. International Journal of Advanced Engineering Technology and Innovative Science. 2015;1(2):1–5. Getahun D, Alemneh T, Akeberegn D, Getabalew M, Zewdie D. Importance of hybrid vigor or heterosis for animal breeding. Biochemistry and Biotechnology Research. 2019;7:1–4. Birchler JA, Yao H, Chudalayandi S. Unraveling the genetic basis of hybrid vigor. Proceedings of the National Academy of Sciences. 2006;103(35):12957-8. Teng J, Wang D, Zhao C, Zhang X, Chen Z, Liu J, et al. Longitudinal genome-wide association studies of milk production traits in Holstein cattle using whole-genome sequence data imputed from medium-density chip data. Journal of Dairy Science. 2023;106(4):2535–50. Bernini F, Mancin E, Sartori C, Mantovani R, Vevey M, Blanchet V, et al. Genome-wide association studies for milk production traits in two autochthonous Aosta cattle breeds. Animal. 2024;18(10):101322. Reynolds EG, Lopdell T, Wang Y, Tiplady KM, Harland CS, Johnson TJ, et al. Non-additive QTL mapping of lactation traits in 124,000 cattle reveals novel recessive loci. Genetics Selection Evolution. 2022;54(1):5. Bisutti V, Pegolo S, Giannuzzi D, Mota L, Vanzin A, Toscano A, et al. The β-casein (CSN2) A2 allelic variant alters milk protein profile and slightly worsens coagulation properties in Holstein cows. Journal of Dairy Science. 2022;105(5):3794–809. Miluchová M, Gábor M, Candrák J. The effect of the genotypes of the CSN2 gene on test-day milk yields in the Slovak Holstein cow. Agriculture. 2023;13(1):154. Li H, Durbin R. Fast and accurate short read alignment with Burrows–Wheeler transform. Bioinformatics. 2009;25(14):1754–60. Han B, Yuan Y, Li Y, Liu L, Sun D. Single nucleotide polymorphisms of NUCB2 and their genetic associations with milk production traits in dairy cows. Genes. 2019;10(6):449. Bolger AM, Lohse M, Usadel B. Trimmomatic: a flexible trimmer for Illumina sequence data. Bioinformatics. 2014;30(15):2114–20. Shen W, Le S, Li Y, Hu F. SeqKit: a cross-platform and ultrafast toolkit for FASTA/Q file manipulation. PloS one. 2016;11(10):e0163962. Dobin A, Davis CA, Schlesinger F, Drenkow J, Zaleski C, Jha S, et al. STAR: ultrafast universal RNA-seq aligner. Bioinformatics. 2013;29(1):15–21. Kim D, Paggi JM, Park C, Bennett C, Salzberg SL. Graph-based genome alignment and genotyping with HISAT2 and HISAT-genotype. Nature biotechnology. 2019;37(8):907–15. Li H, Handsaker B, Wysoker A, Fennell T, Ruan J, Homer N, et al. The sequence alignment/map format and SAMtools. bioinformatics. 2009;25(16):2078–9. Pertea M, Pertea GM, Antonescu CM, Chang T-C, Mendell JT, Salzberg SL. StringTie enables improved reconstruction of a transcriptome from RNA-seq reads. Nature biotechnology. 2015;33(3):290–5. Galili T. dendextend: an R package for visualizing, adjusting and comparing trees of hierarchical clustering. Bioinformatics. 2015;31(22):3718–20. McKenna A, Hanna M, Banks E, Sivachenko A, Cibulskis K, Kernytsky A, et al. The Genome Analysis Toolkit: a MapReduce framework for analyzing next-generation DNA sequencing data. Genome research. 2010;20(9):1297–303. Wang K, Li M, Hakonarson H. ANNOVAR: functional annotation of genetic variants from high-throughput sequencing data. Nucleic acids research. 2010;38(16):e164-e. Danecek P, Auton A, Abecasis G, Albers CA, Banks E, DePristo MA, et al. The variant call format and VCFtools. Bioinformatics. 2011;27(15):2156–8. Browning SR, Browning BL. Rapid and accurate haplotype phasing and missing-data inference for whole-genome association studies by use of localized haplotype clustering. The American Journal of Human Genetics. 2007;81(5):1084–97. Zhang Z, Wang A, Hu H, Wang L, Gong M, Yang Q, et al. The efficient phasing and imputation pipeline of low-coverage whole genome sequencing data using a high‐quality and publicly available reference panel in cattle. Animal Research One Health. 2023;1(1):4–16. Danecek P, Bonfield JK, Liddle J, Marshall J, Ohan V, Pollard MO, et al. Twelve years of SAMtools and BCFtools. Gigascience. 2021;10(2):giab008. Purcell S, Neale B, Todd-Brown K, Thomas L, Ferreira MA, Bender D, et al. PLINK: a tool set for whole-genome association and population-based linkage analyses. The American journal of human genetics. 2007;81(3):559–75. Tamura K, Stecher G, Kumar S. MEGA11: molecular evolutionary genetics analysis version 11. Molecular biology evolution. 2021;38(7):3022–7. Letunic I, Bork P. Interactive Tree of Life (iTOL) v6: recent updates to the phylogenetic tree display and annotation tool. Nucleic Acids Research. 2024:gkae268. Liu Y, Liu X, Zheng Z, Ma T, Liu Y, Long H, et al. Genome-wide analysis of expression QTL (eQTL) and allele-specific expression (ASE) in pig muscle identifies candidate genes for meat quality traits. Genetics Selection Evolution. 2020;52:1–11. Mohammadi P, Castel SE, Brown AA, Lappalainen T. Quantifying the regulatory effect size of cis-acting genetic variation using allelic fold change. Genome research. 2017;27(11):1872–84. Sherman BT, Hao M, Qiu J, Jiao X, Baseler MW, Lane HC, et al. DAVID: a web server for functional enrichment analysis and functional annotation of gene lists (2021 update). Nucleic acids research. 2022;50(W1):W216-W21. Zhou Y, Zhang Z, Bao Z, Li H, Lyu Y, Zan Y, et al. Graph pangenome captures missing heritability and empowers tomato breeding. Nature. 2022;606(7914):527–34. Stegle O, Parts L, Piipari M, Winn J, Durbin R. Using probabilistic estimation of expression residuals (PEER) to obtain increased power and interpretability of gene expression analyses. Nature protocols. 2012;7(3):500–7. Delaneau O, Ongen H, Brown AA, Fort A, Panousis NI, Dermitzakis ET. A complete tool set for molecular QTL discovery and analysis. Nature communications. 2017;8(1):15452. Li YI, Knowles DA, Humphrey J, Barbeira AN, Dickinson SP, Im HK, et al. Annotation-free quantification of RNA splicing using LeafCutter. Nature genetics. 2018;50(1):151–8. Wu D-D, Ding X-D, Wang S, Wójcik JM, Zhang Y, Tokarska M, et al. Pervasive introgression facilitated domestication and adaptation in the Bos species complex. Nature ecology evolution. 2018;2(7):1139–45. Krueger F. Trim Galore!: A wrapper around Cutadapt and FastQC to consistently apply adapter and quality trimming to FastQ files, with extra functionality for RRBS data. Babraham Institute. 2015; https://github.com/FelixKrueger/TrimGalore . Langmead B, Salzberg SL. Fast gapped-read alignment with Bowtie 2. Nature methods. 2012;9(4):357–9. Tarasov A, Vilella AJ, Cuppen E, Nijman IJ, Prins P. Sambamba: fast processing of NGS alignment formats. Bioinformatics. 2015;31(12):2032–4. Zhang Y, Liu T, Meyer CA, Eeckhoute J, Johnson DS, Bernstein BE, et al. Model-based analysis of ChIP-Seq (MACS). Genome biology. 2008;9:1–9. Quinlan AR, Hall IM. BEDTools: a flexible suite of utilities for comparing genomic features. Bioinformatics. 2010;26(6):841–2. Prowse-Wilkins CP, Wang J, Xiang R, Garner JB, Goddard ME, Chamberlain AJ. Putative causal variants are enriched in annotated functional regions from six bovine tissues. Frontiers in genetics. 2021;12:664379. Prowse-Wilkins CP, Lopdell TJ, Xiang R, Vander Jagt CJ, Littlejohn MD, Chamberlain AJ, et al. Genetic variation in histone modifications and gene expression identifies regulatory variants in the mammary gland of cattle. BMC genomics. 2022;23(1):815. Prowse-Wilkins CP, Wang J, Garner JB, Goddard ME, Chamberlain AJ. Allele specific binding of histone modifications and a transcription factor does not predict allele specific expression in correlated ChIP-seq peak-exon pairs. Scientific Reports. 2023;13(1):15596. Villar D, Berthelot C, Aldridge S, Rayner TF, Lukk M, Pignatelli M, et al. Enhancer evolution across 20 mammalian species. Cell. 2015;160(3):554–66. Additional Declarations No competing interests reported. Supplementary Files supplementaryTable.xlsx Supplementary.docx Cite Share Download PDF Status: Under Review Version 1 posted Editorial decision: Revision requested 12 Feb, 2025 Reviews received at journal 02 Feb, 2025 Reviews received at journal 31 Jan, 2025 Reviewers agreed at journal 24 Jan, 2025 Reviewers agreed at journal 24 Jan, 2025 Reviews received at journal 02 Jan, 2025 Reviewers agreed at journal 12 Dec, 2024 Reviewers invited by journal 11 Dec, 2024 Editor assigned by journal 27 Nov, 2024 Submission checks completed at journal 27 Nov, 2024 First submitted to journal 26 Nov, 2024 You are reading this latest preprint version Research Square lets you share your work early, gain feedback from the community, and start making changes to your manuscript prior to peer review in a journal. As a division of Research Square Company, we’re committed to making research communication faster, fairer, and more useful. We do this by developing innovative software and high quality services for the global research community. Our growing team is made up of researchers and industry professionals working together to solve the most critical problems facing scientific publishing. Also discoverable on Platform About Our Team In Review Editorial Policies Advisory Board Help Center Resources Author Services Accessibility API Access RSS feed Manage Cookie Preferences © Research Square 2026 | ISSN 2693-5015 (online) Privacy Policy Terms of Service Do Not Sell My Personal Information {"props":{"pageProps":{"initialData":{"identity":"rs-5530951","acceptedTermsAndConditions":true,"allowDirectSubmit":false,"archivedVersions":[],"articleType":"Research Article","associatedPublications":[],"authors":[{"id":396691766,"identity":"358b2700-6ac7-45b4-8be7-ae315b79b4b8","order_by":0,"name":"Jiaqi Li","email":"","orcid":"","institution":"Xinjiang Agricultural University","correspondingAuthor":false,"prefix":"","firstName":"Jiaqi","middleName":"","lastName":"Li","suffix":""},{"id":396691767,"identity":"4aab18e7-6b6b-40fd-b26c-a7205e0e8173","order_by":1,"name":"Lei Xu","email":"","orcid":"","institution":"Xinjiang Agricultural University","correspondingAuthor":false,"prefix":"","firstName":"Lei","middleName":"","lastName":"Xu","suffix":""},{"id":396691768,"identity":"1e6c745a-53fb-428c-9b08-f5084a9ff7c2","order_by":2,"name":"Xiaoyun Liang","email":"","orcid":"","institution":"Xinjiang Agricultural University","correspondingAuthor":false,"prefix":"","firstName":"Xiaoyun","middleName":"","lastName":"Liang","suffix":""},{"id":396691769,"identity":"d8fb662d-cb02-4fbe-9a4e-a9ffca827417","order_by":3,"name":"Letian Li","email":"","orcid":"","institution":"Xinjiang Agricultural University","correspondingAuthor":false,"prefix":"","firstName":"Letian","middleName":"","lastName":"Li","suffix":""},{"id":396691770,"identity":"2becf4a7-50f0-42c6-a21b-93d53373e92b","order_by":4,"name":"Xixia Huang","email":"","orcid":"","institution":"Xinjiang Agricultural University","correspondingAuthor":false,"prefix":"","firstName":"Xixia","middleName":"","lastName":"Huang","suffix":""},{"id":396691771,"identity":"bfc0ab45-6dc5-4683-a835-66da9b9888ef","order_by":5,"name":"Qiuming Chen","email":"data:image/png;base64,iVBORw0KGgoAAAANSUhEUgAAAZAAAAAyAQMAAABI0h/eAAAABlBMVEX///8AAABVwtN+AAAACXBIWXMAAA7EAAAOxAGVKw4bAAAA3UlEQVRIiWNgGAWjYJCCAx8qbHj4mZkPPkiosCFKB+PBGWfS5CTb2ZINHpxJI0oL82HetkPGBv08ZpIP2w4RVs/fnp1wmOfMgcQNzAxmFQlsB4Ai3Ql4tUicebvh4JyKO4nbmRnSbiTw3AGKnN2AV4uBRO6GA2/OPEvc2cxw7EaCxDOwCGEtvG2HEzccZmwrSDA4TJyWg0AtxgaHmdkYEhKI0AL2CziQm9mYJRIOpPEQ9At/e+7mD+Co5D//8ePPfzZy/O29+LUwMCSgcnkIKMeiZRSMglEwCkYBBgAAeKtVW1Zu2kEAAAAASUVORK5CYII=","orcid":"","institution":"Xinjiang Agricultural University","correspondingAuthor":true,"prefix":"","firstName":"Qiuming","middleName":"","lastName":"Chen","suffix":""}],"badges":[],"createdAt":"2024-11-26 23:53:08","currentVersionCode":1,"declarations":"","doi":"10.21203/rs.3.rs-5530951/v1","doiUrl":"https://doi.org/10.21203/rs.3.rs-5530951/v1","draftVersion":[],"editorialEvents":[],"editorialNote":"","failedWorkflow":false,"files":[{"id":73357321,"identity":"6a4585d3-7ba6-4aa8-9f81-8a44faac04a4","added_by":"auto","created_at":"2025-01-09 08:17:51","extension":"png","order_by":1,"title":"Figure 1","display":"","copyAsset":false,"role":"figure","size":5655929,"visible":true,"origin":"","legend":"\u003cp\u003eComprehensive identification of ASE variants. (\u003cstrong\u003eA\u003c/strong\u003e) Schematic representation of the bioinformatics pipeline employed for the identification of ASE variants across samples. (\u003cstrong\u003eB\u003c/strong\u003e) Overview of the 35 tissue types analyzed in the study, with the corresponding number of samples and identified ASE variants shown in parentheses.\u003c/p\u003e","description":"","filename":"figure1.png","url":"https://assets-eu.researchsquare.com/files/rs-5530951/v1/ff787ab996f8fad1c9699cb5.png"},{"id":73358777,"identity":"fbc7ce9a-ee67-47b6-8789-4b6fa3033656","added_by":"auto","created_at":"2025-01-09 08:25:51","extension":"png","order_by":2,"title":"Figure 2","display":"","copyAsset":false,"role":"figure","size":3650490,"visible":true,"origin":"","legend":"\u003cp\u003eGeneral characteristics of ASE variants. (\u003cstrong\u003eA\u003c/strong\u003e) Fold enrichment analysis of ASE variants across different sequence ontology categories. (\u003cstrong\u003eB\u003c/strong\u003e) Density plot showing the distribution of distances from ASE variants to the TSS of the nearest gene. (\u003cstrong\u003eC\u003c/strong\u003e) Curve depicting the distribution of shared ASE variants across different tissues, illustrating the balance between tissue-specific and common ASE variants. (\u003cstrong\u003eD\u003c/strong\u003e) Scatter plot depicting the relationship between the number of identified ASE variants and sample size within each dataset. (\u003cstrong\u003eE\u003c/strong\u003e) Relationship between the number of ASE variants and heterozygosity for each dataset. (\u003cstrong\u003eF\u003c/strong\u003e) Replication rate analysis of ASE variants within the same tissue across different BioProjects, providing insights into the reproducibility of ASE events. (\u003cstrong\u003eG\u003c/strong\u003e) Density plot illustrating the distribution of effect sizes (allelic fold change, aFC) among ASE variants. (\u003cstrong\u003eH\u003c/strong\u003e) Bar chart showing the number of ASE variants exhibiting opposite directionality across tissues. (\u003cstrong\u003eI\u003c/strong\u003e) Scatter plot analyzing the relationship between the number of ASE variants with opposite directionality and the variance explained by the first PEER factor. (\u003cstrong\u003eJ\u003c/strong\u003e) Box plot comparing the absolute effect sizes of ASE variants with consistent versus opposite directionality. (\u003cstrong\u003eK\u003c/strong\u003e) Distribution of the number of independent ASE variants per gene (r\u003csup\u003e2\u003c/sup\u003e \u0026gt; 0.6).\u003c/p\u003e","description":"","filename":"figure2.png","url":"https://assets-eu.researchsquare.com/files/rs-5530951/v1/b3923c34fea16bfee6f5a84b.png"},{"id":73357346,"identity":"09331936-8274-4847-84d9-0ccd29d3873d","added_by":"auto","created_at":"2025-01-09 08:17:52","extension":"png","order_by":3,"title":"Figure 3","display":"","copyAsset":false,"role":"figure","size":8705290,"visible":true,"origin":"","legend":"\u003cp\u003eDominant effects show closer proximity to ASE events compared to additive effects. (\u003cstrong\u003eA\u003c/strong\u003e) Fold enrichment analysis comparing the distribution of additive and dominant cis-eQTLs and cis-sQTLs across ASE variants. Scatter plot showing the relationship between the effect size (aFC) of shared ASE variants and additive \u003cem\u003ecis\u003c/em\u003e-eQTLs (\u003cstrong\u003eB\u003c/strong\u003e), dominant \u003cem\u003ecis\u003c/em\u003e-eQTLs (\u003cstrong\u003eC\u003c/strong\u003e), degree of dominance (t\u003csub\u003eDom\u003c/sub\u003e/t\u003csub\u003eAdd\u003c/sub\u003e) in \u003cem\u003ecis\u003c/em\u003e-eQTLs (\u003cstrong\u003eD\u003c/strong\u003e), additive \u003cem\u003ecis\u003c/em\u003e-sQTLs (\u003cstrong\u003eE\u003c/strong\u003e), dominant \u003cem\u003ecis\u003c/em\u003e-sQTLs (\u003cstrong\u003eF\u003c/strong\u003e), and degree of dominance (t\u003csub\u003eDom\u003c/sub\u003e/t\u003csub\u003eAdd\u003c/sub\u003e) in \u003cem\u003ecis\u003c/em\u003e-sQTLs (\u003cstrong\u003eG\u003c/strong\u003e).\u003c/p\u003e","description":"","filename":"Figure3.png","url":"https://assets-eu.researchsquare.com/files/rs-5530951/v1/c11d4ef56551b25059a53b84.png"},{"id":73357329,"identity":"ccd2275b-b9d4-4b30-afde-52fdee50c3d0","added_by":"auto","created_at":"2025-01-09 08:17:52","extension":"png","order_by":4,"title":"Figure 4","display":"","copyAsset":false,"role":"figure","size":2771008,"visible":true,"origin":"","legend":"\u003cp\u003eFunctional annotation of ASE variants. (\u003cstrong\u003eA\u003c/strong\u003e) Fold enrichment analysis illustrating the conservation scores of ASE variants. (\u003cstrong\u003eB\u003c/strong\u003e) Proportion of ASE variants located within ATAC-seq peaks. (\u003cstrong\u003eC\u003c/strong\u003e) Proportion of ASE variants located within ChIP-seq peaks. (\u003cstrong\u003eD\u003c/strong\u003e) Fold enrichment analysis of ASE variants within cattle QTL regions from the cattle QTL database.\u003c/p\u003e","description":"","filename":"Figure4.png","url":"https://assets-eu.researchsquare.com/files/rs-5530951/v1/f2c539dc19533bc2f1b72baf.png"},{"id":73358779,"identity":"460cd079-75dc-4c87-8d08-255c63597c0d","added_by":"auto","created_at":"2025-01-09 08:25:52","extension":"png","order_by":5,"title":"Figure 5","display":"","copyAsset":false,"role":"figure","size":2612374,"visible":true,"origin":"","legend":"\u003cp\u003eValidation of shared variants between ASE and QTLs for milk-related traits. (\u003cstrong\u003eA\u003c/strong\u003e) Pipeline for validating ASE variants using multiplex PCR. Manhattan plot for nine milk-related traits using a mixed linear model of additive effects (\u003cstrong\u003eB\u003c/strong\u003e) and dominant effects (\u003cstrong\u003eC\u003c/strong\u003e). The red and black lines indicate the significant and suggestive thresholds, respectively. No significant or suggestive signatures were observed for three traits. Scatter plot showing the relationship between the effect size (aFC) of shared ASE variants and the degree of dominance (t\u003csub\u003eDom\u003c/sub\u003e/t\u003csub\u003eAdd\u003c/sub\u003e) for variants at the significant threshold (\u003cstrong\u003eD\u003c/strong\u003e) and suggestive threshold (\u003cstrong\u003eE\u003c/strong\u003e).\u003c/p\u003e","description":"","filename":"Figure5.png","url":"https://assets-eu.researchsquare.com/files/rs-5530951/v1/ddf8ae69f13b07098f9a8bb1.png"},{"id":73360417,"identity":"95543c9b-0928-45cc-a5e3-80c269200ed0","added_by":"auto","created_at":"2025-01-09 08:42:06","extension":"pdf","order_by":0,"title":"","display":"","copyAsset":false,"role":"manuscript-pdf","size":22192614,"visible":true,"origin":"","legend":"","description":"","filename":"manuscript.pdf","url":"https://assets-eu.researchsquare.com/files/rs-5530951/v1/4b8a5f45-2589-45b2-b458-205d3c5f24cb.pdf"},{"id":73357323,"identity":"52a184b2-7805-459a-8c28-1c84c8e1c59f","added_by":"auto","created_at":"2025-01-09 08:17:51","extension":"xlsx","order_by":0,"title":"","display":"","copyAsset":false,"role":"supplement","size":1003185,"visible":true,"origin":"","legend":"","description":"","filename":"supplementaryTable.xlsx","url":"https://assets-eu.researchsquare.com/files/rs-5530951/v1/3deac9dd7ebe8d2203981e18.xlsx"},{"id":73358780,"identity":"644a31e3-3c39-4d95-a69d-88da40b0e651","added_by":"auto","created_at":"2025-01-09 08:25:52","extension":"docx","order_by":1,"title":"","display":"","copyAsset":false,"role":"supplement","size":9877029,"visible":true,"origin":"","legend":"","description":"","filename":"Supplementary.docx","url":"https://assets-eu.researchsquare.com/files/rs-5530951/v1/09522747f0cac7653138d391.docx"}],"financialInterests":"No competing interests reported.","formattedTitle":"A multi-tissue atlas of allelic-specific expression reveals the characteristics, mechanisms, and relationship with dominant effects in cattle","fulltext":[{"header":"Introduction","content":"\u003cp\u003eOver the past decade, genome-wide association studies (GWAS) have consistently demonstrated that most variants associated with diseases and other complex traits are located in non-coding regions of the genome [\u003cspan citationid=\"CR1\" class=\"CitationRef\"\u003e1\u003c/span\u003e, \u003cspan citationid=\"CR2\" class=\"CitationRef\"\u003e2\u003c/span\u003e]. Understanding how these genetic variants influence gene expression is crucial for uncovering the genetic basis of complex traits. To address this, several Genotype-Tissue Expression (GTEx) projects have been initiated across multiple species [\u003cspan additionalcitationids=\"CR4\" citationid=\"CR3\" class=\"CitationRef\"\u003e3\u003c/span\u003e\u0026ndash;\u003cspan citationid=\"CR5\" class=\"CitationRef\"\u003e5\u003c/span\u003e]. These projects aim to identify expression quantitative trait loci (eQTLs), which are genetic variants that affect gene expression levels. Depending on their proximity to the target genes, eQTLs can be classified into \u003cem\u003ecis\u003c/em\u003e-eQTLs, which act on nearby genes, and \u003cem\u003etrans\u003c/em\u003e-eQTLs, which influence distant genes, sometimes on different chromosomes.\u003c/p\u003e \u003cp\u003eIn these GTEx initiatives, allelic-specific expression (ASE) has been employed as a robust tool for validating \u003cem\u003ecis\u003c/em\u003e-eQTL [\u003cspan citationid=\"CR4\" class=\"CitationRef\"\u003e4\u003c/span\u003e, \u003cspan citationid=\"CR6\" class=\"CitationRef\"\u003e6\u003c/span\u003e, \u003cspan citationid=\"CR7\" class=\"CitationRef\"\u003e7\u003c/span\u003e]. ASE measures deviations from the expected equal expression of two alleles at a heterozygous locus within an individual. The internal consistency in environmental and technical factors between alleles within the same individual makes ASE a powerful approach. By integrating ASE with other approaches, such as eQTL mapping [\u003cspan citationid=\"CR8\" class=\"CitationRef\"\u003e8\u003c/span\u003e], GWAS [\u003cspan citationid=\"CR9\" class=\"CitationRef\"\u003e9\u003c/span\u003e], and selective sweep analysis [\u003cspan citationid=\"CR10\" class=\"CitationRef\"\u003e10\u003c/span\u003e], researchers have identified a smaller yet more plausible set of causal variants that affect complex traits [\u003cspan citationid=\"CR11\" class=\"CitationRef\"\u003e11\u003c/span\u003e]. However, comprehensive ASE analyses across multiple tissues have largely been limited to studies in humans [\u003cspan citationid=\"CR3\" class=\"CitationRef\"\u003e3\u003c/span\u003e, \u003cspan citationid=\"CR12\" class=\"CitationRef\"\u003e12\u003c/span\u003e].\u003c/p\u003e \u003cp\u003eWhile ASE can identify cis-regulatory elements [\u003cspan citationid=\"CR13\" class=\"CitationRef\"\u003e13\u003c/span\u003e], epigenomics remains the primary approach for uncovering these elements. This method has been extensively utilized in projects, such as the Encyclopedia of DNA Elements (ENCODE), which aims to identify causative genetic variants related to human diseases [\u003cspan citationid=\"CR14\" class=\"CitationRef\"\u003e14\u003c/span\u003e], and the Functional Annotation of Animal Genomes (FAANG), which seeks to improve the quantity and quality of agricultural products through molecular selection [\u003cspan citationid=\"CR15\" class=\"CitationRef\"\u003e15\u003c/span\u003e]. Some studies have started integrating ASE and epigenomic data to mutually validate findings [\u003cspan citationid=\"CR16\" class=\"CitationRef\"\u003e16\u003c/span\u003e] or to identify causal variants [\u003cspan citationid=\"CR17\" class=\"CitationRef\"\u003e17\u003c/span\u003e], but the relationship between ASE and epigenomic regulatory elements remains poorly understood.\u003c/p\u003e \u003cp\u003eASE also holds promise in unraveling the genetic basis of heterosis (hybrid vigor) [\u003cspan citationid=\"CR18\" class=\"CitationRef\"\u003e18\u003c/span\u003e, \u003cspan citationid=\"CR19\" class=\"CitationRef\"\u003e19\u003c/span\u003e]. The mechanisms behind heterosis have long been debated, primarily through two genetic models: dominance and overdominance, which are based on the degree of dominance of genetic variants, defined by the ratio of dominant to additive effect [\u003cspan citationid=\"CR20\" class=\"CitationRef\"\u003e20\u003c/span\u003e]. Although the relationship between dominant effects and ASE events has been discussed [\u003cspan citationid=\"CR18\" class=\"CitationRef\"\u003e18\u003c/span\u003e, \u003cspan citationid=\"CR19\" class=\"CitationRef\"\u003e19\u003c/span\u003e], comprehensive and quantitative analyses are still needed.\u003c/p\u003e \u003cp\u003eIn this study, we present a novel pipeline that uniformly integrates 7,532 public RNA-seq datasets to identify ASE variants in cattle. Using these variants, we performed three key analyses: (1) characterization of ASE events, including tissue specificity and directionality; (2) exploration of the relationship between ASE and dominant cis-eQTLs linked to molecular phenotypes of gene expression and dominant QTLs linked to physiological phenotypes of milk-related traits in Holstein cattle; and (3) functional annotation of ASE variants using epigenomic data, QTL databases, and conservation scores of sequences. Our findings not only advance the understanding of the genetic mechanisms underlying ASE events and heterosis but also provide valuable genetic resources for genomics and breeding programs in cattle.\u003c/p\u003e"},{"header":"Results","content":"\u003cdiv id=\"Sec3\" class=\"Section2\"\u003e \u003ch2\u003eASE discovery\u003c/h2\u003e \u003cp\u003eWe downloaded 8,518 public RNA-seq samples from the NCBI SRA database, generating approximately 363\u0026nbsp;billion clean reads. After data cleaning (Materials and Methods), we retained 7,532 samples across 157 subsets, which spanned 35 distinct tissues and 126 NCBI BioProjects (Fig.\u0026nbsp;\u003cspan refid=\"Fig1\" class=\"InternalRef\"\u003e1\u003c/span\u003eA and Additional file 1: Table \u003cspan refid=\"MOESM1\" class=\"InternalRef\"\u003eS1\u003c/span\u003e). Among these subsets, 57 were well-defined as single-breed samples (Additional file 2: Fig. \u003cspan refid=\"MOESM1\" class=\"InternalRef\"\u003eS1\u003c/span\u003e). The median sample size across these tissues was 80, with a range from 20 in ovary to 1,688 in whole blood. Based on the transcripts per million (TPM) values for 25,365 expressed genes (mean TPM\u0026thinsp;\u0026gt;\u0026thinsp;0.1), hierarchical cluster analysis recapitulated tissue types. Unlike previous human GTEx studies [\u003cspan citationid=\"CR3\" class=\"CitationRef\"\u003e3\u003c/span\u003e], which identified blood-related samples as the primary outgroup, our research revealed the inclusion of mammary gland and milk somatic cell within the outgroup, underscoring the significance and similarity between these two tissues (Additional file 2: Fig. \u003cspan refid=\"MOESM2\" class=\"InternalRef\"\u003eS2\u003c/span\u003e). A neighbor-joining (NJ) tree of imputed genotype primarily reflected two ancestries (\u003cem\u003eBos taurus\u003c/em\u003e and \u003cem\u003eBos indicus\u003c/em\u003e) (Additional file 2: Fig. S3). The observed separations of tissue type in gene expression and the distinctions of ancestry in genotype are consistent with findings from the Cattle Genotype-Tissue Expression (CattleGTEx) atlas [\u003cspan citationid=\"CR4\" class=\"CitationRef\"\u003e4\u003c/span\u003e], suggesting the high quality and utility of our data processing for follow-up analyses.\u003c/p\u003e \u003cp\u003e \u003c/p\u003e \u003cp\u003eTo mitigate heterogeneity across different BioProjects, our ASE analysis focused on subsets with at least 20 individuals per tissue within each BioProject. We used the GATK ASEReadCounter module to calculate allele counts at heterozygous sites for each individual. Subsequently, we applied stringent filtering based on read counts and allele ratios, conducted a binomial test with false discovery rate (FDR) adjustment at the individual level, and further filtered by the number of heterozygous sites and the ratio of significant allelic imbalance at the population level. Through these processes, we identified a total of 452,052 ASE variants. After removing redundancies across subsets, we identified 161,059 unique ASE variants (Fig.\u0026nbsp;\u003cspan refid=\"Fig1\" class=\"InternalRef\"\u003e1\u003c/span\u003e). These variants were annotated to 13,136 genes, accounting for 66.1% of all protein-coding genes, which is lower than the 94.7% reported in human eQTL studies [\u003cspan citationid=\"CR21\" class=\"CitationRef\"\u003e21\u003c/span\u003e]. Genes lacking ASE variants were enriched for those lacking expression in our analyzed tissues, such as those involved in keratin filament and acrosomal vesicle (Additional file 2: Table \u003cspan refid=\"MOESM2\" class=\"InternalRef\"\u003eS2\u003c/span\u003e).\u003c/p\u003e \u003cp\u003eASE variants were primarily distributed across exons, introns, 3' untranslated regions (UTRs), and intergenic regions, collectively accounting for 92.33% of the total, with 45,835 in exons, 39,590 in introns, 35,991 in 3' UTRs, and 27,302 in intergenic regions (Additional file 2: Table S3). The distribution pattern of ASE variants shows a similar trend to that of eQTLs in humans [\u003cspan citationid=\"CR21\" class=\"CitationRef\"\u003e21\u003c/span\u003e] and pigs [\u003cspan citationid=\"CR5\" class=\"CitationRef\"\u003e5\u003c/span\u003e], with protein-coding regions (e.g., stop gain, stop loss, nonsynonymous and synonymous mutation) showing the highest enrichment (Fig.\u0026nbsp;\u003cspan refid=\"Fig2\" class=\"InternalRef\"\u003e2\u003c/span\u003eA). This pattern is further supported by the observation that more than 90% of the ASE variants were located downstream of the transcriptional start site (TSS) of the nearest gene (Fig.\u0026nbsp;\u003cspan refid=\"Fig2\" class=\"InternalRef\"\u003e2\u003c/span\u003eB), reflecting the transcriptomic nature of ASE detection. Furthermore, we found a higher enrichment of variants affecting mRNA splicing, consistent with previous eQTL studies in humans [\u003cspan citationid=\"CR21\" class=\"CitationRef\"\u003e21\u003c/span\u003e] and pigs [\u003cspan citationid=\"CR5\" class=\"CitationRef\"\u003e5\u003c/span\u003e]. Additionally, we observed a more significant enrichment in the downstream compared to upstream regions and in the 3\u0026prime; UTR compared to the 5\u0026prime; UTR of protein-coding genes (Fig.\u0026nbsp;\u003cspan refid=\"Fig2\" class=\"InternalRef\"\u003e2\u003c/span\u003eA). This contrasts with findings from humans [\u003cspan citationid=\"CR3\" class=\"CitationRef\"\u003e3\u003c/span\u003e] and cattle [\u003cspan citationid=\"CR4\" class=\"CitationRef\"\u003e4\u003c/span\u003e] eQTL studies, underscoring the distinct characteristics of ASE events.\u003c/p\u003e \u003cp\u003e \u003c/p\u003e \u003cp\u003eAmong the identified ASE variants, 72.8% were unique to a single tissue, 14.2% were found in two tissues, and 5.3% appeared in three tissues (Fig.\u0026nbsp;\u003cspan refid=\"Fig2\" class=\"InternalRef\"\u003e2\u003c/span\u003eC). To reduce heterogeneity across different BioProjects, we compared ASE variants within each BioProject, revealing that on average, 68.4% of ASE variants (ranging from 8.15\u0026ndash;98.7%) were tissue-specific (Additional file 2: Fig. S4 and S5A). Regarding the shared expression of heterozygosity across tissues, an average of 40.0% of ASE variants (ranging from 6.3\u0026ndash;71.3%) are exclusive to one tissue (Additional file 2: Fig. S5B), while an average of 53.3% of ASE variants (ranging from 22.6\u0026ndash;87.9%) were replicated across different tissues (Additional file 2: Fig. S6). The distribution curve of shared ASE variants among tissues exhibited an inverse S-shaped pattern (Fig.\u0026nbsp;\u003cspan refid=\"Fig2\" class=\"InternalRef\"\u003e2\u003c/span\u003eC), contrasting with the U-shaped pattern typically observed in eQTL studies in humans [\u003cspan citationid=\"CR3\" class=\"CitationRef\"\u003e3\u003c/span\u003e], cattle [\u003cspan citationid=\"CR4\" class=\"CitationRef\"\u003e4\u003c/span\u003e], and pigs [\u003cspan citationid=\"CR5\" class=\"CitationRef\"\u003e5\u003c/span\u003e]. These results suggest that ASE events exhibit a higher degree of tissue specificity. Additionally, we calculated the Pearson correlation coefficient of effect size for shared ASE variants across different tissues, yielding a median value of 0.7549 (ranging from \u0026minus;\u0026thinsp;0.6731 to 0.9941) for the 7,585 significant correlations out of 8,342 pairwise comparisons (Additional file 2: Fig. S7). This indicates that most shared ASE variants across tissues are likely driven by the same underlying factors.\u003c/p\u003e \u003cp\u003eThe discovery of ASE variants showed a significant correlation with sample size (Pearson r\u0026thinsp;=\u0026thinsp;0.55; two-sided Student\u0026rsquo;s t-test: \u003cem\u003eP\u003c/em\u003e\u0026thinsp;=\u0026thinsp;6.61 \u0026times; 10\u003csup\u003e\u0026minus;\u0026thinsp;14\u003c/sup\u003e) (Fig.\u0026nbsp;\u003cspan refid=\"Fig2\" class=\"InternalRef\"\u003e2\u003c/span\u003eD), a correlation lower than those reported for eGenes (0.85) and sGenes (0.63) in cattle [\u003cspan citationid=\"CR4\" class=\"CitationRef\"\u003e4\u003c/span\u003e], underscoring the unique nature of ASE events. A similar correlation was observed between the number of ASE variants and heterozygosity (Pearson r\u0026thinsp;=\u0026thinsp;0.54; two-sided Student\u0026rsquo;s t-test: \u003cem\u003eP\u003c/em\u003e\u0026thinsp;=\u0026thinsp;2.09 \u0026times; 10\u003csup\u003e\u0026minus;\u0026thinsp;13\u003c/sup\u003e) (Fig.\u0026nbsp;\u003cspan refid=\"Fig2\" class=\"InternalRef\"\u003e2\u003c/span\u003eE). However, substantial variability in the number of ASE variants was observed among immune-related tissues, further emphasizing the tissue specificity of ASE events.\u003c/p\u003e \u003cp\u003eThe average replication rate of ASE variants between different BioProjects of the same tissue was 49.1%, ranging from 1.9\u0026ndash;100%. Notably, the replication rate was lower in reproductive tissues, suggesting a higher complexity in these tissues (Fig.\u0026nbsp;\u003cspan refid=\"Fig2\" class=\"InternalRef\"\u003e2\u003c/span\u003eF). Despite the lower replication rate, the median Pearson correlation coefficient of effect size for shared ASE variants was relatively high (0.85, ranging from 0.21 to 0.99) in the 707 significant correlations out of 740 pairwise comparisons of the same tissue between different BioProjects (Additional file 2: Fig. S8), indicating that the majority of detected ASE variants were likely genuine.\u003c/p\u003e \u003cp\u003eIn our study, 80% of ASE variants exhibited a more than twofold effect on gene expression (Fig.\u0026nbsp;\u003cspan refid=\"Fig2\" class=\"InternalRef\"\u003e2\u003c/span\u003eG), which is higher than the 22% observed in human \u003cem\u003ecis\u003c/em\u003e-eQTL studies [\u003cspan citationid=\"CR21\" class=\"CitationRef\"\u003e21\u003c/span\u003e]. This underscores the importance and directness of ASE events in identifying regulatory elements. Interestingly, ASE variants in tissues related to reproduction and immune systems exhibited higher effect sizes (Additional file 2: Fig. S9), suggesting a potential heterozygous advantage in reproductive and immune traits.\u003c/p\u003e \u003cp\u003eFurthermore, we identified 23,061 ASE variants with opposite directionality, accounting for 14.3% of the total. Among the 35 tissues examined, immune-related tissues, including milk somatic cells, whole blood, and white cells, exhibited a higher prevalence of ASE variants with opposite directionality (Fig.\u0026nbsp;\u003cspan refid=\"Fig2\" class=\"InternalRef\"\u003e2\u003c/span\u003eH). Notably, the variants rs208685250 of MHC Class I \u003cem\u003eJSP.1\u003c/em\u003e and rs136860823 of MHC Class I \u003cem\u003eBOLA\u003c/em\u003e demonstrated the highest number of subsets with opposite ASE directionality (Additional file 2: Fig. S10). This phenomenon of ASE variants with opposite directionality within the MHC family has also been documented in humans [\u003cspan citationid=\"CR22\" class=\"CitationRef\"\u003e22\u003c/span\u003e]. Gene ontology (GO) enrichment analysis of genes with ASE variants exhibiting opposite directionality highlighted innate immune response as the most significantly enriched biological process (Additional file 2: Table S4). Additionally, we observed a Pearson correlation coefficient of 0.86 between the number of ASE variants with opposite directionality and the variance explained by the first PEER factor (two-sided Student\u0026rsquo;s t-test: \u003cem\u003eP\u003c/em\u003e\u0026thinsp;=\u0026thinsp;3.03 \u0026times; 10\u003csup\u003e\u0026minus;\u0026thinsp;6\u003c/sup\u003e) (Fig.\u0026nbsp;\u003cspan refid=\"Fig2\" class=\"InternalRef\"\u003e2\u003c/span\u003eI), indicating a strong influence of data heterogeneity on ASE directionality. Moreover, the absolute effect size of ASE variants with opposite directionality was significantly lower than that of variants with consistent directionality (two-sided Student\u0026rsquo;s t-test: \u003cem\u003eP\u003c/em\u003e\u0026thinsp;\u0026lt;\u0026thinsp;2.2 \u0026times; 10\u003csup\u003e\u0026minus;\u0026thinsp;16\u003c/sup\u003e) (Fig.\u0026nbsp;\u003cspan refid=\"Fig2\" class=\"InternalRef\"\u003e2\u003c/span\u003eJ). These results suggest that ASE variants with opposite directionality, particularly in immune-related genes, are more resilient.\u003c/p\u003e \u003cp\u003eWe also calculated linkage disequilibrium to examine allelic heterogeneity in gene expression associated with ASE events. Our analysis revealed that 10% of ASE genes contained more than one independent ASE variant (r\u003csup\u003e2\u003c/sup\u003e\u0026thinsp;\u0026gt;\u0026thinsp;0.6) in subsets with sample sizes exceeding 80 (Fig.\u0026nbsp;\u003cspan refid=\"Fig2\" class=\"InternalRef\"\u003e2\u003c/span\u003eK). This proportion is lower than the 46% reported for eGenes [\u003cspan citationid=\"CR4\" class=\"CitationRef\"\u003e4\u003c/span\u003e], suggesting that ASE variants may involve less complex interactions compared to eQTLs.\u003c/p\u003e \u003c/div\u003e\n\u003ch3\u003eRelationship between ASE variants and eQTL\u003c/h3\u003e\n\u003cp\u003eTo better understand the relationship between ASE and eQTL, we performed an extensive analysis of \u003cem\u003ecis\u003c/em\u003e-eQTL and \u003cem\u003ecis\u003c/em\u003e-sQTL using datasets comprising over 80 individuals. This analysis, based on imputed genotype data from RNA-seq, identified an average of 79,814 additive \u003cem\u003ecis\u003c/em\u003e-eQTL variants (ranging from 164 to 315,589), 29,958 dominant \u003cem\u003ecis\u003c/em\u003e-eQTL variants (ranging from 342 to 116,913), 112,686 additive \u003cem\u003ecis\u003c/em\u003e-sQTL variants (ranging from 974 to 375,496), and 50,195 dominant \u003cem\u003ecis\u003c/em\u003e-sQTL variants (ranging from 378 to 220,710) (Additional file 2: Table S5).\u003c/p\u003e \u003cp\u003eAmong the identified ASE variants, we found that, on average, 282 were shared with additive \u003cem\u003ecis\u003c/em\u003e-sQTL variants, 156 with dominant \u003cem\u003ecis\u003c/em\u003e-eQTL variants, 410 with additive \u003cem\u003ecis\u003c/em\u003e-sQTL variants, and 265 with dominant \u003cem\u003ecis\u003c/em\u003e-sQTL variants. These correspond to average fold enrichments of 0.4136 (ranging from 0 to 1.0321), 0.6306 (ranging from 0 to 2.8491), 0.4357 (ranging from 0 to 1.2773), and 0.6089 (ranging from 0 to 1.6658). Although most of these enrichments were not statistically significant, the consistently higher enrichment observed in dominant \u003cem\u003ecis\u003c/em\u003e-eQTLs and \u003cem\u003ecis\u003c/em\u003e-sQTLs compared to their additive counterparts remains noteworthy. Specifically, in 13 out of the 17 subsets where ASE variants overlapped with dominant \u003cem\u003ecis\u003c/em\u003e-eQTL variants, the enrichment of ASE variants was higher within dominant \u003cem\u003ecis\u003c/em\u003e-eQTLs compared to their additive counterparts. Similarly, in 15 out of the 18 subsets where ASE variants overlapped with dominant \u003cem\u003ecis\u003c/em\u003e-sQTL variants, the enrichment was greater in dominant \u003cem\u003ecis\u003c/em\u003e-sQTLs relative to additive \u003cem\u003ecis\u003c/em\u003e-sQTLs. These findings underscore the critical role of dominance effects in ASE regulation and highlight the importance of considering these effects when investigating the genetic mechanisms underlying ASE (Fig.\u0026nbsp;\u003cspan refid=\"Fig3\" class=\"InternalRef\"\u003e3\u003c/span\u003eA and Additional file 2: Table S5).\u003c/p\u003e \u003cp\u003e \u003c/p\u003e \u003cp\u003eTo further investigate the link between ASE and dominance, we calculated the Pearson correlation coefficient between the effect sizes of shared ASE and eQTL variants. The correlation between ASE and dominant \u003cem\u003ecis\u003c/em\u003e-eQTL variants (Pearson r\u0026thinsp;=\u0026thinsp;0.5121; two-sided Student\u0026rsquo;s t-test: \u003cem\u003eP\u003c/em\u003e\u0026thinsp;\u0026lt;\u0026thinsp;2.2 \u0026times; 10\u003csup\u003e\u0026minus;\u0026thinsp;16\u003c/sup\u003e) was notably higher than that between ASE and additive \u003cem\u003ecis\u003c/em\u003e-eQTL variants (Pearson r\u0026thinsp;=\u0026thinsp;0.3752; two-sided Student\u0026rsquo;s t-test: \u003cem\u003eP\u003c/em\u003e\u0026thinsp;\u0026lt;\u0026thinsp;2.2 \u0026times; 10\u003csup\u003e\u0026minus;\u0026thinsp;16\u003c/sup\u003e) (Fig.\u0026nbsp;\u003cspan refid=\"Fig3\" class=\"InternalRef\"\u003e3\u003c/span\u003eB and \u003cspan refid=\"Fig3\" class=\"InternalRef\"\u003e3\u003c/span\u003eC). We also assessed the degree of dominance by computing the ratio of T-statistics (\u003cspan class=\"InlineEquation\"\u003e\u003cspan class=\"mathinline\"\u003e\\(\\:{t}_{Dom}/{t}_{Add}\\)\u003c/span\u003e\u003c/span\u003e) for each eQTL, where \u003cspan class=\"InlineEquation\"\u003e\u003cspan class=\"mathinline\"\u003e\\(\\:{t}_{Dom}\\)\u003c/span\u003e\u003c/span\u003e and \u003cspan class=\"InlineEquation\"\u003e\u003cspan class=\"mathinline\"\u003e\\(\\:{t}_{Add}\\)\u003c/span\u003e\u003c/span\u003e represent the dominant and additive effects, respectively. Interestingly, the correlation between the effect sizes of shared ASE variants and the \u003cspan class=\"InlineEquation\"\u003e\u003cspan class=\"mathinline\"\u003e\\(\\:{t}_{Dom}/{t}_{Add}\\)\u003c/span\u003e\u003c/span\u003e ratio (Pearson r\u0026thinsp;=\u0026thinsp;0.5210; two-sided Student\u0026rsquo;s t-test: \u003cem\u003eP\u003c/em\u003e\u0026thinsp;\u0026lt;\u0026thinsp;2.2 \u0026times; 10\u003csup\u003e\u0026minus;\u0026thinsp;16\u003c/sup\u003e) (Fig.\u0026nbsp;\u003cspan refid=\"Fig3\" class=\"InternalRef\"\u003e3\u003c/span\u003eD) was even higher than the correlation between shared ASE and either dominant or additive \u003cem\u003ecis\u003c/em\u003e-eQTL variants individually. Additionally, in five out of six subsets with significant correlations of ASE effect sizes, including two subsets with a single breed, the \u003cspan class=\"InlineEquation\"\u003e\u003cspan class=\"mathinline\"\u003e\\(\\:{t}_{Dom}/{t}_{Add}\\)\u003c/span\u003e\u003c/span\u003e ratio was higher than the correlations with either additive or dominant \u003cem\u003ecis\u003c/em\u003e-eQTL variants alone (Additional file 2: Fig. S11).\u003c/p\u003e \u003cp\u003eA similar pattern emerged when examining the relationship between shared ASE and \u003cem\u003ecis\u003c/em\u003e-sQTL variants. The correlation between the effect sizes of shared ASE variants and the \u003cspan class=\"InlineEquation\"\u003e\u003cspan class=\"mathinline\"\u003e\\(\\:{t}_{Dom}/{t}_{Add}\\)\u003c/span\u003e\u003c/span\u003e ratio (Pearson r\u0026thinsp;=\u0026thinsp;0.9291; two-sided Student\u0026rsquo;s t-test: \u003cem\u003eP\u003c/em\u003e\u0026thinsp;\u0026lt;\u0026thinsp;2.2 \u0026times; 10\u003csup\u003e\u0026minus;\u0026thinsp;16\u003c/sup\u003e) was higher than that between ASE and dominant \u003cem\u003ecis\u003c/em\u003e-sQTL variants (Pearson r\u0026thinsp;=\u0026thinsp;0.7227; two-sided Student\u0026rsquo;s t-test: \u003cem\u003eP\u003c/em\u003e\u0026thinsp;\u0026lt;\u0026thinsp;2.2 \u0026times; 10\u003csup\u003e\u0026minus;\u0026thinsp;16\u003c/sup\u003e). Moreover, the correlation between ASE and dominant \u003cem\u003ecis\u003c/em\u003e-sQTL variants was higher than that between ASE and additive \u003cem\u003ecis\u003c/em\u003e-sQTL variants (Pearson r\u0026thinsp;=\u0026thinsp;0.7092; two-sided Student\u0026rsquo;s t-test: \u003cem\u003eP\u003c/em\u003e\u0026thinsp;\u0026lt;\u0026thinsp;2.2 \u0026times; 10\u003csup\u003e\u0026minus;\u0026thinsp;16\u003c/sup\u003e) (Fig.\u0026nbsp;\u003cspan refid=\"Fig3\" class=\"InternalRef\"\u003e3\u003c/span\u003eE, \u003cspan refid=\"Fig3\" class=\"InternalRef\"\u003e3\u003c/span\u003eF, and \u003cspan refid=\"Fig3\" class=\"InternalRef\"\u003e3\u003c/span\u003eG). Notably, in all 11 subsets with significant correlations of ASE effect sizes, including four subsets with a single breed, the \u003cspan class=\"InlineEquation\"\u003e\u003cspan class=\"mathinline\"\u003e\\(\\:{t}_{Dom}/{t}_{Add}\\)\u003c/span\u003e\u003c/span\u003e ratio was higher than the correlations with either additive or dominant cis-sQTL variants alone (Additional file 2: Fig. S11).\u003c/p\u003e\n\u003ch3\u003eFunctional annotation of ASE variants\u003c/h3\u003e\n\u003cp\u003eThe neutral theory of molecular evolution posits that most mutations at the molecular level are neutral, with minimal or no phenotypic impact. However, mutations occurring in conserved genomic regions can have significant phenotypic effects [\u003cspan citationid=\"CR23\" class=\"CitationRef\"\u003e23\u003c/span\u003e]. To explore the relationship between ASE variants and sequence conservation, we employed PhastCons and PhyloP scores from the UCSC Genome Browser [\u003cspan citationid=\"CR24\" class=\"CitationRef\"\u003e24\u003c/span\u003e]. Our analysis identified 16,746 ASE variants with phyloP scores greater than 2.0 and 14,307 ASE variants with phastCons scores above 0.8, representing approximately a twofold enrichment in conserved bases (Fig.\u0026nbsp;\u003cspan refid=\"Fig4\" class=\"InternalRef\"\u003e4\u003c/span\u003eA). Additionally, we observed that 93.57% of derived alleles among ASE variants exhibited negative effects (Additional file 2: Fig. S12), suggesting that most of these derived alleles are likely deleterious.\u003c/p\u003e \u003cp\u003e \u003c/p\u003e \u003cp\u003eASE analysis is a powerful approach for identifying \u003cem\u003ecis\u003c/em\u003e-regulatory elements involved in the regulation of diseases and other complex traits [\u003cspan citationid=\"CR13\" class=\"CitationRef\"\u003e13\u003c/span\u003e]. These regulatory elements are predominantly located within accessible chromatin regions, which can be identified using Assay for Transposase Accessible Chromatin Sequencing (ATAC-seq) [\u003cspan citationid=\"CR25\" class=\"CitationRef\"\u003e25\u003c/span\u003e]. To further investigate the relationship between ASE variants and chromatin accessibility, we analyzed ATAC-seq peak signals across 241 samples from 20 tissues obtained from the NCBI SRA database. The median sample size across these tissues was four, with a range from three samples in the hypothalamus, embryonic stem cells, spleen, and testes to 111 samples in the embryo (Additional file 1: Table S6). This analysis revealed that 50,645 ASE variants overlapped with ATAC-seq peaks, accounting for 31.44% of the total ASE variants. The proportion of ASE variants within ATAC-seq peaks ranged from 2.86\u0026ndash;13.39% across the 11 tissues with both RNA-seq and ATAC-seq data (Fig.\u0026nbsp;\u003cspan refid=\"Fig4\" class=\"InternalRef\"\u003e4\u003c/span\u003eB).\u003c/p\u003e \u003cp\u003eIn addition to ATAC-seq, Chromatin Immunoprecipitation Sequencing (ChIP-seq) is another method for globally identifying regulatory elements. ChIP-seq involves immunoprecipitation to selectively enrich DNA fragments associated with specific proteins, such as transcription factors or histone modifications. We examined the relationship between ASE variants and ChIP-seq peak signals across 158 samples from 20 tissues using five different antibodies. The median sample size across these tissues was three, with a range from two samples in the placenta and six other tissues to 103 samples in the mammary gland. The number of samples corresponding to the antibodies CTCF, H3K27ac, H3K27me3, H3K4me1, and H3K4me3 were 49, 90, 50, 153, and 157, respectively (Additional file 1: Table S7). Our analysis identified 90,556 ASE variants located within ChIP-seq peaks, accounting for 56.2% of all ASE variants. Specifically, the number of ASE variants within ChIP-seq peaks for the antibodies CTCF, H3K27ac, H3K27Me3, H3K4Me1, and H3K4Me3 were 45,828, 81,190, 57,010, 79,365, and 55,046, respectively. The proportion of ASE variants in ChIP-seq peaks varied from 0.7\u0026ndash;58.7% across the eight tissues with both RNA-seq and ChIP-seq data. Notably, among the five antibodies, the proportions of ASE variants located within peaks associated with H3K4me1 and H3K27ac, both markers of active enhancers [\u003cspan citationid=\"CR26\" class=\"CitationRef\"\u003e26\u003c/span\u003e], were substantially higher (Fig.\u0026nbsp;\u003cspan refid=\"Fig4\" class=\"InternalRef\"\u003e4\u003c/span\u003eC), indicating that ASE events are more frequently observed in enhancer regions.\u003c/p\u003e \u003cp\u003eFurthermore, we explored the relationship between detected ASE variants and QTL data from the cattle QTL Database [\u003cspan citationid=\"CR27\" class=\"CitationRef\"\u003e27\u003c/span\u003e] to investigate the potential impact of ASE variants on phenotypes. This analysis identified 7,196 ASE variants within QTL regions. When evaluating the enrichment of ASE variants relative to the QTL database, we found that traits with medium or low heritability, such as meat and carcass traits, as well as reproduction and health traits, exhibited higher levels of enrichment. In contrast, traits with medium or high heritability, including milk and exterior traits, showed lower levels of enrichment (Fig.\u0026nbsp;\u003cspan refid=\"Fig4\" class=\"InternalRef\"\u003e4\u003c/span\u003eD).\u003c/p\u003e\n\u003ch3\u003eRelationship between ASE variants and complex traits\u003c/h3\u003e\n\u003cp\u003eThe primary goal of the ASE atlas is to serve as a resource for elucidating the genetic mechanisms underlying complex traits. By focusing on shared variants between ASE and QTLs for milk traits (milk yield, milk fat percentage, milk fat yield, milk protein percentage, milk protein yield, and somatic cell score) in the cattle QTL Database [\u003cspan citationid=\"CR27\" class=\"CitationRef\"\u003e27\u003c/span\u003e], as well as SNPs associated with milk traits reported in a previous study [\u003cspan citationid=\"CR28\" class=\"CitationRef\"\u003e28\u003c/span\u003e], we designed multiplex PCR experiments to validate the relationship between milk traits and ASE variants in Holstein cattle (Fig.\u0026nbsp;\u003cspan refid=\"Fig5\" class=\"InternalRef\"\u003e5\u003c/span\u003eA and Materials and Methods). Of the 161 SNPs that were designed and successfully detected, 155 were retained after filtering for a minor allele frequency threshold of \u0026lt;\u0026thinsp;0.05 in a cohort of 1,052 cows. Using a genome-wide significant threshold (\u003cem\u003eP\u003c/em\u003e\u0026thinsp;\u0026lt;\u0026thinsp;5 \u0026times; 10\u003csup\u003e\u0026minus;\u0026thinsp;8\u003c/sup\u003e) [\u003cspan citationid=\"CR29\" class=\"CitationRef\"\u003e29\u003c/span\u003e] and a mixed linear model, we identified 28 SNPs with additive effects for six traits and 13 SNPs with dominant effects for five traits, including 12 SNPs with both additive and dominant effects (Fig.\u0026nbsp;\u003cspan refid=\"Fig5\" class=\"InternalRef\"\u003e5\u003c/span\u003eB and \u003cspan refid=\"Fig5\" class=\"InternalRef\"\u003e5\u003c/span\u003eC). Interestingly, the most significant locus was located on chromosome 14, a well-known locus responsible for milk traits, with the genes \u003cem\u003eCPSF1, GPAA1\u003c/em\u003e, \u003cem\u003ePLEC\u003c/em\u003e, \u003cem\u003eFAM83H\u003c/em\u003e, and \u003cem\u003eLYNX1\u003c/em\u003e mapped to this locus. Other candidate genes identified in our analysis include \u003cem\u003eCSN2\u003c/em\u003e, \u003cem\u003eGHR\u003c/em\u003e, \u003cem\u003eNUCB2\u003c/em\u003e, \u003cem\u003eELFN2\u003c/em\u003e, \u003cem\u003eTBC1D22A\u003c/em\u003e, \u003cem\u003eLGALS12\u003c/em\u003e, and \u003cem\u003eEP300\u003c/em\u003e (Additional file 2: Table S8).\u003c/p\u003e \u003cp\u003e \u003c/p\u003e \u003cp\u003eSimilar to eQTLs, we also investigated the relationship between the effect size of ASE variants and dominant QTLs, along with the degree of dominance. We found a strong, statistically significant correlation between the effect size of ASE variants and the \u003cspan class=\"InlineEquation\"\u003e\u003cspan class=\"mathinline\"\u003e\\(\\:{t}_{Dom}/{t}_{Add}\\)\u003c/span\u003e\u003c/span\u003e ratio (Pearson r\u0026thinsp;=\u0026thinsp;0.7459; two-sided Student\u0026rsquo;s t-test: \u003cem\u003eP\u003c/em\u003e\u0026thinsp;=\u0026thinsp;0.005) (Fig.\u0026nbsp;\u003cspan refid=\"Fig5\" class=\"InternalRef\"\u003e5\u003c/span\u003eD). Even with a suggestive threshold (\u003cem\u003eP\u003c/em\u003e\u0026thinsp;\u0026lt;\u0026thinsp;1 \u0026times; 10\u003csup\u003e\u0026minus;\u0026thinsp;6\u003c/sup\u003e) (Fig.\u0026nbsp;\u003cspan refid=\"Fig5\" class=\"InternalRef\"\u003e5\u003c/span\u003eE), the correlation remained significant (Pearson r\u0026thinsp;=\u0026thinsp;0.6843; two-sided Student\u0026rsquo;s t-test: \u003cem\u003eP\u003c/em\u003e\u0026thinsp;=\u0026thinsp;0.002). Unfortunately, owing to the limited number of QTLs identified for the same traits, the correlation between ASE effect sizes and dominant QTLs was not statistically significant.\u003c/p\u003e"},{"header":"Discussion","content":"\u003cp\u003eIn this study, we have constructed a comprehensive ASE Atlas, representing a valuable resource of regulatory variants across 35 distinct cattle tissues. Our approach involves an in silico protocol for ASE analysis, encompassing the generation, characterization, and functional annotation using epigenomic data, as well as the analysis of genetic effects (additive and dominant) based on publicly available datasets. While our methodology is both time-efficient and cost-effective, comparable to the CattleGTEx atlas protocol [\u003cspan citationid=\"CR4\" class=\"CitationRef\"\u003e4\u003c/span\u003e], there are potential improvements. For instance, collecting RNA-seq samples from multiple tissues within the same individuals could help mitigate spatiotemporal heterogeneity, thereby enhancing comparability. Furthermore, integrating whole genome sequencing data with RNA-seq would allow for the generation of personalized genomes, which could reduce alignment bias [\u003cspan citationid=\"CR16\" class=\"CitationRef\"\u003e16\u003c/span\u003e]. It would also enable the phasing of haplotypes, thereby increasing the power of ASE analysis [\u003cspan citationid=\"CR12\" class=\"CitationRef\"\u003e12\u003c/span\u003e], and facilitate the identification of intergenic or non-expressed variants, thereby providing a more precise understanding of the relationship between ASE and QTLs related to both physiological and molecular phenotypes. Additionally, matching epigenomic data with RNA-seq would offer deeper insights into the mechanisms driving ASE events.\u003c/p\u003e \u003cp\u003eOur findings highlight several distinctive features of ASE compared to eQTLs [\u003cspan additionalcitationids=\"CR4\" citationid=\"CR3\" class=\"CitationRef\"\u003e3\u003c/span\u003e\u0026ndash;\u003cspan citationid=\"CR5\" class=\"CitationRef\"\u003e5\u003c/span\u003e, \u003cspan citationid=\"CR21\" class=\"CitationRef\"\u003e21\u003c/span\u003e], particularly in terms of localization and tissue specificity. We observed a higher proportion of ASE variants in protein-coding genes and a higher degree of tissue specificity. This distinction likely arises because ASE variants are detected only in expressed sequences, whereas common variants in promotor regions are underrepresented. Specifically, only 1% of ASE variants were located in the promotor regions of protein-coding genes (within 2 kb upstream of the TSS), with 93% found downstream of the TSS. This contrasts with the ~\u0026thinsp;60% of eQTLs found upstream of the TSS [\u003cspan citationid=\"CR3\" class=\"CitationRef\"\u003e3\u003c/span\u003e]. Additionally, the reproducibility of ASE variants was lower than that of pig eQTLs [\u003cspan citationid=\"CR5\" class=\"CitationRef\"\u003e5\u003c/span\u003e] in independent datasets (49% vs. 77%) and across different tissues within the same BioProject (53% vs. 92%). These observations suggest that ASE may exhibit greater specificity than eQTLs.\u003c/p\u003e \u003cp\u003eSupporting this, we found that ASE variants are less frequent in promotor regions of protein-coding genes and more commonly located within ChIP-seq peaks marked by active enhancers (H3K4me1 and H3K27ac) rather than active promoters (H3K4me3) [\u003cspan citationid=\"CR26\" class=\"CitationRef\"\u003e26\u003c/span\u003e]. Moreover, ASE variants were less commonly found in ATAC-seq peaks, which are more typically associated with promoters [\u003cspan citationid=\"CR30\" class=\"CitationRef\"\u003e30\u003c/span\u003e, \u003cspan citationid=\"CR31\" class=\"CitationRef\"\u003e31\u003c/span\u003e], compared to ChIP-seq peaks marked by active enhancers. These results suggest that enhancer mutations may be a primary source of ASE events. The future availability of matched RNA-seq and H3K27ac histone ChIP-seq data will likely facilitate the identification of causative mutations in enhancer RNA [\u003cspan citationid=\"CR32\" class=\"CitationRef\"\u003e32\u003c/span\u003e], potentially driving ASE.\u003c/p\u003e \u003cp\u003eWhen examining the functional enrichment of ASE variants, we observed a significant presence of variants affecting splicing, consistent with our finding that ASE variants are more closely related to sQTLs than eQTLs. The highest enrichment was observed in variants leading to an in-frame stop codon. These findings strongly suggest that mutations altering mRNA splicing are key mechanisms underlying ASE [\u003cspan citationid=\"CR11\" class=\"CitationRef\"\u003e11\u003c/span\u003e]. Similar to eQTLs in humans [\u003cspan citationid=\"CR21\" class=\"CitationRef\"\u003e21\u003c/span\u003e] and pigs [\u003cspan citationid=\"CR5\" class=\"CitationRef\"\u003e5\u003c/span\u003e], ASE variants were also enriched in nonsynonymous, synonymous, stop-loss, and UTR variants, indicating the involvement of post-transcriptional and translational mechanisms, such as altered RNA-binding proteins, RNA stability, RNA editing, and the creation or disruption of upstream initiation codons [\u003cspan citationid=\"CR11\" class=\"CitationRef\"\u003e11\u003c/span\u003e, \u003cspan citationid=\"CR33\" class=\"CitationRef\"\u003e33\u003c/span\u003e].\u003c/p\u003e \u003cp\u003eOur study also showed that ASE variants are more enriched in dominant \u003cem\u003ecis\u003c/em\u003e-eQTLs and \u003cem\u003ecis\u003c/em\u003e-sQTLs compared to additive ones, with a higher correlation coefficient between ASE variant effect sizes and those of dominant \u003cem\u003ecis\u003c/em\u003e-eQTLs/\u003cem\u003ecis\u003c/em\u003e-sQTLs. Previous studies in maize reported that dominant expression patterns were more prevalent than additive patterns in ASE genes [\u003cspan citationid=\"CR19\" class=\"CitationRef\"\u003e19\u003c/span\u003e], and in poplar, 53% of ASE variants exhibited dominant effects on physiological and photosynthetic traits [\u003cspan citationid=\"CR34\" class=\"CitationRef\"\u003e34\u003c/span\u003e]. While research in rice suggests that biased expression of the favorable allele leads to dominant effects [\u003cspan citationid=\"CR18\" class=\"CitationRef\"\u003e18\u003c/span\u003e], our findings indicate a more complex relationship, where the degree of dominance in physiological and molecular phenotypes correlates with ASE effect sizes.\u003c/p\u003e \u003cp\u003eAmong the identified ASE variants, 14.3% exhibited opposite directional effects, particularly in immune-related tissues, genes, and GO terms, such as the MHC family. ASE variants with opposite directional effects in the MHC gene family have also been reported in human T cells [\u003cspan citationid=\"CR22\" class=\"CitationRef\"\u003e22\u003c/span\u003e]. Indeed, immune response eQTLs with opposite directional effects have also been observed in different contexts in humans [\u003cspan citationid=\"CR35\" class=\"CitationRef\"\u003e35\u003c/span\u003e]. Our study also found the number of ASE variants with opposite directional effects was closely related to dataset heterogeneity. Interestingly, ASE variants with opposite directional effects had significantly lower absolute effect sizes than those with consistent directional effects, suggesting that smaller effect sizes may facilitate directional shifts, thereby enhancing resilience in varying environments. This hypothesis is further supported by the relatively high variability in the number of ASE variants in immune-related tissues. Similar findings in maize suggest that ASE variants with opposite directional effects may contribute to adaptation in diverse environments [\u003cspan citationid=\"CR19\" class=\"CitationRef\"\u003e19\u003c/span\u003e], challenging the hypothesis that direction-shifting ASE causes overdominance in rice [\u003cspan citationid=\"CR18\" class=\"CitationRef\"\u003e18\u003c/span\u003e].\u003c/p\u003e \u003cp\u003eAlthough ASE has been widely used to explore the mechanism of heterosis [\u003cspan citationid=\"CR18\" class=\"CitationRef\"\u003e18\u003c/span\u003e, \u003cspan citationid=\"CR19\" class=\"CitationRef\"\u003e19\u003c/span\u003e], to our knowledge, this is the first study to report a linear correlation between ASE and the degree of dominance implicated in heterosis models [\u003cspan citationid=\"CR36\" class=\"CitationRef\"\u003e36\u003c/span\u003e, \u003cspan citationid=\"CR37\" class=\"CitationRef\"\u003e37\u003c/span\u003e]. We also found that tissues related to immune and reproductive functions, which are associated with high levels of heterosis [\u003cspan citationid=\"CR38\" class=\"CitationRef\"\u003e38\u003c/span\u003e], exhibited higher ASE effect sizes. Consistently, ASE variants were more enriched in functional traits with high levels of heterosis and low heritability [\u003cspan citationid=\"CR39\" class=\"CitationRef\"\u003e39\u003c/span\u003e] in the cattle QTL database. The predominance of derived alleles with negative effects supports the hypothesis that the dominance of advantageous ancestral alleles complements the deleterious effects of derived alleles [\u003cspan citationid=\"CR37\" class=\"CitationRef\"\u003e37\u003c/span\u003e]. Based on our discovery of a linear relationship between ASE and dominance, along with the hypothesis that overdominance results from allele interactions at a single locus [\u003cspan citationid=\"CR40\" class=\"CitationRef\"\u003e40\u003c/span\u003e], we propose that overdominance in physiological phenotypes stems from the interaction of two alleles at the molecular level, with certain loci manifesting overdominant effects.\u003c/p\u003e \u003cp\u003eTo validate our hypothesis inferred from the molecular phenotype of gene expression, we also observed the linear correlation between ASE effect size and the degree of dominance in QTLs affecting milk-related traits. Notably, the well-known QTL harboring \u003cem\u003eCPSF1\u003c/em\u003e on chromosome 14, affecting milk yield [\u003cspan citationid=\"CR41\" class=\"CitationRef\"\u003e41\u003c/span\u003e, \u003cspan citationid=\"CR42\" class=\"CitationRef\"\u003e42\u003c/span\u003e], exhibited both additive and dominant effects in our study. In fact, this dominant QTL was also observed in a previous study [\u003cspan citationid=\"CR43\" class=\"CitationRef\"\u003e43\u003c/span\u003e]. Other well-established candidate genes associated with milk traits, including \u003cem\u003eCSN2\u003c/em\u003e [\u003cspan citationid=\"CR41\" class=\"CitationRef\"\u003e41\u003c/span\u003e, \u003cspan citationid=\"CR44\" class=\"CitationRef\"\u003e44\u003c/span\u003e, \u003cspan citationid=\"CR45\" class=\"CitationRef\"\u003e45\u003c/span\u003e], \u003cem\u003eGHR\u003c/em\u003e [\u003cspan citationid=\"CR43\" class=\"CitationRef\"\u003e43\u003c/span\u003e, \u003cspan citationid=\"CR46\" class=\"CitationRef\"\u003e46\u003c/span\u003e], and \u003cem\u003eNUCB2\u003c/em\u003e [\u003cspan citationid=\"CR47\" class=\"CitationRef\"\u003e47\u003c/span\u003e], also exhibited dominant effects. Additionally, genes with limited or no previous reports in association with milk traits, such as \u003cem\u003eEP300, TBC1D22A\u003c/em\u003e, \u003cem\u003eELFN2\u003c/em\u003e, and \u003cem\u003eLGALS12\u003c/em\u003e, also displayed dominant effects. These results supported the hypothesis that dominance is pervasive in mammals [\u003cspan citationid=\"CR36\" class=\"CitationRef\"\u003e36\u003c/span\u003e]. Compared with previous studies, and thanks to the homogeneity of ASE data in RNA-seq, the variants we identified in association with milk traits are closer to causal.\u003c/p\u003e"},{"header":"Conclusions","content":"\u003cp\u003eOur comprehensive ASE atlas offers a valuable resource of regulatory variants for uncovering the genetic mechanisms underlying complex traits and enhancing economic traits through molecular selection. The characterization and functional annotation of ASE variants deepen our understanding of ASE formation and aid in identifying causal variants associated with complex traits. The observed correlation between ASE and dominance effects provides new insights into the genetic basis of heterosis, with promising applications in cattle production.\u003c/p\u003e"},{"header":"Materials and Methods","content":"\u003cdiv id=\"Sec10\" class=\"Section2\"\u003e \u003ch2\u003eData origin, alignment, and clustering\u003c/h2\u003e \u003cp\u003eA total of 33,999 RNA-seq entries were retrieved from the NCBI SRA database (4 July 2023) by searching for \u0026ldquo;cattle\u0026rdquo; and selecting \u0026ldquo;RNA\u0026rdquo; as the source. An in-house Perl script was used to refine this dataset by filtering for biological samples with the assay type of \u0026ldquo;RNA-seq\u0026rdquo; and the organisms of \u0026ldquo;\u003cem\u003eBos taurus\u003c/em\u003e, \u003cem\u003eBos indicus\u003c/em\u003e, or \u003cem\u003eBos taurus \u0026times; Bos indicus\u003c/em\u003e.\u0026rdquo; Further filtering was applied to include only samples sequenced on either the ILLUMINA or BGISEQ platforms. We also required a minimum of 20 samples per well-defined tissue type within each BioProject. This curation process yielded a final set of 8,518 transcriptomic datasets.\u003c/p\u003e \u003cp\u003eThe raw reads were processed using Trimmomatic v0.39 [\u003cspan citationid=\"CR48\" class=\"CitationRef\"\u003e48\u003c/span\u003e] with the parameters \u0026ldquo;LEADING:3 TRAILING:3 SLIDINGWINDOW:4:15 MINLEN:36\u0026rdquo; to ensure high-quality reads. The maximum read lengths were determined using seqkit v2.5.1 [\u003cspan citationid=\"CR49\" class=\"CitationRef\"\u003e49\u003c/span\u003e]. The filtered reads were aligned to the \u003cem\u003eBos taurus\u003c/em\u003e reference genome assembly (ARS-UCD1.2_Btau5.0.1Y) using STAR v2.7.10b [\u003cspan citationid=\"CR50\" class=\"CitationRef\"\u003e50\u003c/span\u003e] with the parameters \u0026ldquo;--outFilterMultimapNmax 1 --outFilterIntronMotifs RemoveNoncanonical Unannotated --outFilterMismatchNmax 10 --outSAMstrandField intronMotif --outSJfilterReads Unique.\u0026rdquo; Unaligned reads were further aligned using HISAT2 v2.2.1 [\u003cspan citationid=\"CR51\" class=\"CitationRef\"\u003e51\u003c/span\u003e]. The alignment files from STAR and HISAT2 were merged using Picard v3.0.0, and alignment statistics were computed using Samtools v1.17 (stat command) [\u003cspan citationid=\"CR52\" class=\"CitationRef\"\u003e52\u003c/span\u003e]. Samples were retained for downstream analysis if they met the criteria of a read mapping ratio\u0026thinsp;\u0026ge;\u0026thinsp;0.6 and a minimum of 50,000,000 mapped bases. We further filtered for sample size, requiring a minimum sample size of \u0026ge;\u0026thinsp;20 per clearly defined tissue within each BioProject, resulting in a final dataset of 7,532 samples.\u003c/p\u003e \u003cp\u003eGene expression levels were quantified as TPM using Stringtie v2.2.1 [\u003cspan citationid=\"CR53\" class=\"CitationRef\"\u003e53\u003c/span\u003e], with the GFF annotation file of the reference assembly. Hierarchical clustering of autosomal gene expression levels with an average TPM greater than 0.1 was performed using the R v4.2.3 package dendextend [\u003cspan citationid=\"CR54\" class=\"CitationRef\"\u003e54\u003c/span\u003e], employing Euclidean distance.\u003c/p\u003e \u003c/div\u003e \u003cdiv id=\"Sec11\" class=\"Section2\"\u003e \u003ch2\u003eSNP calling and imputation\u003c/h2\u003e \u003cp\u003eTo ensure high-quality variant calling, reads with a mapping quality score below 1 and unmapped reads were excluded using Samtools v1.17 (view command) [\u003cspan citationid=\"CR52\" class=\"CitationRef\"\u003e52\u003c/span\u003e]. The filtered alignment files were then processed with Picard tools, utilizing the ReorderSam, SortSam, and AddOrReplaceReadGroups commands. The SplitNCigarReads function in GATK v4.4.0.0 [\u003cspan citationid=\"CR55\" class=\"CitationRef\"\u003e55\u003c/span\u003e] was used to split reads at exon junctions and hard-clip overhanging sequences within intronic regions.\u003c/p\u003e \u003cp\u003eBase quality recalibration was performed using GATK\u0026rsquo;s BaseRecalibrator and ApplyBQSR modules. Variant calling was conducted using GATK\u0026rsquo;s HaplotypeCaller, followed by joint genotyping using the GenomicsDBImport and GenotypeGVCFs modules. SNPs were filtered using the VariantFiltration tool with the expression \u0026ldquo;QD\u0026thinsp;\u0026lt;\u0026thinsp;2.0 || FS\u0026thinsp;\u0026gt;\u0026thinsp;30.0.\u0026rdquo; The SelectVariants module of GATK was used to extract biallelic SNPs with a minor allele frequency of 0.01. Clean SNPs across all autosomes were merged using the MergeVcfs module of GATK and subsequently annotated using ANNOVAR v2016-02-01 [\u003cspan citationid=\"CR56\" class=\"CitationRef\"\u003e56\u003c/span\u003e].\u003c/p\u003e \u003cp\u003eTo further ensure data quality, SNPs were additionally filtered using VCFtools v0.1.17 [\u003cspan citationid=\"CR57\" class=\"CitationRef\"\u003e57\u003c/span\u003e] with the parameters \u0026ldquo;--maf 0.05 --max-missing 0.7 --minDP 10 --minGQ 20.\u0026rdquo; This filtering was applied both to the entire dataset of 7,532 samples and to 19 subsets, each containing more than 80 samples. Genotype imputation was then performed using Beagle v27Jan18.7e1 [\u003cspan citationid=\"CR58\" class=\"CitationRef\"\u003e58\u003c/span\u003e], leveraging a reference SNP panel [\u003cspan citationid=\"CR59\" class=\"CitationRef\"\u003e59\u003c/span\u003e]. After imputation, SNPs were filtered using BCFtools v1.17 [\u003cspan citationid=\"CR60\" class=\"CitationRef\"\u003e60\u003c/span\u003e] with thresholds of MAF\u0026thinsp;\u0026gt;\u0026thinsp;0.05 and DR2\u0026thinsp;\u0026ge;\u0026thinsp;0.8 to maintain high data quality. If imputed genotypes from the entire dataset of 7,532 samples were absent in the 19 subsets, they were incorporated into the respective subsets for further eQTL and related analyses.\u003c/p\u003e \u003cp\u003eHamming distances were calculated on the imputed genotypes of the entire dataset using PLINK v1.90b7.1 [\u003cspan citationid=\"CR61\" class=\"CitationRef\"\u003e61\u003c/span\u003e]. A neighbor-joining (NJ) tree was constructed based on these distances using MEGA11 [\u003cspan citationid=\"CR62\" class=\"CitationRef\"\u003e62\u003c/span\u003e] and then visualized using iTOL v6 [\u003cspan citationid=\"CR63\" class=\"CitationRef\"\u003e63\u003c/span\u003e]. Linkage disequilibrium (LD) analysis for each dataset was also conducted using PLINK v1.90b7.1.\u003c/p\u003e \u003c/div\u003e \u003cdiv id=\"Sec12\" class=\"Section2\"\u003e \u003ch2\u003eASE detection\u003c/h2\u003e \u003cp\u003eThe procedure for detecting ASE was adapted from a previous study [\u003cspan citationid=\"CR64\" class=\"CitationRef\"\u003e64\u003c/span\u003e] with slight modifications. Based on SNP data from all 7,532 samples, the ASEReadCounter module of GATK v4.4.0.0 was employed to count the reads corresponding to each allele at heterozygous sites for each individual. To be included in the analysis, both the reference and alternative allele counts were required to exceed three, and the minor allele read count ratio had to be greater than 0.01. A binomial test was then applied to these sites, followed by Benjamini-Hochberg multiple testing correction to control the FDR. Significant allelic imbalance at the individual level was determined by an FDR threshold of 0.05. At the population level, we focused on sites where at least six individuals were heterozygous. Sites showing significant allelic imbalance (FDR\u0026thinsp;\u0026lt;\u0026thinsp;0.05) in over 90% of these individuals were retained for further analysis. The effect size of allelic imbalance was quantified using the log allelic fold change (aFC), calculated according to a modified model from a previous study [\u003cspan citationid=\"CR65\" class=\"CitationRef\"\u003e65\u003c/span\u003e]. The calculations were performed as follows: (1) when the direction of allelic imbalance was consistent across individuals:\u003cdiv id=\"Equa\" class=\"Equation\"\u003e\u003cdiv format=\"TEX\" class=\"mathdisplay\" id=\"FileID_Equa\" name=\"EquationSource\"\u003e\n$$\\:{\\delta\\:}_{\\text{1,0}}=\\begin{array}{c}median\\\\\\:n=1...N\\end{array}\\frac{{c}_{1,n}}{{c}_{0,n}}$$\u003c/div\u003e\u003c/div\u003e\u003c/p\u003e \u003cp\u003eThe reported effect size:\u003cdiv id=\"Equb\" class=\"Equation\"\u003e\u003cdiv format=\"TEX\" class=\"mathdisplay\" id=\"FileID_Equb\" name=\"EquationSource\"\u003e\n$$\\:{s}_{\\text{1,0}}={{log}}_{2}{\\delta\\:}_{\\text{1,0}}$$\u003c/div\u003e\u003c/div\u003e\u003c/p\u003e \u003cp\u003e(2) when the directions of allelic imbalance varied among individuals:\u003cdiv id=\"Equc\" class=\"Equation\"\u003e\u003cdiv format=\"TEX\" class=\"mathdisplay\" id=\"FileID_Equc\" name=\"EquationSource\"\u003e\n$$\\:{\\delta\\:}_{\\text{1,0}}=\\left|{log}_{2}\\left(\\frac{{c}_{1,n}}{{c}_{0,n}}\\right)\\right|$$\u003c/div\u003e\u003c/div\u003e\u003c/p\u003e \u003cp\u003eThe reported effect size:\u003cdiv id=\"Equd\" class=\"Equation\"\u003e\u003cdiv format=\"TEX\" class=\"mathdisplay\" id=\"FileID_Equd\" name=\"EquationSource\"\u003e\n$$\\:{s}_{\\text{1,0}}=\\begin{array}{c}median\\\\\\:n=1...N\\end{array}{\\delta\\:}_{\\text{1,0}}$$\u003c/div\u003e\u003c/div\u003e\u003c/p\u003e \u003cp\u003eIn these equations, \u003cspan class=\"InlineEquation\"\u003e\u003cspan class=\"mathinline\"\u003e\\(\\:{c}_{1,n}\\)\u003c/span\u003e\u003c/span\u003e represents the read count for the alternative allele, \u003cspan class=\"InlineEquation\"\u003e\u003cspan class=\"mathinline\"\u003e\\(\\:{c}_{0,n}\\)\u003c/span\u003e\u003c/span\u003e represents the read count of the reference allele, \u003cspan class=\"InlineEquation\"\u003e\u003cspan class=\"mathinline\"\u003e\\(\\:n\\)\u003c/span\u003e\u003c/span\u003e is the index of the heterozygous site, and \u003cspan class=\"InlineEquation\"\u003e\u003cspan class=\"mathinline\"\u003e\\(\\:N\\)\u003c/span\u003e\u003c/span\u003e is the number of heterozygous sites showing significant allelic imbalance in the population. The distance from ASE variants to the TSS was defined based on annotations from ANNOVAR and the GFF annotation file of the reference assembly, measuring the distance from each ASE variant to the TSS of the nearest genes. GO enrichment analysis was conducted using the DAVID server [\u003cspan citationid=\"CR66\" class=\"CitationRef\"\u003e66\u003c/span\u003e].\u003c/p\u003e \u003c/div\u003e \u003cdiv id=\"Sec13\" class=\"Section2\"\u003e \u003ch2\u003eCovariate analysis for eQTL discovery\u003c/h2\u003e \u003cp\u003eFor the eQTL analysis, we employed a covariate analysis pipeline based on the methodology outlined in the CattleGTEx project [\u003cspan citationid=\"CR4\" class=\"CitationRef\"\u003e4\u003c/span\u003e]. First, gene expression levels for 25,365 expressed genes, each with a mean TPM value greater than 0.1, were normalized using the quantile-quantile normalization method [\u003cspan citationid=\"CR67\" class=\"CitationRef\"\u003e67\u003c/span\u003e]. To identify hidden factors contributing to transcriptome-wide variation in gene expression, Bayesian methods were employed to estimate latent covariates using PEER v1.0 [\u003cspan citationid=\"CR68\" class=\"CitationRef\"\u003e68\u003c/span\u003e]. The first ten PEER factors, where the cumulative posterior variances reached or nearly reached a plateau [\u003cspan citationid=\"CR4\" class=\"CitationRef\"\u003e4\u003c/span\u003e], were included as covariates. Additionally, to account for population structure in eQTL analysis, principal component (PC) analysis was performed based on the imputed genotypes using PLINK v1.90b7.1 [\u003cspan citationid=\"CR61\" class=\"CitationRef\"\u003e61\u003c/span\u003e]. Following the recommendations of CattleGTEx [\u003cspan citationid=\"CR4\" class=\"CitationRef\"\u003e4\u003c/span\u003e], the number of PCs included as covariates was determined by sample size: the first three PCs for sample sizes less than 150, the first five PCs for sample sizes between 150 and 249, and the first ten PCs for sample sizes of 250 or greater.\u003c/p\u003e \u003cp\u003e \u003cb\u003eCis\u003c/b\u003e \u003cb\u003e-eQTL mapping\u003c/b\u003e \u003c/p\u003e \u003cp\u003eFor \u003cem\u003ecis\u003c/em\u003e-eQTL mapping, normalized gene expression levels were used as the molecular phenotype, focusing on genes with a mean TPM value greater than 0.1. The analysis was conducted across 19 subsets, each containing more than 80 samples, using QTLtools v1.2 [\u003cspan citationid=\"CR69\" class=\"CitationRef\"\u003e69\u003c/span\u003e] with the parameters \u0026ldquo;--nominal 0.01.\u0026rdquo; To account for confounding factors, PEER factors and principal components were included as covariates in the cis-eQTL detection. A significance threshold of \u003cem\u003eP\u003c/em\u003e\u0026thinsp;\u0026lt;\u0026thinsp;1 \u0026times; 10\u003csup\u003e\u0026minus;\u0026thinsp;6\u003c/sup\u003e was applied for SNP associations. When evaluating additive effects, genotypes were coded as follows: homozygous reference (0), heterozygous (1), and homozygous alternative (2). For dominant effects, the genotypes were coded as homozygous reference (0), heterozygous (1), and homozygous alternative (0). The effect size of \u003cem\u003ecis\u003c/em\u003e-eQTL was described using log allelic fold change (aFC), calculated according to a modified linear regression model [\u003cspan citationid=\"CR65\" class=\"CitationRef\"\u003e65\u003c/span\u003e]:\u003cdiv id=\"Eque\" class=\"Equation\"\u003e\u003cdiv format=\"TEX\" class=\"mathdisplay\" id=\"FileID_Eque\" name=\"EquationSource\"\u003e\n$$\\:y={\\beta\\:}_{0}+{\\beta\\:}_{1}{X}_{1}+{\\beta\\:}_{2}{X}_{2}+\\epsilon\\:$$\u003c/div\u003e\u003c/div\u003e\u003c/p\u003e \u003cp\u003eWhere \u003cspan class=\"InlineEquation\"\u003e\u003cspan class=\"mathinline\"\u003e\\(\\:y\\)\u003c/span\u003e\u003c/span\u003e denotes the normalized molecular phenotype, \u003cspan class=\"InlineEquation\"\u003e\u003cspan class=\"mathinline\"\u003e\\(\\:{\\beta\\:}_{0}\\)\u003c/span\u003e\u003c/span\u003e represents the intercept, \u003cspan class=\"InlineEquation\"\u003e\u003cspan class=\"mathinline\"\u003e\\(\\:{\\beta\\:}_{1}\\)\u003c/span\u003e\u003c/span\u003e captures the effect of SNP markers, \u003cspan class=\"InlineEquation\"\u003e\u003cspan class=\"mathinline\"\u003e\\(\\:{X}_{1}\\)\u003c/span\u003e\u003c/span\u003e represents the maker genotypes, \u003cspan class=\"InlineEquation\"\u003e\u003cspan class=\"mathinline\"\u003e\\(\\:{\\beta\\:}_{2}\\)\u003c/span\u003e\u003c/span\u003e accounts for the effect of covariates, \u003cspan class=\"InlineEquation\"\u003e\u003cspan class=\"mathinline\"\u003e\\(\\:{X}_{2}\\)\u003c/span\u003e\u003c/span\u003e represents the covariates, \u003cspan class=\"InlineEquation\"\u003e\u003cspan class=\"mathinline\"\u003e\\(\\:\\epsilon\\:\\)\u003c/span\u003e\u003c/span\u003e represents the random residuals. The effect size (\u003cspan class=\"InlineEquation\"\u003e\u003cspan class=\"mathinline\"\u003e\\(\\:{s}_{\\text{1,0}}\\)\u003c/span\u003e\u003c/span\u003e) was reported as:\u003cdiv id=\"Equf\" class=\"Equation\"\u003e\u003cdiv format=\"TEX\" class=\"mathdisplay\" id=\"FileID_Equf\" name=\"EquationSource\"\u003e\n$$\\:{\\delta\\:}_{\\text{1,0}}=\\frac{\\:2{\\beta\\:}_{1}}{{\\beta\\:}_{0}}+1$$\u003c/div\u003e\u003c/div\u003e\u003cdiv id=\"Equg\" class=\"Equation\"\u003e\u003cdiv format=\"TEX\" class=\"mathdisplay\" id=\"FileID_Equg\" name=\"EquationSource\"\u003e\n$$\\:{s}_{\\text{1,0}}={{log}}_{2}{\\delta\\:}_{\\text{1,0}}$$\u003c/div\u003e\u003c/div\u003e\u003c/p\u003e \u003cp\u003eThe degree of dominance was assessed by calculating the ratios \u003cspan class=\"InlineEquation\"\u003e\u003cspan class=\"mathinline\"\u003e\\(\\:{t}_{Dom}/{t}_{Add}\\)\u003c/span\u003e\u003c/span\u003e according to established methods [\u003cspan citationid=\"CR20\" class=\"CitationRef\"\u003e20\u003c/span\u003e]. The T-statistics for additive and dominant effects, \u003cspan class=\"InlineEquation\"\u003e\u003cspan class=\"mathinline\"\u003e\\(\\:{t}_{Add}\\)\u003c/span\u003e\u003c/span\u003e and \u003cspan class=\"InlineEquation\"\u003e\u003cspan class=\"mathinline\"\u003e\\(\\:{t}_{Dom}\\)\u003c/span\u003e\u003c/span\u003e, were computed as follows:\u003cdiv id=\"Equh\" class=\"Equation\"\u003e\u003cdiv format=\"TEX\" class=\"mathdisplay\" id=\"FileID_Equh\" name=\"EquationSource\"\u003e\n$$\\:{t}_{Add}=\\frac{{\\beta\\:}_{Add}}{se\\left({\\beta\\:}_{Add}\\right)}\\:\\:\\:\\:\\:\\:\\:\\:\\:\\:\\:{t}_{Dom}=\\frac{{\\beta\\:}_{Dom}}{se\\left({\\beta\\:}_{Dom}\\right)}\\:$$\u003c/div\u003e\u003c/div\u003e\u003c/p\u003e \u003cp\u003eHere, \u003cspan class=\"InlineEquation\"\u003e\u003cspan class=\"mathinline\"\u003e\\(\\:{\\beta\\:}_{Add}\\)\u003c/span\u003e\u003c/span\u003e and \u003cspan class=\"InlineEquation\"\u003e\u003cspan class=\"mathinline\"\u003e\\(\\:{\\beta\\:}_{Dom}\\)\u003c/span\u003e\u003c/span\u003e represent the additive and dominant effects, respectively, while \u003cspan class=\"InlineEquation\"\u003e\u003cspan class=\"mathinline\"\u003e\\(\\:se\\left({\\beta\\:}_{Add}\\right)\\)\u003c/span\u003e\u003c/span\u003e and \u003cspan class=\"InlineEquation\"\u003e\u003cspan class=\"mathinline\"\u003e\\(\\:se\\left({\\beta\\:}_{Dom}\\right)\\)\u003c/span\u003e\u003c/span\u003e denote their corresponding standard errors. If the effect of the minor allele is negative, the degree of dominance is defined as \u003cspan class=\"InlineEquation\"\u003e\u003cspan class=\"mathinline\"\u003e\\(\\:-{t}_{Dom}/{t}_{Add}\\)\u003c/span\u003e\u003c/span\u003e; otherwise, it is defined as \u003cspan class=\"InlineEquation\"\u003e\u003cspan class=\"mathinline\"\u003e\\(\\:{t}_{Dom}/{t}_{Add}\\)\u003c/span\u003e\u003c/span\u003e.\u003c/p\u003e \u003cp\u003e \u003cb\u003eCis\u003c/b\u003e \u003cb\u003e-sQTL mapping\u003c/b\u003e \u003c/p\u003e \u003cp\u003eFor \u003cem\u003ecis\u003c/em\u003e-sQTL mapping, splicing junctions were identified from alignment files using the Leafcutter v0.2.7 tool [\u003cspan citationid=\"CR70\" class=\"CitationRef\"\u003e70\u003c/span\u003e]. Initially, the bam2junc.sh script was used to convert each individual\u0026rsquo;s BAM file into a junction file. The identified junction files were then clustered across the population using the leafcutter_cluster.py script, with a minimum of 50 reads per junction and a maximum intron length of 500,000 base pairs. Phenotype tables required for downstream splicing analysis were prepared using the prepare_phenotype_table.py script from the Leafcutter tool and subsequently normalized using quantile-quantile normalization. The covariates and analytical methods for \u003cem\u003ecis\u003c/em\u003e-sQTL mapping were identical to those employed for \u003cem\u003ecis\u003c/em\u003e-eQTL mapping.\u003c/p\u003e \u003c/div\u003e \u003cdiv id=\"Sec14\" class=\"Section2\"\u003e \u003ch2\u003eOverlap analysis of ASE variants with conservation sites\u003c/h2\u003e \u003cp\u003eBased on the multiple alignment format (maf) file from UCSC [\u003cspan citationid=\"CR24\" class=\"CitationRef\"\u003e24\u003c/span\u003e] (\u003cspan class=\"ExternalRef\"\u003e\u003cspan class=\"RefSource\"\u003ehttps://hgdownload.soe.ucsc.edu/goldenPath/hg38/multiz470way/maf/\u003c/span\u003e\u003cspan address=\"https://hgdownload.soe.ucsc.edu/goldenPath/hg38/multiz470way/maf/\" targettype=\"URL\" class=\"RefTarget\"\u003e\u003c/span\u003e\u003c/span\u003e), we established the correspondence between human autosomal physical positions and the cattle ARS-UCD1.2 assembly using an in-house Perl script. Conservation scores for humans were retrieved from UCSC (\u003cspan class=\"ExternalRef\"\u003e\u003cspan class=\"RefSource\"\u003ehttps://hgdownload.soe.ucsc.edu/goldenPath/hg38/phyloP470way/hg38.470way.phyloP/\u003c/span\u003e\u003cspan address=\"https://hgdownload.soe.ucsc.edu/goldenPath/hg38/phyloP470way/hg38.470way.phyloP/\" targettype=\"URL\" class=\"RefTarget\"\u003e\u003c/span\u003e\u003c/span\u003e and \u003cspan class=\"ExternalRef\"\u003e\u003cspan class=\"RefSource\"\u003ehttps://hgdownload.soe.ucsc.edu/goldenPath/hg38/phastCons470way/hg38.470way.phastCons/\u003c/span\u003e\u003cspan address=\"https://hgdownload.soe.ucsc.edu/goldenPath/hg38/phastCons470way/hg38.470way.phastCons/\" targettype=\"URL\" class=\"RefTarget\"\u003e\u003c/span\u003e\u003c/span\u003e) were subsequently translated to corresponding cattle positions using the established mappings. We then extracted the relevant conservation scores for ASE variants in cattle through another custom Perl script.\u003c/p\u003e \u003cp\u003eTo infer ancestral alleles, we obtained the whole-genome sequences of five outgroup species (American Bison, Banteng, Yak, European Bison, and Gaur) [\u003cspan citationid=\"CR71\" class=\"CitationRef\"\u003e71\u003c/span\u003e] (Additional file 1: Table S9). After performing quality control similar to the RNA-seq data processing, we aligned the trimmed reads to the Bos taurus reference assembly using bwa-mem v0.7.17-r1188 [\u003cspan citationid=\"CR46\" class=\"CitationRef\"\u003e46\u003c/span\u003e] with default parameters. The alignment files were sorted by coordinate using Picard's SortSam tool, and duplicate reads were marked and removed using Picard's MarkDuplicates tool. Variant calling was conducted using BCFtools v1.17 mpileup with parameters \u0026ldquo;-q 30 -C 50 -Q 20 -B.\u0026rdquo; The ancestral allele was inferred using a consensus sequence from at least two of the outgroup species through an in-house Perl script.\u003c/p\u003e \u003c/div\u003e \u003cdiv id=\"Sec15\" class=\"Section2\"\u003e \u003ch2\u003eATAC-seq analysis\u003c/h2\u003e \u003cp\u003eWe retrieved 528 entries of run information from the NCBI SRA database (14 December 2023) by searching for \"cattle ATAC-seq,\u0026rdquo; specifying \"illumina\" as the platform, \u0026ldquo;DNA\u0026rdquo; as the source, \u0026ldquo;EpiGenomics\u0026rdquo; as the strategy, and \u0026ldquo;\u003cem\u003eBos taurus\u003c/em\u003e\u0026rdquo; as the organism. An in-house Perl script was used to further refine the dataset, applying criteria including a minimum sample size of three and clearly defined tissue information within each BioProject, resulting in 241 ATAC-seq data. The raw reads were quality-trimmed, and adapters were removed using Trim Galore v0.6.10 [\u003cspan citationid=\"CR72\" class=\"CitationRef\"\u003e72\u003c/span\u003e] with the parameters \u0026ldquo;--q 25 --phred33 --stringency 3 --length 35 -e 0.1.\u0026rdquo; The trimmed reads were aligned to the cattle reference genome (ARS-UCD1.2) using Bowtie2 v2.5.3 [\u003cspan citationid=\"CR73\" class=\"CitationRef\"\u003e73\u003c/span\u003e] with the parameters \u0026ldquo;--very-sensitive -X 2000.\u0026rdquo; Reads mapped to the sex chromosomes (X and Y) were filtered out. Duplicate reads were marked and removed using Sambamba v1.0.1 [\u003cspan citationid=\"CR74\" class=\"CitationRef\"\u003e74\u003c/span\u003e], followed by sorting the BAM files by coordinate and read name using Samtools v1.17 [\u003cspan citationid=\"CR52\" class=\"CitationRef\"\u003e52\u003c/span\u003e]. Peaks representing open chromatin regions were called using MACS3 v3.0.0 [\u003cspan citationid=\"CR75\" class=\"CitationRef\"\u003e75\u003c/span\u003e] with the parameters \u0026ldquo;--shift 75 --extsize 150 --nomodel --call-summits --nolambda --keep-dup all -q 0.01\u0026rdquo;. We then used the intersect command of bedtools v2.26.0 [\u003cspan citationid=\"CR76\" class=\"CitationRef\"\u003e76\u003c/span\u003e] to obtain overlapping peak regions from the same tissue type, which were subsequently utilized to identify ASE variants located within ATAC-seq peaks.\u003c/p\u003e \u003c/div\u003e \u003cdiv id=\"Sec16\" class=\"Section2\"\u003e \u003ch2\u003eChIP-seq analysis\u003c/h2\u003e \u003cp\u003eA total of 5,251 entries of run information were retrieved from the NCBI SRA database (6 May 2024) by searching for \"cattle ChIP-seq,\u0026rdquo; specifying \"illumina\" as the platform, \u0026ldquo;DNA\u0026rdquo; as the source, \u0026ldquo;EpiGenomics\u0026rdquo; as the strategy, and \u0026ldquo;\u003cem\u003eBos taurus\u003c/em\u003e\u0026rdquo; as the organism. Using an in-house Perl script, we selected samples from four BioProjects (PRJEB41939 [\u003cspan citationid=\"CR77\" class=\"CitationRef\"\u003e77\u003c/span\u003e], PRJEB52456 [\u003cspan citationid=\"CR78\" class=\"CitationRef\"\u003e78\u003c/span\u003e], PRJEB53044 [\u003cspan citationid=\"CR79\" class=\"CitationRef\"\u003e79\u003c/span\u003e] and PRJEB6906 [\u003cspan citationid=\"CR80\" class=\"CitationRef\"\u003e80\u003c/span\u003e]) that had at least two samples per tissue per antibody, resulting in 499 ChIP-seq datasets. Both input and antibody-treated samples were initially processed using Trim Galore v0.6.10 with the parameters \u0026ldquo;--q 25 --phred33 --length 25 -e 0.1 --stringency 4 --paired\u0026rdquo; to remove low-quality bases and adapters. The trimmed reads were aligned to the reference genome (ARS-UCD1.2) using Bowtie2 v2.5.3, and the resulting SAM files were converted to BAM format using Samtools v1.17. Duplicate reads were marked and removed using Sambamba v1.0.1. The deduplicated BAM files were then sorted using Samtools v1.17. To identify significant DNA-protein interaction regions, peaks were called using MACS3 v3.0.0 [\u003cspan citationid=\"CR75\" class=\"CitationRef\"\u003e75\u003c/span\u003e] with a stringent cutoff (-q 0.01) to ensure high-confidence detection. Antibody-treated samples were compared against input controls to filter out background signals and enhance the specificity of the detected peaks. As with ATAC-seq, we used the intersect command of bedtools v2.26.0 [\u003cspan citationid=\"CR76\" class=\"CitationRef\"\u003e76\u003c/span\u003e] to find overlapping peak regions from the same tissue and antibody type and then identified ASE variants located within these ChIP-seq peaks.\u003c/p\u003e \u003c/div\u003e \u003cdiv id=\"Sec17\" class=\"Section2\"\u003e \u003ch2\u003eValidation of shared variants between ASE and QTL for milk traits\u003c/h2\u003e \u003cp\u003eAmong the 8,733 SNPs associated with milk traits reported by a previous study [\u003cspan citationid=\"CR28\" class=\"CitationRef\"\u003e28\u003c/span\u003e], 156 overlapped with ASE variants. Additionally, of the 15,958 SNPs linked to milk traits in the cattle QTL database, 104 overlapped with ASE variants. After merging these two sets, we retained 248 SNPs. Filtering out those with a minor allele frequency\u0026thinsp;\u0026lt;\u0026thinsp;0.05 in Holstein cattle based on a reference panel from a previous study [\u003cspan citationid=\"CR59\" class=\"CitationRef\"\u003e59\u003c/span\u003e] reduced the number to 208 SNPs. Further pruning using PLINK with the parameter \u0026ldquo;--indep-pairwise 5 1 0.8\u0026rdquo; in the reference panel resulted in a final set of 171 SNPs. An additional 10 SNPs were discarded due to being flanked by repeat sequences as determined by Adsen Biotechnology Co., Ltd. (Urumchi, China) during genotyping.\u003c/p\u003e \u003cp\u003eA total of 1,052 Chinese Holstein cows from Xinjiang were used in this study. Nine milk-related traits were considered: milk yield (kg), milk fat percentage (%), milk fat yield (kg), milk protein percentage (%), milk protein yield (kg), somatic cell count (10\u003csup\u003e4\u003c/sup\u003e/ml), milk lactose percentage (%), total solids percentage (%), and urea nitrogen concentration (mg/dl). The total number of test-day records was 8,839. Descriptive statistics of the test-day traits are provided in Additional file 2: Table S10.\u003c/p\u003e \u003cp\u003eWhole blood was collected from the jugular vein of each individual using a blood collection needle by an experienced veterinarian and then placed in a blood collection tube containing EDTA. The genomic DNA was extracted using the standard phenol-chloroform procedure, and then its quantity and quality were assessed using agarose gel electrophoresis and a Qubit fluorometer (Invitrogen, Carlsbad, USA), respectively. The qualified DNA was transported to Adsen Biotechnology Co., Ltd. for multiplex PCR experiments, following sequencing and variant detection (Additional file 2: Supplementary Materials and Methods).\u003c/p\u003e \u003cp\u003eWe used the following mixed linear model, implemented with the lmerTest package in R [\u003cspan citationid=\"CR6\" class=\"CitationRef\"\u003e6\u003c/span\u003e].\u003cdiv id=\"Equi\" class=\"Equation\"\u003e\u003cdiv format=\"TEX\" class=\"mathdisplay\" id=\"FileID_Equi\" name=\"EquationSource\"\u003e\n$$\\:y={\\beta\\:}_{0}+{\\beta\\:}_{1}{X}_{1}+{\\beta\\:}_{2}{X}_{2}+{\\beta\\:}_{3}{X}_{3}+\\epsilon\\:$$\u003c/div\u003e\u003c/div\u003e\u003c/p\u003e \u003cp\u003eWhere \u003cspan class=\"InlineEquation\"\u003e\u003cspan class=\"mathinline\"\u003e\\(\\:y\\)\u003c/span\u003e\u003c/span\u003e denotes the phenotype, \u003cspan class=\"InlineEquation\"\u003e\u003cspan class=\"mathinline\"\u003e\\(\\:{\\beta\\:}_{0}\\)\u003c/span\u003e\u003c/span\u003e denotes the intercept of the linear regression, \u003cspan class=\"InlineEquation\"\u003e\u003cspan class=\"mathinline\"\u003e\\(\\:{\\beta\\:}_{1}\\)\u003c/span\u003e\u003c/span\u003e represents the effects of SNP markers, and \u003cspan class=\"InlineEquation\"\u003e\u003cspan class=\"mathinline\"\u003e\\(\\:{X}_{1}\\)\u003c/span\u003e\u003c/span\u003e denotes the maker genotypes. The term \u003cspan class=\"InlineEquation\"\u003e\u003cspan class=\"mathinline\"\u003e\\(\\:{\\beta\\:}_{2}\\)\u003c/span\u003e\u003c/span\u003e represents the fixed effects of herd, parity, and milk month, while \u003cspan class=\"InlineEquation\"\u003e\u003cspan class=\"mathinline\"\u003e\\(\\:{X}_{2}\\)\u003c/span\u003e\u003c/span\u003e denotes these factors. The term \u003cspan class=\"InlineEquation\"\u003e\u003cspan class=\"mathinline\"\u003e\\(\\:{\\beta\\:}_{3}\\)\u003c/span\u003e\u003c/span\u003e represents the random effect of test day, with \u003cspan class=\"InlineEquation\"\u003e\u003cspan class=\"mathinline\"\u003e\\(\\:{X}_{3}\\)\u003c/span\u003e\u003c/span\u003e denoting the test day, and \u003cspan class=\"InlineEquation\"\u003e\u003cspan class=\"mathinline\"\u003e\\(\\:\\epsilon\\:\\)\u003c/span\u003e\u003c/span\u003e represents random residuals. Descriptive statistics for the fixed and random effects are provided in Additional file 2: Table S11. Similar to \u003cem\u003ecis\u003c/em\u003e-eQTL analysis, for examining additive effects, genotypes were coded as follows: homozygous reference (0), heterozygous (1), and homozygous alternative (2). For dominant effects, genotypes were coded as homozygous reference (0), heterozygous (1), and homozygous alternative (0). The ratios of \u003cspan class=\"InlineEquation\"\u003e\u003cspan class=\"mathinline\"\u003e\\(\\:{t}_{Dom}/{t}_{Add}\\)\u003c/span\u003e\u003c/span\u003e, used to measure the degree of dominance, were calculated as described previously [\u003cspan citationid=\"CR20\" class=\"CitationRef\"\u003e20\u003c/span\u003e] and were consistent with the \u003cem\u003ecis\u003c/em\u003e-eQTL analysis.\u003c/p\u003e \u003c/div\u003e"},{"header":"Abbreviations","content":"\u003cdiv class=\"DefinitionList\"\u003e \u003cdiv class=\"DefinitionListEntry\"\u003e \u003cdiv class=\"Term\"\u003eASE\u003c/div\u003e \u003cdiv class=\"Description\"\u003e \u003cp\u003eAllelic-specific expression\u003c/p\u003e \u003c/div\u003e \u003c/div\u003e \u003cdiv class=\"DefinitionListEntry\"\u003e \u003cdiv class=\"Term\"\u003eQTL\u003c/div\u003e \u003cdiv class=\"Description\"\u003e \u003cp\u003equantitative trait loci\u003c/p\u003e \u003c/div\u003e \u003c/div\u003e \u003cdiv class=\"DefinitionListEntry\"\u003e \u003cdiv class=\"Term\"\u003eeQTL\u003c/div\u003e \u003cdiv class=\"Description\"\u003e \u003cp\u003eexpression QTL\u003c/p\u003e \u003c/div\u003e \u003c/div\u003e \u003cdiv class=\"DefinitionListEntry\"\u003e \u003cdiv class=\"Term\"\u003esQTL\u003c/div\u003e \u003cdiv class=\"Description\"\u003e \u003cp\u003esplicing QTL\u003c/p\u003e \u003c/div\u003e \u003c/div\u003e \u003cdiv class=\"DefinitionListEntry\"\u003e \u003cdiv class=\"Term\"\u003eGWAS\u003c/div\u003e \u003cdiv class=\"Description\"\u003e \u003cp\u003egenome-wide association study\u003c/p\u003e \u003c/div\u003e \u003c/div\u003e \u003cdiv class=\"DefinitionListEntry\"\u003e \u003cdiv class=\"Term\"\u003eGTEx\u003c/div\u003e \u003cdiv class=\"Description\"\u003e \u003cp\u003eGenotype-Tissue Expression\u003c/p\u003e \u003c/div\u003e \u003c/div\u003e \u003cdiv class=\"DefinitionListEntry\"\u003e \u003cdiv class=\"Term\"\u003eENCODE\u003c/div\u003e \u003cdiv class=\"Description\"\u003e \u003cp\u003eEncyclopedia of DNA Elements\u003c/p\u003e \u003c/div\u003e \u003c/div\u003e \u003cdiv class=\"DefinitionListEntry\"\u003e \u003cdiv class=\"Term\"\u003eFAANG\u003c/div\u003e \u003cdiv class=\"Description\"\u003e \u003cp\u003eFunctional Annotation of Animal Genomes\u003c/p\u003e \u003c/div\u003e \u003c/div\u003e \u003cdiv class=\"DefinitionListEntry\"\u003e \u003cdiv class=\"Term\"\u003eTPM\u003c/div\u003e \u003cdiv class=\"Description\"\u003e \u003cp\u003etranscripts per million\u003c/p\u003e \u003c/div\u003e \u003c/div\u003e \u003cdiv class=\"DefinitionListEntry\"\u003e \u003cdiv class=\"Term\"\u003eCattleGTEx\u003c/div\u003e \u003cdiv class=\"Description\"\u003e \u003cp\u003eCattle Genotype-Tissue Expression\u003c/p\u003e \u003c/div\u003e \u003c/div\u003e \u003cdiv class=\"DefinitionListEntry\"\u003e \u003cdiv class=\"Term\"\u003eUTR\u003c/div\u003e \u003cdiv class=\"Description\"\u003e \u003cp\u003euntranslated region\u003c/p\u003e \u003c/div\u003e \u003c/div\u003e \u003cdiv class=\"DefinitionListEntry\"\u003e \u003cdiv class=\"Term\"\u003eTSS\u003c/div\u003e \u003cdiv class=\"Description\"\u003e \u003cp\u003etranscriptional start site\u003c/p\u003e \u003c/div\u003e \u003c/div\u003e \u003cdiv class=\"DefinitionListEntry\"\u003e \u003cdiv class=\"Term\"\u003eATAC-seq\u003c/div\u003e \u003cdiv class=\"Description\"\u003e \u003cp\u003eAssay for Transposase Accessible Chromatin sequencing\u003c/p\u003e \u003c/div\u003e \u003c/div\u003e \u003cdiv class=\"DefinitionListEntry\"\u003e \u003cdiv class=\"Term\"\u003eChIP-seq\u003c/div\u003e \u003cdiv class=\"Description\"\u003e \u003cp\u003eChromatin Immunoprecipitation Sequencing\u003c/p\u003e \u003c/div\u003e \u003c/div\u003e \u003cdiv class=\"DefinitionListEntry\"\u003e \u003cdiv class=\"Term\"\u003eLD\u003c/div\u003e \u003cdiv class=\"Description\"\u003e \u003cp\u003eLinkage disequilibrium\u003c/p\u003e \u003c/div\u003e \u003c/div\u003e \u003cdiv class=\"DefinitionListEntry\"\u003e \u003cdiv class=\"Term\"\u003eMAF\u003c/div\u003e \u003cdiv class=\"Description\"\u003e \u003cp\u003eminor allele frequency\u003c/p\u003e \u003c/div\u003e \u003c/div\u003e \u003cdiv class=\"DefinitionListEntry\"\u003e \u003cdiv class=\"Term\"\u003eNJ\u003c/div\u003e \u003cdiv class=\"Description\"\u003e \u003cp\u003eneighbor-joining\u003c/p\u003e \u003c/div\u003e \u003c/div\u003e \u003cdiv class=\"DefinitionListEntry\"\u003e \u003cdiv class=\"Term\"\u003eFDR\u003c/div\u003e \u003cdiv class=\"Description\"\u003e \u003cp\u003efalse discovery rate\u003c/p\u003e \u003c/div\u003e \u003c/div\u003e \u003cdiv class=\"DefinitionListEntry\"\u003e \u003cdiv class=\"Term\"\u003eaFC\u003c/div\u003e \u003cdiv class=\"Description\"\u003e \u003cp\u003eallelic fold change\u003c/p\u003e \u003c/div\u003e \u003c/div\u003e \u003cdiv class=\"DefinitionListEntry\"\u003e \u003cdiv class=\"Term\"\u003eGO\u003c/div\u003e \u003cdiv class=\"Description\"\u003e \u003cp\u003eGene Ontology\u003c/p\u003e \u003c/div\u003e \u003c/div\u003e \u003cdiv class=\"DefinitionListEntry\"\u003e \u003cdiv class=\"Term\"\u003ePC\u003c/div\u003e \u003cdiv class=\"Description\"\u003e \u003cp\u003eprincipal component.\u003c/p\u003e \u003c/div\u003e \u003c/div\u003e \u003c/div\u003e"},{"header":"Declarations","content":"\u003cp\u003e\u003cstrong\u003eAcknowledgements\u003c/strong\u003e\u003c/p\u003e\n\u003cp\u003eWe would like to express our sincere gratitude to Professor Tom Druet from University of Li\u0026egrave;ge, Professor Yu Wang from Northwest A\u0026amp;F University and Associate Professor Han Xu from Anhui Agricultural University for their valuable suggestions and insightful contributions to this study.\u003c/p\u003e\n\u003cp\u003e\u003cstrong\u003eAuthors\u0026rsquo; contributions\u003c/strong\u003e\u003c/p\u003e\n\u003cp\u003eQ.C. conceptualized and designed the study. Q.C. and J.L. conducted the allele-specific expression analysis. X.L. carried out the ATAC-seq and ChIP-seq analyses. Q.C. and L.X. performed mixed linear model analysis on milk related traits. L.L. prepared the schematic diagram of cattle tissues. L.X. coordinated with the Holstein farm for phenotype data collection and blood sampling. X.H. secured funding and supervised the entire study. Q.C. wrote the manuscript.\u0026nbsp;\u003c/p\u003e\n\u003cp\u003e\u003cstrong\u003eFunding\u003c/strong\u003e\u003c/p\u003e\n\u003cp\u003eThis study was financially supported by the National Key R\u0026amp;D Program of China (2021YFD1200903).\u003c/p\u003e\n\u003cp\u003e\u003cstrong\u003eData Availability\u003c/strong\u003e\u003c/p\u003e\n\u003cp\u003eAll raw data analyzed in this study are publicly available for download without restrictions from the NCBI SRA database (https://www.ncbi.nlm.nih.gov/sra/). Details of RNA-seq, ATAC-seq, ChIP-seq and whole genome sequence can be found in Additional file 1: Table S1, S6, S7 and S9, respectively. All the computational scripts and codes for RNA-seq, ATAC-seq. and ChIP-seq data quality control, gene expression normalization, ASE identification, SNP detection, genotype imputation, \u003cem\u003ecis\u003c/em\u003e-eQTL and \u003cem\u003ecis\u003c/em\u003e-sQTL mapping, and functional annotation are available at the github website (https://github.com/Qiuming1986/ASE-in-cattle). All the original results of ASE, eQTL, sQTL, ATAC-seq, and Chip-seq have been deposited in the Figshare database with the following digital object identifier: 10.6084/m9.figshare.26808418.\u003c/p\u003e\n\u003cp\u003e\u003cstrong\u003eEthics approval and consent to participate\u003c/strong\u003e\u003c/p\u003e\n\u003cp\u003eAll animal procedures were conducted in accordance with the Regulations for the Administration of Affairs Concerning Experimental Animals of China and were approved by the Animal Care Committee of Xinjiang Agricultural University, which oversees the ethical use of animals in research at the university.\u003c/p\u003e\n\u003cp\u003e\u003cstrong\u003eCompeting interests\u003c/strong\u003e\u003c/p\u003e\n\u003cp\u003eThe authors declare that they have no competing interests.\u003c/p\u003e"},{"header":"References","content":"\u003col\u003e\u003cli\u003e\u003cspan\u003eWard LD, Kellis M. Interpreting noncoding genetic variation in complex traits and human disease. Nature biotechnology. 2012;30(11):1095\u0026ndash;106.\u003c/span\u003e\u003c/li\u003e \u003cli\u003e\u003cspan\u003eMaurano MT, Humbert R, Rynes E, Thurman RE, Haugen E, Wang H, et al. Systematic localization of common disease-associated variation in regulatory DNA. Science. 2012;337(6099):1190\u0026ndash;5.\u003c/span\u003e\u003c/li\u003e \u003cli\u003e\u003cspan\u003eConsortium TG, Ardlie KG, DeLuca DS, Segr\u0026egrave; AV, Sullivan TJ, Young TR, et al. The Genotype-Tissue Expression (GTEx) pilot analysis: multitissue gene regulation in humans. Science. 2015;348(6235):648\u0026ndash;60.\u003c/span\u003e\u003c/li\u003e \u003cli\u003e\u003cspan\u003eLiu S, Gao Y, Canela-Xandri O, Wang S, Yu Y, Cai W, et al. A multi-tissue atlas of regulatory variants in cattle. Nat Genet. 2022;54(9):1438\u0026ndash;47. Epub 2022/08/12. doi: \u003cspan class=\"ExternalRef\"\u003e\u003cspan class=\"RefSource\"\u003e10.1038/s41588-022-01153-5\u003c/span\u003e\u003cspan address=\"10.1038/s41588-022-01153-5\" targettype=\"DOI\" class=\"RefTarget\"\u003e\u003c/span\u003e\u003c/span\u003e. PubMed PMID: 35953587; PubMed Central PMCID: PMCPMC7613894.\u003c/span\u003e\u003c/li\u003e \u003cli\u003e\u003cspan\u003eTeng J, Gao Y, Yin H, Bai Z, Liu S, Zeng H, et al. A compendium of genetic regulatory effects across pig tissues. Nature Genetics. 2024:1\u0026ndash;12.\u003c/span\u003e\u003c/li\u003e \u003cli\u003e\u003cspan\u003eKuznetsova A, Brockhoff PB, Christensen RHB. lmerTest package: tests in linear mixed effects models. Journal of statistical software. 2017;82(13).\u003c/span\u003e\u003c/li\u003e \u003cli\u003e\u003cspan\u003eMapel XM, Kadri NK, Leonard AS, He Q, Lloret-Villas A, Bhati M, et al. Molecular quantitative trait loci in reproductive tissues impact male fertility in cattle. nature communications. 2024;15(1):674.\u003c/span\u003e\u003c/li\u003e \u003cli\u003e\u003cspan\u003eZou J, Hormozdiari F, Jew B, Castel SE, Lappalainen T, Ernst J, et al. Leveraging allelic imbalance to refine fine-mapping for eQTL studies. PLoS genetics. 2019;15(12):e1008481.\u003c/span\u003e\u003c/li\u003e \u003cli\u003e\u003cspan\u003eCavalli M, Pan G, Nord H, Arzt EW, Wallerman O, Wadelius C. Allele-specific transcription factor binding in liver and cervix cells unveils many likely drivers of GWAS signals. Genomics. 2016;107(6):248\u0026ndash;54.\u003c/span\u003e\u003c/li\u003e \u003cli\u003e\u003cspan\u003eMagris G, Jurman I, Fornasiero A, Paparelli E, Schwope R, Marroni F, et al. The genomes of 204 Vitis vinifera accessions reveal the origin of European wine grapes. Nature Communications. 2021;12(1):7240.\u003c/span\u003e\u003c/li\u003e \u003cli\u003e\u003cspan\u003eCleary S, Seoighe C. Perspectives on allele-specific expression. Annual Review of Biomedical Data Science. 2021;4(1):101\u0026ndash;22.\u003c/span\u003e\u003c/li\u003e \u003cli\u003e\u003cspan\u003eCastel SE, Aguet F, Mohammadi P, Ardlie KG, Lappalainen T. A vast resource of allelic expression data spanning human tissues. Genome biology. 2020;21:1\u0026ndash;12.\u003c/span\u003e\u003c/li\u003e \u003cli\u003e\u003cspan\u003eDelbare SY, Clark AG. Allele-specific expression elucidates cis-regulatory logic. PLoS Genetics. 2018;14(11):e1007690.\u003c/span\u003e\u003c/li\u003e \u003cli\u003e\u003cspan\u003eMoore JE, Purcaro MJ, Pratt HE, Epstein CB, Shoresh N, Adrian J, et al. Expanded encyclopaedias of DNA elements in the human and mouse genomes. Nature. 2020;583(7818):699\u0026ndash;710.\u003c/span\u003e\u003c/li\u003e \u003cli\u003e\u003cspan\u003eClark EL, Archibald AL, Daetwyler HD, Groenen MA, Harrison PW, Houston RD, et al. From FAANG to fork: application of highly annotated genomes to improve farmed animal production. Genome Biology. 2020;21:1\u0026ndash;9.\u003c/span\u003e\u003c/li\u003e \u003cli\u003e\u003cspan\u003eQuan J, Yang M, Wang X, Cai G, Ding R, Zhuang Z, et al. Multi-omic characterization of allele-specific regulatory variation in hybrid pigs. Nature Communications. 2024;15(1):5587.\u003c/span\u003e\u003c/li\u003e \u003cli\u003e\u003cspan\u003eKundu K, Tardaguila M, Mann AL, Watt S, Ponstingl H, Vasquez L, et al. Genetic associations at regulatory phenotypes improve fine-mapping of causal variants for 12 immune-mediated diseases. Nature genetics. 2022;54(3):251\u0026ndash;62.\u003c/span\u003e\u003c/li\u003e \u003cli\u003e\u003cspan\u003eShao L, Xing F, Xu C, Zhang Q, Che J, Wang X, et al. Patterns of genome-wide allele-specific expression in hybrid rice and the implications on the genetic basis of heterosis. Proceedings of the National Academy of Sciences. 2019;116(12):5653-8.\u003c/span\u003e\u003c/li\u003e \u003cli\u003e\u003cspan\u003eZhan W, Cui L, Yang S, Zhang K, Zhang Y, Yang J. Natural variations of heterosis-related allele-specific expression genes in promoter regions lead to allele-specific expression in maize. BMC genomics. 2024;25(1):476.\u003c/span\u003e\u003c/li\u003e \u003cli\u003e\u003cspan\u003eCui L, Yang B, Pontikos N, Mott R, Huang L. ADDO: a comprehensive toolkit to detect, classify and visualize additive and non-additive quantitative trait loci. Bioinformatics. 2020;36(5):1517\u0026ndash;21.\u003c/span\u003e\u003c/li\u003e \u003cli\u003e\u003cspan\u003eConsortium TG, Aguet F, Anand S, Ardlie KG, Gabriel S, Getz GA, et al. The GTEx Consortium atlas of genetic regulatory effects across human tissues. Science. 2020;369(6509):1318\u0026ndash;30.\u003c/span\u003e\u003c/li\u003e \u003cli\u003e\u003cspan\u003eGutierrez-Arcelus M, Baglaenko Y, Arora J, Hannes S, Luo Y, Amariuta T, et al. Allele-specific expression changes dynamically during T cell activation in HLA and other autoimmune loci. Nature genetics. 2020;52(3):247\u0026ndash;53.\u003c/span\u003e\u003c/li\u003e \u003cli\u003e\u003cspan\u003eEyre-Walker A, Keightley PD. The distribution of fitness effects of new mutations. Nature Reviews Genetics. 2007;8(8):610\u0026ndash;8.\u003c/span\u003e\u003c/li\u003e \u003cli\u003e\u003cspan\u003eRaney BJ, Barber GP, Benet-Pag\u0026egrave;s A, Casper J, Clawson H, Cline MS, et al. The UCSC Genome Browser database: 2024 update. Nucleic Acids Research. 2024;52(D1):D1082-D8.\u003c/span\u003e\u003c/li\u003e \u003cli\u003e\u003cspan\u003eLu Z, Hofmeister BT, Vollmers C, DuBois RM, Schmitz RJ. Combining ATAC-seq with nuclei sorting for discovery of cis-regulatory regions in plant genomes. Nucleic acids research. 2017;45(6):e41-e.\u003c/span\u003e\u003c/li\u003e \u003cli\u003e\u003cspan\u003eShlyueva D, Stampfel G, Stark A. Transcriptional enhancers: from properties to genome-wide predictions. Nature Reviews Genetics. 2014;15(4):272\u0026ndash;86.\u003c/span\u003e\u003c/li\u003e \u003cli\u003e\u003cspan\u003eHu Z-L, Park CA, Reecy JM. Bringing the animal QTLdb and CorrDB into the future: meeting new challenges and providing updated services. Nucleic acids research. 2022;50(D1):D956-D61.\u003c/span\u003e\u003c/li\u003e \u003cli\u003e\u003cspan\u003eJiang J, Cole JB, Freebern E, Da Y, VanRaden PM, Ma L. Functional annotation and Bayesian fine-mapping reveals candidate genes for important agronomic traits in Holstein bulls. Communications biology. 2019;2(1):212.\u003c/span\u003e\u003c/li\u003e \u003cli\u003e\u003cspan\u003eSchaid DJ, Chen W, Larson NB. From genome-wide associations to candidate causal variants by statistical fine-mapping. Nature Reviews Genetics. 2018;19(8):491\u0026ndash;504.\u003c/span\u003e\u003c/li\u003e \u003cli\u003e\u003cspan\u003eYuan C, Tang L, Lopdell T, Petrov VA, Oget-Ebrad C, Moreira GCM, et al. An organism-wide ATAC-seq peak catalog for the bovine and its use to identify regulatory variants. Genome Research. 2023;33(10):1848\u0026ndash;64.\u003c/span\u003e\u003c/li\u003e \u003cli\u003e\u003cspan\u003eAlexandre PA, Naval-S\u0026aacute;nchez M, Menzies M, Nguyen LT, Porto-Neto LR, Fortes MR, et al. Chromatin accessibility and regulatory vocabulary across indicine cattle tissues. Genome biology. 2021;22:1\u0026ndash;20.\u003c/span\u003e\u003c/li\u003e \u003cli\u003e\u003cspan\u003eWang C, Chen C, Lei B, Qin S, Zhang Y, Li K, et al. Constructing eRNA-mediated gene regulatory networks to explore the genetic basis of muscle and fat-relevant traits in pigs. Genetics Selection Evolution. 2024;56(1):1\u0026ndash;21.\u003c/span\u003e\u003c/li\u003e \u003cli\u003e\u003cspan\u003eFlynn ED, Lappalainen T. Functional characterization of genetic variant effects on expression. Annual Review of Biomedical Data Science. 2022;5(1):119\u0026ndash;39.\u003c/span\u003e\u003c/li\u003e \u003cli\u003e\u003cspan\u003eXuan A, Song Y, Bu C, Chen P, El-Kassaby YA, Zhang D. Changes in DNA methylation in response to 6-benzylaminopurine affect allele-specific gene expression in Populus tomentosa. International Journal of Molecular Sciences. 2020;21(6):2117.\u003c/span\u003e\u003c/li\u003e \u003cli\u003e\u003cspan\u003eKim-Hellmuth S, Bechheim M, P\u0026uuml;tz B, Mohammadi P, N\u0026eacute;d\u0026eacute;lec Y, Giangreco N, et al. Genetic regulatory effects modified by immune activation contribute to autoimmune disease associations. Nature communications. 2017;8(1):1\u0026ndash;10.\u003c/span\u003e\u003c/li\u003e \u003cli\u003e\u003cspan\u003eCui L, Yang B, Xiao S, Gao J, Baud A, Graham D, et al. Dominance is common in mammals and is associated with trans-acting gene expression and alternative splicing. Genome Biology. 2023;24(1):215.\u003c/span\u003e\u003c/li\u003e \u003cli\u003e\u003cspan\u003eChen ZJ. Genomic and epigenetic insights into the molecular bases of heterosis. Nature Reviews Genetics. 2013;14(7):471\u0026ndash;82.\u003c/span\u003e\u003c/li\u003e \u003cli\u003e\u003cspan\u003eWakchaure R, Ganguly S, Praveen PK, Sharma S, Kumar A, Mahajan T, et al. Importance of heterosis in animals: a review. International Journal of Advanced Engineering Technology and Innovative Science. 2015;1(2):1\u0026ndash;5.\u003c/span\u003e\u003c/li\u003e \u003cli\u003e\u003cspan\u003eGetahun D, Alemneh T, Akeberegn D, Getabalew M, Zewdie D. Importance of hybrid vigor or heterosis for animal breeding. Biochemistry and Biotechnology Research. 2019;7:1\u0026ndash;4.\u003c/span\u003e\u003c/li\u003e \u003cli\u003e\u003cspan\u003eBirchler JA, Yao H, Chudalayandi S. Unraveling the genetic basis of hybrid vigor. Proceedings of the National Academy of Sciences. 2006;103(35):12957-8.\u003c/span\u003e\u003c/li\u003e \u003cli\u003e\u003cspan\u003eTeng J, Wang D, Zhao C, Zhang X, Chen Z, Liu J, et al. Longitudinal genome-wide association studies of milk production traits in Holstein cattle using whole-genome sequence data imputed from medium-density chip data. Journal of Dairy Science. 2023;106(4):2535\u0026ndash;50.\u003c/span\u003e\u003c/li\u003e \u003cli\u003e\u003cspan\u003eBernini F, Mancin E, Sartori C, Mantovani R, Vevey M, Blanchet V, et al. Genome-wide association studies for milk production traits in two autochthonous Aosta cattle breeds. Animal. 2024;18(10):101322.\u003c/span\u003e\u003c/li\u003e \u003cli\u003e\u003cspan\u003eReynolds EG, Lopdell T, Wang Y, Tiplady KM, Harland CS, Johnson TJ, et al. Non-additive QTL mapping of lactation traits in 124,000 cattle reveals novel recessive loci. Genetics Selection Evolution. 2022;54(1):5.\u003c/span\u003e\u003c/li\u003e \u003cli\u003e\u003cspan\u003eBisutti V, Pegolo S, Giannuzzi D, Mota L, Vanzin A, Toscano A, et al. The β-casein (CSN2) A2 allelic variant alters milk protein profile and slightly worsens coagulation properties in Holstein cows. Journal of Dairy Science. 2022;105(5):3794\u0026ndash;809.\u003c/span\u003e\u003c/li\u003e \u003cli\u003e\u003cspan\u003eMiluchov\u0026aacute; M, G\u0026aacute;bor M, Candr\u0026aacute;k J. The effect of the genotypes of the CSN2 gene on test-day milk yields in the Slovak Holstein cow. Agriculture. 2023;13(1):154.\u003c/span\u003e\u003c/li\u003e \u003cli\u003e\u003cspan\u003eLi H, Durbin R. Fast and accurate short read alignment with Burrows\u0026ndash;Wheeler transform. Bioinformatics. 2009;25(14):1754\u0026ndash;60.\u003c/span\u003e\u003c/li\u003e \u003cli\u003e\u003cspan\u003eHan B, Yuan Y, Li Y, Liu L, Sun D. Single nucleotide polymorphisms of NUCB2 and their genetic associations with milk production traits in dairy cows. Genes. 2019;10(6):449.\u003c/span\u003e\u003c/li\u003e \u003cli\u003e\u003cspan\u003eBolger AM, Lohse M, Usadel B. Trimmomatic: a flexible trimmer for Illumina sequence data. Bioinformatics. 2014;30(15):2114\u0026ndash;20.\u003c/span\u003e\u003c/li\u003e \u003cli\u003e\u003cspan\u003eShen W, Le S, Li Y, Hu F. SeqKit: a cross-platform and ultrafast toolkit for FASTA/Q file manipulation. PloS one. 2016;11(10):e0163962.\u003c/span\u003e\u003c/li\u003e \u003cli\u003e\u003cspan\u003eDobin A, Davis CA, Schlesinger F, Drenkow J, Zaleski C, Jha S, et al. STAR: ultrafast universal RNA-seq aligner. Bioinformatics. 2013;29(1):15\u0026ndash;21.\u003c/span\u003e\u003c/li\u003e \u003cli\u003e\u003cspan\u003eKim D, Paggi JM, Park C, Bennett C, Salzberg SL. Graph-based genome alignment and genotyping with HISAT2 and HISAT-genotype. Nature biotechnology. 2019;37(8):907\u0026ndash;15.\u003c/span\u003e\u003c/li\u003e \u003cli\u003e\u003cspan\u003eLi H, Handsaker B, Wysoker A, Fennell T, Ruan J, Homer N, et al. The sequence alignment/map format and SAMtools. bioinformatics. 2009;25(16):2078\u0026ndash;9.\u003c/span\u003e\u003c/li\u003e \u003cli\u003e\u003cspan\u003ePertea M, Pertea GM, Antonescu CM, Chang T-C, Mendell JT, Salzberg SL. StringTie enables improved reconstruction of a transcriptome from RNA-seq reads. Nature biotechnology. 2015;33(3):290\u0026ndash;5.\u003c/span\u003e\u003c/li\u003e \u003cli\u003e\u003cspan\u003eGalili T. dendextend: an R package for visualizing, adjusting and comparing trees of hierarchical clustering. Bioinformatics. 2015;31(22):3718\u0026ndash;20.\u003c/span\u003e\u003c/li\u003e \u003cli\u003e\u003cspan\u003eMcKenna A, Hanna M, Banks E, Sivachenko A, Cibulskis K, Kernytsky A, et al. The Genome Analysis Toolkit: a MapReduce framework for analyzing next-generation DNA sequencing data. Genome research. 2010;20(9):1297\u0026ndash;303.\u003c/span\u003e\u003c/li\u003e \u003cli\u003e\u003cspan\u003eWang K, Li M, Hakonarson H. ANNOVAR: functional annotation of genetic variants from high-throughput sequencing data. Nucleic acids research. 2010;38(16):e164-e.\u003c/span\u003e\u003c/li\u003e \u003cli\u003e\u003cspan\u003eDanecek P, Auton A, Abecasis G, Albers CA, Banks E, DePristo MA, et al. The variant call format and VCFtools. Bioinformatics. 2011;27(15):2156\u0026ndash;8.\u003c/span\u003e\u003c/li\u003e \u003cli\u003e\u003cspan\u003eBrowning SR, Browning BL. Rapid and accurate haplotype phasing and missing-data inference for whole-genome association studies by use of localized haplotype clustering. The American Journal of Human Genetics. 2007;81(5):1084\u0026ndash;97.\u003c/span\u003e\u003c/li\u003e \u003cli\u003e\u003cspan\u003eZhang Z, Wang A, Hu H, Wang L, Gong M, Yang Q, et al. The efficient phasing and imputation pipeline of low-coverage whole genome sequencing data using a high‐quality and publicly available reference panel in cattle. Animal Research One Health. 2023;1(1):4\u0026ndash;16.\u003c/span\u003e\u003c/li\u003e \u003cli\u003e\u003cspan\u003eDanecek P, Bonfield JK, Liddle J, Marshall J, Ohan V, Pollard MO, et al. Twelve years of SAMtools and BCFtools. Gigascience. 2021;10(2):giab008.\u003c/span\u003e\u003c/li\u003e \u003cli\u003e\u003cspan\u003ePurcell S, Neale B, Todd-Brown K, Thomas L, Ferreira MA, Bender D, et al. PLINK: a tool set for whole-genome association and population-based linkage analyses. The American journal of human genetics. 2007;81(3):559\u0026ndash;75.\u003c/span\u003e\u003c/li\u003e \u003cli\u003e\u003cspan\u003eTamura K, Stecher G, Kumar S. MEGA11: molecular evolutionary genetics analysis version 11. Molecular biology evolution. 2021;38(7):3022\u0026ndash;7.\u003c/span\u003e\u003c/li\u003e \u003cli\u003e\u003cspan\u003eLetunic I, Bork P. Interactive Tree of Life (iTOL) v6: recent updates to the phylogenetic tree display and annotation tool. Nucleic Acids Research. 2024:gkae268.\u003c/span\u003e\u003c/li\u003e \u003cli\u003e\u003cspan\u003eLiu Y, Liu X, Zheng Z, Ma T, Liu Y, Long H, et al. Genome-wide analysis of expression QTL (eQTL) and allele-specific expression (ASE) in pig muscle identifies candidate genes for meat quality traits. Genetics Selection Evolution. 2020;52:1\u0026ndash;11.\u003c/span\u003e\u003c/li\u003e \u003cli\u003e\u003cspan\u003eMohammadi P, Castel SE, Brown AA, Lappalainen T. Quantifying the regulatory effect size of cis-acting genetic variation using allelic fold change. Genome research. 2017;27(11):1872\u0026ndash;84.\u003c/span\u003e\u003c/li\u003e \u003cli\u003e\u003cspan\u003eSherman BT, Hao M, Qiu J, Jiao X, Baseler MW, Lane HC, et al. DAVID: a web server for functional enrichment analysis and functional annotation of gene lists (2021 update). Nucleic acids research. 2022;50(W1):W216-W21.\u003c/span\u003e\u003c/li\u003e \u003cli\u003e\u003cspan\u003eZhou Y, Zhang Z, Bao Z, Li H, Lyu Y, Zan Y, et al. Graph pangenome captures missing heritability and empowers tomato breeding. Nature. 2022;606(7914):527\u0026ndash;34.\u003c/span\u003e\u003c/li\u003e \u003cli\u003e\u003cspan\u003eStegle O, Parts L, Piipari M, Winn J, Durbin R. Using probabilistic estimation of expression residuals (PEER) to obtain increased power and interpretability of gene expression analyses. Nature protocols. 2012;7(3):500\u0026ndash;7.\u003c/span\u003e\u003c/li\u003e \u003cli\u003e\u003cspan\u003eDelaneau O, Ongen H, Brown AA, Fort A, Panousis NI, Dermitzakis ET. A complete tool set for molecular QTL discovery and analysis. Nature communications. 2017;8(1):15452.\u003c/span\u003e\u003c/li\u003e \u003cli\u003e\u003cspan\u003eLi YI, Knowles DA, Humphrey J, Barbeira AN, Dickinson SP, Im HK, et al. Annotation-free quantification of RNA splicing using LeafCutter. Nature genetics. 2018;50(1):151\u0026ndash;8.\u003c/span\u003e\u003c/li\u003e \u003cli\u003e\u003cspan\u003eWu D-D, Ding X-D, Wang S, W\u0026oacute;jcik JM, Zhang Y, Tokarska M, et al. Pervasive introgression facilitated domestication and adaptation in the Bos species complex. Nature ecology evolution. 2018;2(7):1139\u0026ndash;45.\u003c/span\u003e\u003c/li\u003e \u003cli\u003e\u003cspan\u003eKrueger F. Trim Galore!: A wrapper around Cutadapt and FastQC to consistently apply adapter and quality trimming to FastQ files, with extra functionality for RRBS data. Babraham Institute. 2015;\u003cspan class=\"ExternalRef\"\u003e\u003cspan class=\"RefSource\"\u003ehttps://github.com/FelixKrueger/TrimGalore\u003c/span\u003e\u003cspan address=\"https://github.com/FelixKrueger/TrimGalore\" targettype=\"URL\" class=\"RefTarget\"\u003e\u003c/span\u003e\u003c/span\u003e.\u003c/span\u003e\u003c/li\u003e \u003cli\u003e\u003cspan\u003eLangmead B, Salzberg SL. Fast gapped-read alignment with Bowtie 2. Nature methods. 2012;9(4):357\u0026ndash;9.\u003c/span\u003e\u003c/li\u003e \u003cli\u003e\u003cspan\u003eTarasov A, Vilella AJ, Cuppen E, Nijman IJ, Prins P. Sambamba: fast processing of NGS alignment formats. Bioinformatics. 2015;31(12):2032\u0026ndash;4.\u003c/span\u003e\u003c/li\u003e \u003cli\u003e\u003cspan\u003eZhang Y, Liu T, Meyer CA, Eeckhoute J, Johnson DS, Bernstein BE, et al. Model-based analysis of ChIP-Seq (MACS). Genome biology. 2008;9:1\u0026ndash;9.\u003c/span\u003e\u003c/li\u003e \u003cli\u003e\u003cspan\u003eQuinlan AR, Hall IM. BEDTools: a flexible suite of utilities for comparing genomic features. Bioinformatics. 2010;26(6):841\u0026ndash;2.\u003c/span\u003e\u003c/li\u003e \u003cli\u003e\u003cspan\u003eProwse-Wilkins CP, Wang J, Xiang R, Garner JB, Goddard ME, Chamberlain AJ. Putative causal variants are enriched in annotated functional regions from six bovine tissues. Frontiers in genetics. 2021;12:664379.\u003c/span\u003e\u003c/li\u003e \u003cli\u003e\u003cspan\u003eProwse-Wilkins CP, Lopdell TJ, Xiang R, Vander Jagt CJ, Littlejohn MD, Chamberlain AJ, et al. Genetic variation in histone modifications and gene expression identifies regulatory variants in the mammary gland of cattle. BMC genomics. 2022;23(1):815.\u003c/span\u003e\u003c/li\u003e \u003cli\u003e\u003cspan\u003eProwse-Wilkins CP, Wang J, Garner JB, Goddard ME, Chamberlain AJ. Allele specific binding of histone modifications and a transcription factor does not predict allele specific expression in correlated ChIP-seq peak-exon pairs. Scientific Reports. 2023;13(1):15596.\u003c/span\u003e\u003c/li\u003e \u003cli\u003e\u003cspan\u003eVillar D, Berthelot C, Aldridge S, Rayner TF, Lukk M, Pignatelli M, et al. Enhancer evolution across 20 mammalian species. Cell. 2015;160(3):554\u0026ndash;66.\u003c/span\u003e\u003c/li\u003e\u003c/ol\u003e"}],"fulltextSource":"","fullText":"","funders":[],"hasAdminPriorityOnWorkflow":false,"hasManuscriptDocX":true,"hasOptedInToPreprint":true,"hasPassedJournalQc":"","hasAnyPriority":false,"hideJournal":false,"highlight":"","institution":"","isAcceptedByJournal":true,"isAuthorSuppliedPdf":false,"isDeskRejected":"","isHiddenFromSearch":false,"isInQc":false,"isInWorkflow":false,"isPdf":false,"isPdfUpToDate":true,"isWithdrawnOrRetracted":false,"journal":{"display":true,"email":"[email protected]","identity":"bmc-biology","isNatureJournal":false,"hasQc":true,"allowDirectSubmit":false,"externalIdentity":"","sideBox":"Learn more about [BMC Biology](https://bmcbiol.biomedcentral.com/)","snPcode":"12915","submissionUrl":"https://submission.springernature.com/new-submission/12915/3","title":"BMC Biology","twitterHandle":"","acdcEnabled":true,"dfaEnabled":true,"editorialSystem":"stoa","reportingPortfolio":"BMC Series","inReviewEnabled":true,"inReviewRevisionsEnabled":true},"keywords":"","lastPublishedDoi":"10.21203/rs.3.rs-5530951/v1","lastPublishedDoiUrl":"https://doi.org/10.21203/rs.3.rs-5530951/v1","license":{"name":"CC BY 4.0","url":"https://creativecommons.org/licenses/by/4.0/"},"manuscriptAbstract":"\u003cp\u003e\u003cstrong\u003eBackground\u003c/strong\u003e\u003c/p\u003e\n\u003cp\u003eAllele-specific expression (ASE) analysis is a crucial tool for validating expression quantitative trait loci (eQTLs), identifying causal variants associated with complex traits, and investigating the genetic mechanisms underlying heterosis. In this study, we characterized ASE variants across 35 tissues using 7,532 publicly available RNA-seq datasets. Additionally, we explored the mechanisms driving ASE through integration with epigenomic data and examined the relationship between ASE and dominance effects on gene expression and milk-related traits in Holstein cattle.\u003c/p\u003e\n\u003cp\u003e\u003cstrong\u003eResults\u003c/strong\u003e\u003c/p\u003e\n\u003cp\u003eASE variants exhibited stronger tissue specificity and lower reproducibility compared to eQTLs. Interestingly, variants with opposite directional effects demonstrated greater resilience across diverse environments. Functional annotation revealed that ASE variants were predominantly located in enhancer regions during transcription, rather than promoter regions. Furthermore, ASE variants were implicated in post-transcriptional and translational processes, including mutations affecting mRNA splicing and triggering nonsense-mediated decay. Analysis of eQTLs, splicing QTLs (sQTLs), and validated QTLs associated with milk-related traits in Holstein cattle, coupled with enrichment analysis in QTL databases and effect size evaluation, indicated that ASE variants were more closely aligned with dominant effects than additive effects, particularly in reproductive and immune-related tissues/traits, which exhibited higher levels of heterosis.\u003c/p\u003e\n\u003cp\u003e\u003cstrong\u003eConclusions\u003c/strong\u003e\u003c/p\u003e\n\u003cp\u003eOur findings not only enhance our understanding of the genetic mechanisms underlying heterosis and ASE formation but also provide a valuable resource of regulatory variants that can be leveraged to improve economic traits through molecular breeding or the strategic exploitation of heterosis.\u003c/p\u003e","manuscriptTitle":"A multi-tissue atlas of allelic-specific expression reveals the characteristics, mechanisms, and relationship with dominant effects in cattle","msid":"","msnumber":"","nonDraftVersions":[{"code":1,"date":"2025-01-09 08:17:46","doi":"10.21203/rs.3.rs-5530951/v1","editorialEvents":[{"type":"communityComments","content":0},{"type":"decision","content":"Revision requested","date":"2025-02-12T23:49:23+00:00","index":"","fulltext":""},{"type":"editorInvitedReview","content":"","date":"2025-02-02T07:54:01+00:00","index":"hide","fulltext":""},{"type":"editorInvitedReview","content":"","date":"2025-01-31T07:56:21+00:00","index":"hide","fulltext":""},{"type":"reviewerAgreed","content":"127416212843488670167388464173046249417","date":"2025-01-25T02:36:13+00:00","index":"hide","fulltext":""},{"type":"reviewerAgreed","content":"178420038241610038633466231013536868962","date":"2025-01-25T00:12:50+00:00","index":"hide","fulltext":""},{"type":"editorInvitedReview","content":"","date":"2025-01-02T09:53:04+00:00","index":"hide","fulltext":""},{"type":"reviewerAgreed","content":"183787388704084829668280668692572621199","date":"2024-12-12T15:23:40+00:00","index":"hide","fulltext":""},{"type":"reviewersInvited","content":"","date":"2024-12-11T17:00:34+00:00","index":"","fulltext":""},{"type":"editorAssigned","content":"","date":"2024-11-27T16:27:58+00:00","index":"","fulltext":""},{"type":"checksComplete","content":"","date":"2024-11-27T09:12:36+00:00","index":"","fulltext":""},{"type":"submitted","content":"BMC Biology","date":"2024-11-26T23:47:18+00:00","index":"","fulltext":""}],"status":"published","journal":{"display":true,"email":"[email protected]","identity":"bmc-biology","isNatureJournal":false,"hasQc":true,"allowDirectSubmit":false,"externalIdentity":"","sideBox":"Learn more about [BMC Biology](https://bmcbiol.biomedcentral.com/)","snPcode":"12915","submissionUrl":"https://submission.springernature.com/new-submission/12915/3","title":"BMC Biology","twitterHandle":"","acdcEnabled":true,"dfaEnabled":true,"editorialSystem":"stoa","reportingPortfolio":"BMC Series","inReviewEnabled":true,"inReviewRevisionsEnabled":true}}],"origin":"","ownerIdentity":"0e22485f-97c0-4dc8-85e5-a8570bc5e1f9","owner":[],"postedDate":"January 9th, 2025","published":true,"recentEditorialEvents":[],"rejectedJournal":[],"revision":"","amendment":"","status":"under-review","subjectAreas":[],"tags":[],"updatedAt":"2025-05-21T13:38:45+00:00","versionOfRecord":[],"versionCreatedAt":"2025-01-09 08:17:46","video":"","vorDoi":"","vorDoiUrl":"","workflowStages":[]},"version":"v1","identity":"rs-5530951","journalConfig":"researchsquare"},"__N_SSP":true},"page":"/article/[identity]/[[...version]]","query":{"redirect":"/article/rs-5530951","identity":"rs-5530951","version":["v1"]},"buildId":"XKTyCvWXoU3ODBz1xrDgd","isFallback":false,"isExperimentalCompile":false,"dynamicIds":[84888],"gssp":true,"scriptLoader":[]}