3.1. Expression QTL mapping
3.1.1. Expression cis and DE cis eQTL analysis
LSM2 binds to other members of the ubiquitous and multifunctional family Sm-like (LSM) in order to form RNA-processing complexes. These complexes are involved in processes such as stabilization of the spliceosomal U6 snRNA, mRNA decay and guide site-specific pseudouridylation of rRNA [5]. The hetero-heptameric complex LSM2-8 promotes small stable RNAs and pre-mRNAs nuclear processing meanwhile the LSM1-7 complex stimulates mRNA decapping and decay in the cytoplasm. The LSM2 subunit seems to be crucial for determining cellular location of this complex [6]. Lu et al. [7] identified two missense polymorphisms in SOAT1 associated with plasma cholesterol and triglyceride levels in mice. These SNPs increase enzyme activity despite similar gene expression levels. Although SOAT1 has been associated with cholesterol and triglyceride levels in mice, expression of this gene was not identified in the DE analysis in the present research.
3.1.2. Expression trans and DE trans eQTL analysis, and master regulators
The 27 master regulators identified in the eQTL analysis could contribute to gene expression control by promoting cell signaling or by direct transcriptional activation/repression mechanisms. The most important master regulators are described below.
Neurotrophin 3 (NTF3) was identified as master regulator in the present analysis since rs207649022 was able to explain variation in the expression of 76 genes (Table 1), 69.7% of which were DE genes (Figure 3B). The Neurotrophic Factor gene family regulates myoblast and muscle fiber differentiation. It also coordinates muscle innervation and functional differentiation of neuromuscular junctions [8]. Mice with only one functional copy of the NTF3 gene showed smaller cross-sectional fiber area and more densely distributed muscle fibers [9]. Upregulation of NTF3, stimulated by the transcription factor POU3F2, is present during neuronal differentiation [10]. The neocortex has multiple layers originated by cell fate restriction of cortical progenitors and NTF3 induces cell fate switches by controlling a feedback signal between postmitotic neurons and progenitors. Therefore, changes in NTF3 expression can modulate the amount of tissue present in each neocortex layer [11]. This gene was identified in a previous study as highly associated with cooking loss [12] pointing out that markers inside this locus are able to explain variation in both the phenotypic meat quality and the gene expression associated with meat quality. This provides positional and functional support for the potential role of NTF3 on meat quality. These effects seem not to be due to cis regulation on NTF3 giving that this gene was not expressed in skeletal muscle in the present population. However, NTF3 could be involved in earlier expression regulation based on its biological function. A Functional Annotation Clustering analysis for the NTF3 cluster indicated that the master regulator NTF3 could be involved in regulation of specific mechanisms and pathways related to Mitochondrion, Transit peptide and Mitochondrion inner membrane (Additional file 4).
The Phosphodiesterase 8B (PDE8B) was expressed in longissimus dorsi muscle in the present population. A SNP in this gene was found associated with expression levels of 27 genes, 74% of which were identified as DE genes related to meat quality. PDE8B hydrolyzes adenosine 3',5'-cyclic monophosphate (cAMP), and it is involved in a number of signal transduction pathways and physiological processes such as cell growth, cell differentiation, transcription and expression [13]. The cAMP-degrading enzymes have tight subcellular localization, and control local cAMP and signal compartmentalization. Zhang et al. [14] reported that the phosphorylated cAMP-response Element Binding Protein (CREB) binds to approximately 4,000 promoter sites in vivo, depending on the methylation level of cAMP response elements near the promoter which highlights the role CREB plays in regulation of expression for a multitude of genes. CREB has a crucial role in neuronal membrane-to-nucleus signal transduction, regulate bone and cartilage remodeling by regulating the Matrix Metalloproteinase 1 (MMP1) and 13 (MMP13), stimulates promoter activity of adiponectin gene in adipocytes and inhibits the Cholesterol 7 α-Hydroxylase (CYP7A1) gene, which is involved in lipid homeostasis and bile acid synthesis in hepatocytes [15–17].
Expression of 36 genes was associated with rs211476449, a marker located in the Glutamate Decarboxylase 1 (GAD1) master regulator. This gene encodes a plasma membrane associated protein and downregulation of GAD1 has a detrimental effect on the cortical c-aminobutyric acid (GABA) system. This system has a paramount function on cellular proliferation and differentiation in multiple tissues [18].
FAT Atypical Cadherin 4 (FAT4) encodes a calcium-dependent cell adhesion transmembrane protein and it is involved in the Hippo signaling pathway. This pathway regulates organ size and tissue organization in vertebrates; when this gene is disrupted it alters oriented cell divisions. Inactivation of FAT4 promotes tumorigenesis in mammary epithelial cells and tumor progression is suppressed by FAT4 re-expression [19,20]. FAT is also involved in mammalian neurogenesis since downregulation of FAT in embryonic neuroepithelium in mice reduces the proportion of differentiated neurons [21]. Expression of 34 genes was associated with BTB_00676236, a marker located close to the master regulator FAT4.
The expression of 26 genes was associated with rs208227436, a polymorphism located in the master regulator Zinc Finger Protein 804A (ZNF804A). One polymorphism in the nuclear transmembrane protein ZNF804A gene is associated with schizophrenia and bipolar disorder; however, no function was identified for this gene. Since this protein has a zinc-finger domain, it was suggested that ZNF804A is involved in DNA binding and transcription regulation [22]. A direct interaction was reported between ZNF804A and the chromatin proximal to the promoter regions of several genes, some of them associated with cell adhesion and the Transforming Growth Factor-β (TGF-β) signaling pathway, probably influencing cell growth and differentiation.
The ENSBTAG00000035487 gene encodes a homologous of the Olfactory Receptor 4A47 (OR4A47) and the SNP rs109630111, was associated with expression of 36 genes. During perception of smell, olfactory receptors trigger signal transduction pathways by recognizing odorant molecules and these pathways might regulate apoptotic cycle of olfactory sensory neurons. Because these genes are expressed in non-olfactory related tissues, it is suggested that OR4A47 could have different biological functions [23]. Neuhaus et al. [24] and Ranzani et al [23] found expression of the olfactory receptor OR2C3 in human melanomas and documented that activation of OR51E2, another member of the olfactory receptor family, increases intracellular Ca2+, triggers some Mitogen-Activated Protein Kinases (MAPK) and inhibits cell proliferation.
Another master regulator, Paired Box 8 (PAX8) encodes a transcription factor involved in thyroid cell differentiation, and kidney and gonadal development. During kidney and kidney related tumor development, both transcription factors WT1 and PAX8 are co-expressed and since WT1 promoter has one PAX-binding site with enhancer activity, PAX8 is suggested as able to drive WT1 expression [25]. The plasma membrane associated protein Pleckstrin and Sec7 Domain Containing 4 (PSD4) activates ARF6, regulating its interaction with specific plasma membrane subdomains and controlling reorganization of the actin cytoskeletal structure [26]. A total of 37 genes belong to the PAX8 cluster by being associated with rs209448226.
The SNP rs208451702 is located in RUNX1 Translocation Partner 1 (RUNX1T1 or Myeloid Translocation Gene on 8q22-MTG8) and explains expression variability in 24 other loci. It has transcriptional corepressor activity by recruiting other corepressors, and by interacting with DNA-binding transcription factors and histone-modifying enzymes. MTGs were recognized in acute myelogenous leukemia and mutations in this gene are associated with colon, breast and lung carcinoma, and have negative effects on the WNT and Notch signaling [27].
The rs377935001 marker is located in Tricopeptide Repeat Domain 25 (TTC25) gene and it was associated with expression of 34 genes. Deficient TTC25 mice lack the cilium outer dynein arms and its docking complex. The cilium is involved in motility and sensory-related processes such as left-right axis patterning, having a crucial developmental and homeostatic role [28,29].
The cytoskeleton associated protein Keratin 7 (KRT7) cluster shares some trans regulated genes with CSAD and was associated with expression of other 25 genes. The KRT7 is part of a keratin superfamily; KRT7 knockout mice have upregulated proliferation in bladder urothelium, KRT18 downregulation in bladder and KRT20 upregulation in superficial urothelial cells [30].
Lysine Demethylase 4A (KDM4A) cluster has 32 regulated genes associated with rs135786834; KDM4A encodes a histone lysine demethylase able to modify trimethylated H3-K9/K36 to dimethylated products, contributing to gene expression, cellular differentiation and cancer development [31]. Histone H3K9 methylation is used for silencing muscle specific genes in proliferating myoblasts and their derepression is required to initiate muscle differentiation; expression of a KDM4A isoform named DN-JMJD2A is upregulated during differentiation of myoblasts into myotubes promoting myotube formation and transcriptionally activating muscle-specific genes such as MyoD [32]. Hu et al. [33] reported that IOX1, an KDM4A inhibitor, is able to stall proliferation, migration and cell cycle progression by modulating cyclin D1 and p21 expression in angiotensin II stimulated vascular smooth muscle cells. This stalling process is mediated by promoter methylation. The only DE master regulator identified in the present analysis was KDM4A and this master regulator harbors rs135786834, a SNP associated with expression of 32 genes by trans association. Therefore, changes in expression of the master regulator KDM4A did not show evidence of promoting expression associated with meat quality inside its cluster.
Expression of 62 genes was associated with rs378343630, a marker located in the Transmembrane 4 L Six Family Member 1 (TM4SF1) master regulator. This gene encodes a plasma transmembrane protein and belongs to a gene family involved in signal transduction processes; thus, it modulates development, growth and motility [34]. The TM4SF1 protein physically interacts with membrane and some cytoskeleton-associated proteins to form cell projections named ‘nanopodia’ [35], which are described as frequently identified in multiple types of cancer. This gene is highly expressed in pancreatic cancer cells and stimulates metastasis by upregulating the Discoidin Domain Receptor Tyrosine Kinase 1 (DDR1) gene, Matrix Metallopeptidase 2 (MMP2) and Matrix Metallopeptidase 9 (MMP9) [36]. In liver, TM4SF1 reduced apoptosis and promoted cell migration by upregulating MMP-2, MMP-9 and VEGF, and downregulating Caspase-3 and Caspase-9 [34]. Upregulation of miR-9 produces downregulation of TM4SF1, MMP2, MMP9 and VEGF in colorectal carcinoma inhibiting cell migration and invasion [37]. In esophageal cancer stem-like cells, downregulation of miR-141 increases TM4SF1 expression, self-renewal ability and carcinogenicity, and promotes cell invasion and migration [38]. The Functional Annotation Clustering analysis for TM4SF1 found overrepresentation of the transcription, DNA-templated term (Additional file 4); thus, TM4SF1 could be involved in regulation of specific mechanisms and pathways associated with transcription in longissimus dorsi muscle.
The master regulator GPR98 was identified as associated with tenderness from the sensory panel in the present population by Leal-Gutiérrez et al. (2019) and rs110618957, a polymorphism harbored by this gene, is able to explain variability in expression of 34 genes. GPR98 encodes a receptor associated with some bone related features in human and mice [39].
Pierce et al. [2] theorized that a high proportion of trans associations were caused by cis effects. However, no cis QTL was identified in any expression or splicing master regulator. The previous result seems to reveal that in the present analysis trans effects might contribute more to phenotypic variation related to meat quality than cis effects.
3.1.3. Multigenic effects based on the eQTL analysis
A total of 126 markers were identified as able to explain variation in SLC43A1 gene expression. This gene encodes a Na+-independent neutral amino acid transporter and it is upregulated in prostate cancer and hepatocarcinoma cells [40,41]. Forty three SNPs were associated with ULK2 expression. ULK2 is a serine/threonine protein kinase required for autophagy under nutrient-deprived conditions [42] and downregulation of ULK2 activates mTOR c1 signaling, promoting cell proliferation rates [43]. The MYL1 gene encodes a fast-twitch regulatory light chain of myosin in skeletal muscle; downregulation of MYL1 alters myocyte morphology and muscle structure, and generates congenital myopathy in zebrafish [44]; expression of MYL1 exclusively starts in fast-twitch cells during fast fibre differentiation [45]. A total of 40 and 36 polymorphisms were associated with MYL1 and PHF14 expression, respectively. PHF14 is ubiquitously expressed and its protein has multiple PHD fingers, a domain present in chromatin-binding proteins which is able to recognize particular epigenetic marks on histone tails. Knockout for PHF14 generates neonatal lethality and severe structural changes in multiple organs especially lungs in mice, being PHF14 an epigenetic regulator required for normal development of multiple organs and associated with some types of cancer [46,47]. Thirty six SNPs were associated with ENO3 expression and upregulation of ENO3 was evident in metastatic liver and muscle tissue [48]. These genes might have a tight expression regulation since multiple genomic regions are able to explain their expression.
3.2. Splicing QTL mapping
3.2.1. Splicing cis and DE cis analysis
The TTN gene harbors a DE cis sQTL (p-value = 2.0x10-7) and encodes a central sarcomeric protein. Some TTN mutations are associated with skeletal-muscle diseases such as tibial muscular dystrophy [49]. Fernandez-Marmiesse et al. [50] identified a non-sense mutation in a TTN exon only present in a fetal skeletal isoform and associated with a neuromuscular disorder; histologically, this mutation promotes sarcomeric deposition of a filamentous material.
A DE cis sQTL (p-value = 5.1x10-7) was identified in the TEK gene. The tyrosine kinase TEK encodes a receptor for Angiopoietin-1 (ANGPT1), and this signaling pathway is critical for migration, sprouting and survival of endothelial cells; TEK activates the SHC Adaptor Protein 1 (SHC1), a protein involved in triggering the Ras/mitogen-activated protein kinase pathway, regulating migration and endothelial organization induced by ANGPT1 [51]. An ANGPT1-TEK antagonist named Angiopoietin-2 (ANGPT2) is expressed in regions of vascular remodeling in mice and humans and its upregulation is able to stall blood vessel formation in mouse embryos [52].
3.2.2. Trans and DE trans splicing QTL analysis, and master regulators
The non-coding RNAs, Small Nucleolar RNA, H/ACA Box 3A (SNORA3) and Small Nucleolar RNA, H/ACA Box 19 (SNORA19) are small nucleolar RNA molecules (snoRNA); snoRNAs modulate stability, folding and interaction with proteins and more recently, functions such as mRNA editing, alternative splicing and posttranscriptional gene silencing were discovered [53]. However, no clear function about SNORA3 and SNORA19 was described. Exon expression of 34 exons from 17 genes and 20 exons from 15 genes were associated with the polymorphisms rs209617551 (SNORA3) and BTB_01634267 (SNORA19), respectively.
Expression of 33 exons from 21 genes was associated with the SNP rs381222773, located in Phosphodiesterase 9A (PDE9A). This gene encodes a metal ion-dependent enzyme. The subcellular location of PDE9A is the plasma membrane in almost all organs but in bladder, where it is a cytoplasmic protein. Several PDE9A isoforms are described in human and mouse brain, and their expression and subcellular compartmentalization are age and tissue dependent [54]. PDE9A regulates nuclear- and membrane-proximal cGMP, regulating cellular signaling, and being involved in multiple biological and metabolic processes. It is also associated with several human neurological disorders by modulating intraneuronal signal transmission [55]. Upregulation of PDE9A is noticed in hypertrophic cardiomyocytes and during heart failure, but this symptomology is attenuated by PDE9A downregulation [56].
Expression of 23 exons form 19 genes was associated with rs382101207, a SNP located in Ring Finger Protein 20 (RNF20). This protein stimulates global H2B ubiquitylation at lysine 120 and promotes activator-dependent transcription [57]. Upregulation of RNF20 stimulates H2B monoubiquitination and methylation at H3K4 and H3K79; it promotes expression of Homeobox genes, a group of transcription factors [58]. RNF20 regulates expression of H2A and H2B histones, p53, several protooncogenes, epidermal growth factor and promotes cell migration and tumorigenesis [59]. The RNF20/RNF20 (Bre1 complex) is documented as a tumor suppressor by upregulating a set of tumor suppressor genes and by contributing to genomic stability maintenance. Bre1 deficient cells presents a high frequency of DNA double-strand breaks (DSB), having abundant aberrant RNA-DNA structures (R-loops) as indicators of replication stress and genomic instability [58]. H2B ubiquitylation is required during embryonic stem cell differentiation, process induced by stimulation of RNF20 activity and stalling of the deubiquitinase USP44 [60].
3.2.3. Multigenic effects based on the sQTL analysis
The large number of sQTLs identified in genes like TTN (324) and NEB (63) could be related to gene size, since these genes were 275 and 219 kb long, respectively, which would increase the probability of being involved in trans regulation. On the other hand, some genes relatively small in size: TCEB2 (9.9 kb) and USF2 (3.9 kb) also had a large number of sQTLs (43 and 33, respectively) indicating a complex splicing regulation.
A total of 324 and 67 polymorphisms were associated with TTN and NEB ratio exon counts, respectively. TTN and NEB are involved in assembly and mechanical activity of striated muscles. Both proteins are large sarcomere filament-binding proteins uniformly expressed in skeletal muscle and multiple splicing events in the bovine homologous have been described; in human brain, NEB acts as an actin filament stabilizer, regulates neuronal length, it is involved in myofibrillogenesis, modulates thin filament length and allows proper muscle contraction [61]. NEB deficient mice present muscular weakness, they die 8-11 days after birth, and show altered regulation of calcium homeostasis and glycogen metabolism associated genes [61]. Mouse NEB has 166 exons and four alternative splicing regions, two of them in exons 127 and 128; in most muscles, young mice had higher expression of the isoform including the exon 127 while the isoform including the exon 128 was highly expressed in older mice, showing that there could exist isoform switching related to muscle maturation [62].
Thirty-three sQTLs were identified in CAMP Responsive Element Binding Protein 5 (CREB5); CREB5 is upregulated in several types of human cancer since it modulates cell cycle. Upregulation of CREB5 promotes liver, colorectal and epithelial ovarian cancer cell proliferation and is an indicative of poor prognosis.
A total of 33 sQTLs were identified for the Upstream Transcription Factor 2, C-Fos Interacting (USF2) gene. Cell cycle control and modulation of apoptosis are some of the cellular processes where USF2 is involved. Christensen et al. [63] and Jaiswal & Narayan [64] reported that miR-362-3p modulate expression of USF2, a transcription factor able to interact with E-boxes located in the promoter region of Adenomatous Polyposis Coli (APC); APC is able to promote cell cycle arrest and apoptosis, being its malfunctioning associated with colorectal carcinogenesis. DNA-binding assays for a polymorfism associated with thyroid cancer susceptibility located in the 5’UTR of FOXE1, showed that the alternative allele is able to recruit the transcription factor USF2 [65].
3.3. Gene expression and splicing regulation mechanisms by plasma and organelle associated proteins
Cell cytoskeleton provides cellular mechanical constraints and extracellular matrix stiffness [66]. However, structural proteins are involved in multiple biological processes different from the organizational ones, with signaling and cell fate being some of the most important ones. Cell signaling is crucial since it orchestrates cellular responses to different microenvironmental stimulus, and transcription repression-activation and splicing regulation are influenced by signaling proteins. A number of receptors, transmembrane linkers, cytoskeletal fibers and membrane associated transcription factors were previously associated with transcription repression-activation.
The OR4A47, GPR98, PDE9A, OR13F1 and SYT14 master regulators are also described as transmembrane proteins and this type of molecule is involved in cell signaling processes. Pandey et al. [67] reported that estrogen can signal using diverse receptors, the G Protein-Coupled Estrogen Receptor 1 (GPR30) being one of them. Stimulation of GPR30 by estrogen activates a transcription factor network that upregulates Cellular Communication Network Factor 2 (CCN2), promoting proliferation and cell migration. The master regulators GAD1 and TM4SF1 encode transmembrane linkers and the integrin family has some well described transmembrane proteins. Integrins can modulate signal transduction cascades involved in cell survival, proliferation, differentiation and organ development [68]. The dimer ITGA1-ITGB1 can stall Epidermal Growth Factor Receptor (EGFR) signaling by stimulating Protein Tyrosine Phosphatase, Non-Receptor Type 2 (PTPN2). The cytoplasmic domain of ITGA1 interacts with PTPN2 and decreases EGFR phosphorylation after Epidermal Growth Factor (EGF) stimulation [69].
The cytoskeletal proteins KRT7, FAT4, MYH14 and DNAH7 were identified as master regulators, showing that cytoskeletal proteins might drive transcription regulation promoting cellular mechanisms such as growth and apoptosis. Flouriot et al. [66] reported that actin network can regulate Myocardin Related Transcription Factor A (MRTFA) subcellular localization, a protein involved in growth-quiescence switch. High F/G actin ratio or mutant MRTFA cells showed higher global biosynthetic activity and open chromatin state, associated with extensive histone modifications. In Drosophila, Hippo tumor suppressor pathway controls organ size, and proteins such as Yorkie (human homologous Yes Associated Protein 1 -YAP), a transcriptional coactivator, and Hpo and Warts kinases (human homologous Serine/Threonine Kinase 3 -STK3- and Large Tumor Suppressor Kinase 1 -LATS1, respectively) belong to this pathway. YAP is negatively regulated by STK3 and LATS1. F-actin accumulation promotes overgrowth in Drosophila imaginal discs through modulating the activity of the Hippo pathway [70]. Zhao et al. [71] reported that cell detachment promoted by cytoskeleton reorganization induces phosphorylation and inactivation of YAP, a process required for anoikis, a type of anchorage-dependent apoptosis.
The expression master regulator ZNF804A can be classified as a membrane associated transcription factor. A better described member of this group is CAMP Responsive Element Binding Protein 3 Like 1 (CREB3L1). This gene encodes an endoplasmic reticulum transmembrane stress transducer with a basic leucine zipper (bZIP) transcription factor domain. As a response to endoplasmic reticulum stress in osteoblasts, CREB3L1 is subject of regulated intramembrane proteolysis and its cleaved bZIP domain translocates into the nucleus, where initiates gene transcription [72,73]. The promotor region of the Collagen Type I Alpha 1 Chain (COL1A1) has an unfolded protein response element (UPRE)-like sequence in its promoter region, allowing COL1A1 expression controlled by the membrane associated protein CREB3L1 [73].