Genome variants altering gene splicing in bovine are extensively shared between tissues

Background Mammalian phenotypes are shaped by numerous genome variants, many of which may regulate gene transcription or exon splicing. We aim to identify splicing quantitative trait nucleotides (sQTN) in bovine, an important economic and model species, then compare sQTN with well-reported gene expression geQTN and lesser-reported exon expression eeQTN in different tissues. Results Using whole genome and RNA sequence data from four tissues of over 200 cattle, sQTN identified using exon inclusion ratios were verified by matching their effects on adjacent intron excision ratios. sQTN had the highest percentage of intronic variants and the lowest percentage of intergenic variants, compared to eeQTN and geQTN. For sQTN and eeQTN, there were more cases where the same variant was detected (FDR<0.01) in multiple tissues, compared to geQTN. To verify such expression QTN sharing between tissues, variants surrounding (±1Mb) the exon or gene were used to build local genomic relationship matrices and estimated genetic correlations between tissues. For many exons, the splicing and expression level was determined by the same cis additive genetic variance in different tissues. Thus, an effective meta-analysis combining information from three tissues increased power to detect and validate sQTN. sQTN and eeQTN together were more enriched for variants associated with cattle complex traits, compared to geQTN. Several potentially causal mutations were identified, including an sQTN at Chr6:87392580 within the 5th exon of kappa casein (CSN3) affecting milk production traits. Conclusions We have substantially improved the knowledge of genome regulation for an important mammalian species and provided novel insights into the analytics of expression QTN.


Background
Cattle are an important source of meat and dairy products for humans worldwide. Also, cattle can be used as clinical models to study genetic causes of human diseases [1]. To improve productivity, health performance and efficiency of cattle, traditional selective breeding has been widely used. In the last decade, the genomic selection technique, originally developed in cattle breeding, has further increased the rate of genetic improvement of complex traits of all livestock species [2,3]. However, genomic selection commonly uses random single nucleotide polymorphisms (SNPs) and, with a few exceptions, we do not know the genes or the sites causing genetic variation in important traits. Knowledge of the genes involved and polymorphic sites would increase our understanding of the biology and may further increase the rate of genetic improvement [4].
Many of the sequence variants that affect complex traits (quantitative trait nucleotides or QTN) are not coding variants and are presumed to influence the regulation of gene expression, that is to be expression QTN (eQTN) [5]. A regulatory mutation might affect expression of the gene in general, which commonly refers to a gene eQTN or geQTN. In humans, geQTN show significant enrichments for mutations associated with diseases and complex traits [5,6].
A regulatory mutation might change the expression of one particular or several exons within a gene, which can be defined as an exon eQTN or eeQTN. Furthermore, after transcription, RNAs are spliced by intron removal and exon ligation to create various mature transcripts. A mutation altering the RNA splicing by changing the balance between various transcripts from the gene, is then defined as a splicing QTN or sQTN. If a mutation changes the expression ratio of an exon to the gene, it implies that it alters the splicing pattern. In humans, sQTN were studied by inferring the individual splicing ratio from RNA sequence data [7]. More recently, sQTN identified using intron information extracted from RNA sequence data were demonstrated to have fundamental links with human diseases [8,9].
Despite vastly reported human regulatory QTN information, the knowledge of large animal regulatory mutations is minimal. In this study, we aim to map bovine cis splicing QTN using counts of genes, exons and introns of RNA sequence data from hundreds of animals and multiple tissues along with imputed whole genome sequences. We examined the extent to which sQTN affected splicing in white blood cells, milk cells, liver and muscle transcriptomes and the extent to which sQTN are enriched for variants affecting complex traits. We also used the counts of genes and exons to map cis eeQTN and geQTN, then analysed their relationships with sQTN in different tissues.

Results
In total, we analysed 378 transcriptomes of 19 tissue types from 214 cattle generated from four experiments covering major dairy and beef cattle breeds (Table 1, Figure 1a and Additional file 1: Supplementary Methods). Overall, samples from the same or similar tissues clustered together rather than clustering by experiments, based on exon expression levels ( Figure 1a). Consistent with previous reports [10], milk cells and mammary gland transcriptomes were closely related.

Differential splicing between tissues and breeds
We primarily defined a differentially spliced gene as a gene which contained exons whose inclusion ratios (exon expression divided by gene expression) were significantly affected by tissue or breeds (FDR<0.1). To verify the significantly spliced exons, we imposed a requirement that at least one adjacent intron had an excision ratio [8,9] that was also significantly (FDR<0.1) affected by tissue or breeds (See methods). The overall FDR of such exon splicing events was considered as 0.01 as the exon and intron analyses were independent. The agreement of genes displaying differential splicing from exon and intron analyses were shown and examined in 2 × 2 tables by Chi-square tests (Additional file 2: Supplementary Table S1). Overall, the results from exon and intron analyses were significantly overlapped.
Using data from all experiments (Table 1), there were 8,657 genes with at least one differential spliced exon caused by differences between tissue types. A list of these genes with the significances of differential splicing for the exons and introns were shown in Additional file 3: Supplementary Table S2. The top 10% of these significantly differentially spliced genes had a GO term enrichment (FDR<0.01) of 'regulation of cellular process', suggesting very general roles of these genes in cell functionalities. A much smaller number of genes (148) had differential splicing events in the milk cell transcriptome between breeds (  Figure 1b).

cis splicing quantitative trait nucleotides (sQTN)
The sQTN mapping was based on data from 312 transcriptomes data generated from experiments III and IV, covering four major tissues (Table 1). In total 207 individuals had imputed whole genome sequence data and in experiment III, 105 genotyped cattle had both white blood and milk cell transcriptome data (Table 1). Similar to differential splicing analyses described above, a cis-acting sQTN was defined as a SNP significantly (FDR<0.1) affecting the inclusion ratio of the exon (up to 1Mb away) and significantly (FDR<0.1) affecting at least one adjacent introns' excision ratio [8,9] . When analysed separately, the sQTN found by exon and intron analyses were also significantly overlapped (Additional file 2: Supplementary The significant sQTN in the white blood and milk cells were mapped to 929 and 283 genes, respectively ( Table 2). For instance, in the milk cell transcriptome, the fifth exon of CSN3 (6:87392578-87392750), which as described above was differentially spliced in Holstein and Jersey cattle (Figure 1b), was significantly affected by a sQTN (Chr6:87392580, Additional file 5: Supplementary Table S4). The B allele of this SNP increased the expression of the 5 th exon and had a higher allele frequency among Holstein cattle than Jerseys cattle (0.79 vs 0.02), which explains the difference between the breeds. This sQTN was also predicted to be a splice site by Ensemble [14]. Much smaller numbers of significant sQTN were detected in liver (11,544 SNPs) and muscle (5,783 SNPs) (Figure 2c,d). This was probably due to the smaller sample size (Table 1) and lower sequence depth of liver and muscle from the experiment IV (Additional file 1: Supplementary Methods) than that of blood and milk from the experiment III ( Figure 2a,b, Table 2).

Shared genetic influences between cis QTN types
Within each tissue, the sharing of SNPs between all three eQTN types was significantly more than expected by chance ( Figure 4). However, in the white blood and milk cells, which had relatively large sample size (n>=105, Table 1), the largest absolute amount of SNP sharing appeared to be between sQTN and eeQTN ( Figure 4). This was followed by the amount of SNP sharing between eeQTN and geQTN ( Figure 4a,b). In liver and muscle tissue which had relatively small sample size (n<=41), the largest absolute amount of SNP sharing was between eeQTN and geQTN, followed by the amount of SNP sharing between sQTN and eeQTN ( Figure 4c,d).
To further examine the relationship between sQTN and eeQTN, a two by two

Shared genetic influences between tissues
Within each QTN type between different tissues, the majority of the significant QTN sharing were observed between white blood and milk cells and between liver and muscle ( Figure 5a).
This not unexpected since the most of the white blood and milk cells came from the same lactating cows and the muscle and liver from different growing Angus bulls. The largest amount of cross-tissue QTN sharing was observed in eeQTN, followed by sQTN and geQTN ( Figure 5). Where a SNP significantly affected expression in two tissues, the direction of effect was usually the same in both tissues (Additional file 10: Supplementary Figure S1).
The correlation of QTN effects between white blood cells and milk cells (Additional file 10: Supplementary Figure S1a,c,e) was stronger than that between liver and muscle (Additional file 10: Supplementary Figure S1b,d,f). The SNP level sharing between white blood cells and milk cells and between liver and muscle were also evident at the exon (splicing quantitative trait loci sQTL and exon expression eeQTL) and gene (geQTL) level ( Figure 5b).
The QTN sharing between tissues was further examined for all types of eQTN by using a less stringent p-value (p<0.05) to test their effect (Additional file 10: Supplementary Figure S2).
This showed that the QTN sharing between tissues were stronger for sQTN and eeQTN Often both splicing events and exon expression within a gene were highly correlated between white blood and milk cells, for instance DDX19B, CTSD and EFF1A1 ( Figure 5c). In liver and muscle, exons from HLA-DQA1 encoding major histocompatibility complexes [15] also showed significant genetic correlations between tissues based on both exon expression and splicing. There were more cases of eeQTN than sQTN and geQTN and so there were more estimates of genetic correlations between white blood and milk cells in Figure 5c. The genetic correlations between eeQTN in white blood and milk cells show a range from +1 to -1 although most are close to +1. Exons with negative genetic correlations of expression between white blood cells and milk cells were mapped to SART1 [15], a post-transcriptional regulator in epithelial tissues and TTC4 with potential to mediate protein-protein interactions [15]. These negative genetic correlations imply that there are mutations that increase the expression of the exon in milk cells but decrease it in white blood cells. An exon within S100A10, a cell cycle progress regulator, showed negative genetic correlation of expression between liver and muscle.
Genetic correlations between exon expression in two tissues can be different between exons within the same gene. For example, only the 2 nd exon (5:93,942,055 -93,942,195) of MGST1 (which affects dairy cattle milk fat yield [13,16]) had a significant genetic correlation of expression between white blood cells and milk cells (Figure 5c). This was largely due to a few eeQTN with very strong effects on the expression levels of the 2 nd exon in the milk cells and a similar but smaller effect on the 2 nd exon expression in white blood cells (Additional file 10: Supplementary Figure S3).

Multi-transcriptome meta-analysis to increase power of expression QTN detection
Based on shared genetic effects of all types of eQTN across tissues, a multi-transcriptome meta-analysis was introduced to increase the power to detect sQTN, eeQTN and geQTN ( Figure 6, Table 3). For sQTN, eeQTN and geQTN that had significant effects (p < 0.05) in all of white blood cells, milk cells and muscle transcriptomes, their standardised effects (signed t values) in each transcriptome were simply combined and tested for significance against a χ 2 distribution with 1 degree of freedom. Overall, the multi-transcriptome metaanalysis based on summary statistics substantially increased the power of QTN detection ( Figure 6). The significance of multi-transcriptome QTN were compared with their significance in the liver transcriptome ( Figure 6, Table 3). For a criteria where the expression QTN had both multi-transcriptome meta-analysis p < 1×10 -5 and liver transcriptome analysis p < 0.05, all types of expression QTN had significant overlaps of the SNPs between the metaanalysis and the single transcriptome analysis in liver. In fact most of the significant sQTN, eeQTN and geQTN detected by the meta-analysis were also detected by the liver analysis but at a much higher p-value (Table 3).

Overlap between QTN and QTL for dairy and beef traits
We examined whether cis sQTN, eeQTN and geQTN were significantly enriched amongst SNPs associated with cattle economic traits. Pleiotropic SNPs significantly (FDR<0.01) associated with more than one of 24 dairy traits [17] and with more than one of 16 beef traits [18] were used to overlap with detected sQTN, eeQTN and geQTN ( present on the high density SNP chip data we analysed for dairy cattle pleiotropy [17] (Additional file 9: Supplementary Table S8). However, 53 significant milk cell eeQTN identified by the current study overlapped with the top 200 SNPs from Littlejohn et al [16] (Additional file 9: Supplementary Table S8), which was significantly more than expected by random chance. The overlapping milk cell eeQTN included the SNP suggested as a potential causal candidate (Chr5:93945738) [16], which significantly affected the expression level of the third exon (5:93939150-93939244) of MGST1 in milk cells (Additional file 9: Supplementary

Discussion
We performed a systematic analysis of transcriptomic cis QTN (<=1Mb) in multiple tissues centred around RNA splicing events, using a large number of RNA and whole genome sequence data from an important domestic animal species. Overall, differential splicing between tissues is ubiquitous and between breeds is common. Differential splicing between individuals due to SNPs (sQTN) occurs for many genes and are enriched with cattle complex trait QTL. Within each tissue, all cis QTN types showed significant overlap. However, SNPs causing an sQTN are also likely to cause eeQTN and, to a lesser extent, geQTN. Between tissues, while all QTN types showed significant overlap between white blood cells and milk cells and between liver and muscle, the strongest cross-tissue sharing appeared to be at the exon level (sQTN and eeQTN). This is supported by many significant tissue pair genetic correlations. Such cross-tissue QTN sharing allowed the multi-transcriptome meta-analysis of QTN effects which substantially increases power to detect significant QTN.
The majority of significant sQTN were detected from white blood and milk cells (Figure 2a,b) which also overlapped with SNP chip based complex trait QTL (Figure 7), compared to sQTN detected from liver and muscle. This is probably due to the larger sample size for white blood and milk cells than for liver and muscle (Table 1) Table S6) and this SNP is associated with shear force in multiple taurine breeds [19]. This implies that a splicing QTN for CAPN1 may contribute to beef tenderness variation.
In the milk cell transcriptome, a significant sQTN (Chr6:87392580, Figure 2a) with predicted splicing function [14] within the fifth exon (6:87,392,578-87,392,750) of CSN3 causes differential splicing between Holstein and Jerseys (Figure 1b). Variants within CSN3 have long been found to be associated with milk traits [20,21] but only recently have potential causal variants been prioritised [12]. The milk cell sQTN 6:87392580 had a perfect linkage disequilibrium (r =1) with the variant 6:87390576 which has been suggested as a possible causal variant for effects on milk protein yield and percentage [12,13] (Additional file 7: Supplementary Table S6). Given it is at a splicing site, it is likely that 6:87392580 is a causal variant contributing to milk production in dairy cattle by altering exon splicing.
Compared to identified bovine cis geQTN, cis sQTN tended to be closer to the transcription starting site (TSS) and had highest concentrations of intronic SNPs (Figure 3). In humans, cis sQTN [8,22] also had stronger enrichments of intron SNPs than other types of QTN.
However, reports of the distance between human QTN and TSS appear be inconsistent.
While no difference in enrichment of SNPs near TSS between sQTN and geQTN were found by the human GTEx project [7], A more recent study [8] found that human geQTN were more enriched near TSS than sQTN. Our results appear to stand in between the results of GTEx project and the later findings from Li et. al. [8], where cattle sQTN were slightly more closer to TSS compared to geQTN. However, this difference is not significant in all tissues ( Figure   3a). On the other hand, overall significant overlap between sQTN and geQTN was found in this study ( Figure 4) and by the human GTEx project [7]. However, Li et. al. [8] found that human cis sQTN were independent of geQTN. These inconsistent observations were likely to be due to a number of differences between these studies, including definition of sQTN, choice of tissues and populations and computational procedures. Also, these inconsistent observations also suggest that we are still at the very early stage of understanding of sQTN features.
Within each studied bovine tissue, the largest amount of SNP overlap between QTN types were found either between exon expression eeQTN and sQTN or between eeQTN and geQTN ( Figure 4). Further, the largest amount of enrichments of cattle pleiotropic SNPs were found for eeQTN, followed by sQTN and geQTN. The white blood cell eeQTN showed particularly strong enrichments of pleiotropic SNPs for dairy and beef cattle. In a large scale human blood cell QTN study [23], eeQTN also showed the strongest enrichments of GWAS variants, followed by sQTN and geQTN. Thus, focusing on exon-level QTN, including eeQTN and sQTN, could increase the chance of finding regulatory variants for complex traits, as proposed by Guan et. al. [24].
A hypothesis to explain these results is that mutations in regulatory DNA may increase the expression of one or more transcripts from a gene. If they increase expression of one transcript then they may be detected that as an eeQTN for the exons in that transcript, as an sQTN for exons spliced out of that transcript or as a geQTN if this transcript forms a large part of the total transcription from the gene. Thus, there is expected to be overlap between eeQTN, geQTN and sQTN, but at least sQTN and eeQTN should overlap and this is what we found ( Figure 4, Additional file 6: Supplementary Table S5). It appears that eeQTN detect the largest proportion of these regulatory polymorphisms provided sequencing depth is high.
In humans, significant cross-tissue sharing of sQTN was reported [7] and so was geQTN in different studies [7,25]. In our study of cattle, the strongest evidence of QTN sharing appeared to be at the exon level. This included sQTN and eeQTN sharing between white blood and milk cells and between liver and muscle ( Figure 5). When extending the examination of QTN to include those with small effects (p < 0.05, Additional file 10: Supplementary Figure S2), the exon-level QTN cross-tissue sharing is also the greatest. We Based on the sharing of QTN between tissues, a multi-transcriptome meta-analysis which simply combined QTN effects to substantially increase the power ( Figure 6) was introduced.
Using this approach, combined QTN effects of white blood cells, milk cells and muscle were validated in the liver (Table 3). This also demonstrated the significant extent of QTN sharing across tissues. Previously, Flutre et. al. [25] combined data from human fibroblasts, lymphoblastoid cell lines and T-cells and found that up to 88% of geQTN were shared across tissues at FDR<0.05 level. In the current study, when the FDR threshold was lowered to 0.05, the meta-analysis combining SNP effects from tissues of white blood cells, milk cells and

Methods
Sample collection. For Experiment I, the sampling of 18 tissues from one lactating Holstein cow followed procedures described by Chamberlain et al [11]. For Experiment II and III, the sampling and processing of all tissues including white blood and milk cells was detailed semitendinosus muscle and 35 liver from Angus bulls was previously described by [27,28].

RNA seq data. For Experiment I, RNA extraction and sequencing was followed
Chamberlain et al [11]. For Experiment II and III, the RNA extraction and sequencing procedures was detailed in Additional file 1: Supplementary Methods. For Experiment IV, RNA extraction and sequencing is previously described by Khansefid et al. [28]. For all experiments sequence quality was checked and were aligned to the Ensembl UMD3.1 bovine genome assembly using TopHat2 [26]. The RNA sequence data processing and results were detailed in Additional file 1: Supplementary Methods.
Whole genome seq data. Experiments III and IV had whole genome sequence genotypes imputed from the SNP chip genotypes using FImpute [29] based on the 1000 bull genomes project [30]. 50K Illumina genotypes were used for imputation for experiment III with previous protocols [13,31]. For experiment IV, 800K and 50K Illumina genotypes were used with procedures following [18]. SNPs were filtered for minor allele frequency > 0.01 and resulted in 14,302,604 and 13,632,145 SNPs used in the analysis in experiment III and IV, respectively. There were 10,242,837 SNPs shared between experiment III and IV.
Gene/exon analysis. Gene count data were generated by Python package HTSeq [32] using default settings. The exon count data were generated by Bioconductor package featureCounts [33] in R v3.3.2 [34]. Genes and exons with count per million >0 in more than 40% individuals in each tissue type were used for all following analyses. The exon-based tissue principal components analysis used DEseq2 [35]. The exon inclusion as calculated as the exon to gene expression ratio and the intron excision was calculated as the intron to global intron cluster expression ratio, processed by leafcutter package [8,9]. Those exons and introns with ratio values <0.001 were removed and the remaining ratio values were transformed to log 2 scale, then underwent exon/intron -wise quantile normalisation and individual -wise zscore standardisation [36].
Gene differential splicing. Both exon inclusions and intron excisions were analysed and used in combination for gene differential splicing for (1) the overall tissue effects and (2) the breed effects. Primarily, differential splicing was defined for the gene containing exons whose inclusion ratios were significantly (FDR<0.1) affected by the tissue or breed variable.
To be called as significantly spliced exons, they were required to have at least one adjacent intron whose excision ratios were also significantly (FDR<0.1) affected by the tissue or breed variable. The tissue effects were analysed in a linear mixed model in lme4 [37] in R as: Where y = exon inclusion or intron excision ratios, b i = the animal random effects (i=214), x j = the experiments (j=4), t k = tissue type (k=19), e= random residual terms. The breed effects for the milk transcriptome data were analysed in a linear model in R as: Where y = exon inclusion or intron excision ratios in the milk cell transcriptome, breed l = breeds (l=2, Holstein and Jersey). p values of exons/introns for the tissue effects in equation (1) or the breed effects in equation (2) were used to calculate the false discovery rates (FDR) using qvalue [38] in R. As the exon and intron analyses were relatively independent, the overall FDR threshold of such detected exon/intron group was considered as 0.1× 01 = 0.01.
cis expression splicing QTN. Similar to differential splicing analysis described above, a significant (FDR<0.01) cis splicing QTN was demanded to satisfy two conditions simultaneously: (1) a SNP, within or up to ± 1Mb away from the exon, significantly (FDR<0.1) affected its inclusion ratio and (2) the same SNP significantly (FDR<0.1) affected at least one event of the excision of the intron next to the same exon at the same significance level. Both individual exon inclusion and intron excision values were used as phenotype to map associated QTN with widely used [7] Matrix eQTL [39] package in R. For each cell type of the experiment III (white blood and milk cells) and experiment IV (liver and muscle), SNPs ± 1Mb from the exon or intron were tested for regressions with the exon inclusion or intron excision phenotype. For milk cell transcriptome, breed was fitted as a covariate.
To compare cis sQTN with exon expression cis eeQTN and cis gene expression geQTN, the expression count data were normalised by voom [40] estimating mean-variance relationship to calculate observation-level weighted expression values. Normalised expression values of exons and genes were used as phenotype to map cis QTN (within ± 1Mb) at FDR <0.01 level as described above.
SNP annotation. The gene transcription start site coordinates were downloaded from Ensembl (http://www.ensembl.org) and the absolute difference between the position of a SNP and the transcription start site of the gene were calculated for SNP with significant cis effects.
The SNP functional categories were generated using predictions from Ensembl Variant Effect Predictor [14] in conjunction with NGS-SNP [41]. All analysed SNPs were assigned a functional category.

Dairy and beef cattle pleiotropic QTL. To test the significance of overlap between cis QTN
and SNPs associated with cattle phenotype, meta-analyses of dairy and beef cattle pleiotropy were performed using single-trait GWAS results from Xiang et al [17] and Bolormaa et al [18]. HD 800K SNP chip genotype were used for trait GWAS. 24 dairy cattle traits with matching phenotype in 9,662 bulls and cows and 16 beef cattle traits with animal numbers >2,000 were selected. Briefly, the multi-trait χ 2 statistic for the ith SNP was calculated based on its signed t values generated from each single trait GWAS [18]: For dairy cattle, the meta-analysis was based on the weighted SNP effects t w combining SNP effects calculated separately in bulls and cows. The t w accounting for the phenotypic error differences between bulls and cows [13] was calculated as: Where the weighted SNP t value t w was the quotient of the weighted SNP effects B w and the weighted effect error se w . B bull and se bull were the SNP effects and error obtained from singletrait GWAS in bulls and B cow and se cow were the SNP effects and error of cows. Those SNPs which had meta-analysis FDR<0.01 were chosen to be compared with cis QTN. The lead SNP loci were defined as ±1 Mb from the lead SNPs identified in the previous analysis [17] and [18].
The significances of overlap were compared with the expected number using the Fisher's exact test (p) implemented in GeneOverlap [42] in R. The union number of whole genome sequence SNPs with MAF>0.01 in each breed and the bovine high density chip SNPs was used as the background. If QTN categories from different breeds of dairy and beef cattle were tested for the overlap, the number of common SNPs between breeds were used.
The cross-tissue sharing of SNPs were confirmed by bivariate GREML analysis using GCTA [43]. For an exon or a gene of interest, its inclusion ratios or expression levels in two different tissues were treated as two different phenotype, tr 1 and tr 2 . The SNPs within 1Mb of this exon or gene were used to make a genomic relationship matrix representing the local polygenic component a. This allowed linear mixed modelling of the additive genetic variances of tr 1 , var g (tr 1 ) and of tr 2 , var g (tr 2 ) and the additive genetic covariance between t 1 and t 2 , cov g (tr 1 ,tr 2 ) using GREML [43]. This then allowed the estimation of genetic  (5). Genetic correlations were also tested for their significances of being different from 0 or 1, by fixing the correlation value to 0 and 1 using GCTA [43].

Declarations Ethics approval and consent to participate
Not applicable.

Consent for publication
Not applicable.

Availability of data and materials
The RNA sequence data for experiment I was published [11] (NCBI Sequence Read Archive, SRA, accession SRP042639); For experiment II: SRA accession SRP111067; For experiment III: SRA accession PRJNA305942; For experiment IV: SRA accession PRJNA392196; The imputed whole genome sequence data was part of the 1000 bull genome project [30].

Competing interests
The authors declare that they have no competing interests.   and muscle tissue (d). A significant sQTN was defined as a SNP affecting the exon inclusion ratio and also at least one excision of an adjacent intron at the same significance level. The input SNPs had significance p<0.0001. sQTN in all tissues with their associated genes and significance were given in Supplementary Table S4.