Characterization of novel lncRNA muscle expression profiles associated with meat quality in beef cattle

Abstract The aim of this study was to identify novel lncRNA differentially expressed (DE) between divergent animals for beef tenderness and marbling traits in Nellore cattle. Longissimus thoracis muscle samples from the 20 most extreme bulls (of 80 bulls set) for tenderness, tender (n = 10) and tough (n = 10) groups, and marbling trait, high (n = 10) and low (n = 10) groups were used to perform transcriptomic analysis using RNA‐Sequencing. For tenderness, 29 lncRNA were DE (p‐value ≤ 0.01) in tough beef animals in relation to tender beef animals. We observed that genic lncRNAs, for example, lncRNA_595.1, were overlapping exonic part of the PICK gene, while lncRNA_3097.2 and lncRNA_3129.5 overlapped intronic part of the genes GADL1 and PSMD6. The lncRNA associated with PICK1, GADL1, and PMD6 genes were enriched in the pathways associated with the ionotropic glutamate receptor, gamma‐aminobutyric acid synthesis, and the ubiquitin–proteasome pathway. For marbling, 50 lncRNA were DE (p‐value ≤ 0.01) in high marbling group compared with low marbling animals. The genic lncRNAs, such as lncRNA_3191.1, were overlapped exonic part of the ITGAL gene, and the lncRNA_512.1, lncRNA_3721.1, and lncRNA_41.4 overlapped intronic parts of the KRAS and MASP1 genes. The KRAS and ITGAL genes were enriched in pathways associated with integrin signaling, which is involved in intracellular signals in response to the extracellular matrix, including cell form, mobility, and mediates progression through the cell cycle. In addition, the lincRNAs identified to marbling trait were associated with several genes related to calcium binding, muscle hypertrophy, skeletal muscle, lipase, and oxidative stress response pathways that seem to play a role important in the physiological processes related to meat quality. These findings bring new insights to better understand the biology mechanisms involved in the gene regulation of these traits, which will be valuable for a further investigation of the interactions between lncRNA and mRNAs, and of how these interactions may affect meat quality traits.


| INTRODUC TI ON
Meat quality is defined by the product nutritional composition and by a group of subjective consumer perceptions as visual appearance, smell, firmness, juiciness, tenderness, and flavor (FAO, 2014;Maltin et al., 2003). Results from meat consumption studies indicate that quality is even more important than the product price to consumers (Henchion et al., 2014). Among the many factors that affect meat quality, the main one reported by consumers is tenderness (Fletcher, 2002), and another major factor is marbling, which affects both flavor and juiciness (Williams, 2008). These are complex traits, controlled by many genes and greatly influenced by environmental factors (Leal-Gutiérrez & Mateescu, 2019). Moreover, these traits are expensive and difficult to measure, making the application of traditional selection, based-phenotype and pedigree, even more difficult, since it requires the slaughter of the animals, leading to an increase in the generation interval and a decrease in genetic gain .
An approach that can collaborate to increase the knowledge about the genetic regulation of these traits is transcriptomics. The understanding of the individual genetic mechanisms behind transcript expression profiles involved in the variation of meat quality traits (e.g., tenderness and marbling) are still unclear, and further studies about regulatory elements which control gene expression of complex traits are necessary. The long noncoding RNA (lncRNAs) have a significant role in a wide variety of important biological processes, such as gene expression regulation and control of translation or genomic imprinting (Wucher et al., 2017). The lncRNAs are transcripts longer than 200 nucleotides and without any protein-coding capabilities (Etebari et al., 2015;Gupta et al., 2019). Over the last decade, several lncRNAs have been identified and characterized, and this became possible because of the implementation of whole transcriptome sequencing technologies (Yao et al., 2019). However, few studies on the overall expression patterns of lncRNAs in Longissimus thoracis muscle have been conducted (Sun et al., 2016;Zhang et al., 2020), and they are even scarcer for meat quality traits. One of those few studies reported on this field was performed by Jiang et al. (2020) that investigated the transcriptome profiling of lncRNA related to fat tissues in Qinchuan cattle. Another study, by , investigated the molecular and expression characteristics of a novel lncRNA (lncFAM200B) along with its crucial genetic variations. They reported that the lncFAM200B expression trend is positively correlated with MyoG and Myf5 expression in myoblast proliferation. However, no studies on lncRNA expression profile for meat tenderness in beef cattle have been reported. There is still much to discover about the lncRNA functions involved in cellular regulation pathways in meat quality traits. Thus, the objective of this study was to identify novel long genic and intergenic noncoding RNA, differentially expressed in L. thoracis muscle of Nellore cattle divergent for tenderness and marbling phenotypes. These results would provide new insights about the lncRNA function within a context of meat quality traits and enhance scientific data for further investigation in beef cattle breeding.

| Samples set and RNA-seq data preparation
Animals were selected from a population of 80 Nellore bulls belonging to the Capivara farm located in Sao Paulo state, Brazil, which participates in the Nellore Qualitas breeding program and belonged to the same contemporary group (Muniz et al., 2021). Animals were slaughtered with an average age of 24 months, in a commercial slaughterhouse. Samples of the L. thoracis muscle were collected between the 12th and 13th ribs of the left half carcass. Two samples were obtained at two times: one sample at the time of slaughter, for RNA extraction, and another sample 24 h after slaughter for meat quality evaluation by methodologies described by Fonseca et al. (2017, Santos et al. (2020), and Muniz et al.
(2020) that used the same dataset.
Tenderness was evaluated by Warner-Bratzler shear force (WBSF), following the methodology proposed by Bratzler (1949), using the mechanical salter WBSF device. The visual marbling scores were evaluated according to USDA Quality Grade methodology (2000), which uses a marbling classification scale with values ranging from 1 to 9, being 1 value applied to absent marbling (devoid) and 9 corresponding to excessive marbling (moderately abundant).
However, in Brazil, this scale ranges from zero (devoid) to six (moderate), because of the low degree of marbling in the Brazilian herd Fonseca, Suárez-Vega, et al., 2020).
According to the phenotypes scores obtained, the samples were sorted, and the 20 most extreme animals, for each trait, were chosen, composing two groups of 10 animals each (HIGH (n = 10) and LOW (n = 10) values) for marbling score and tenderness. The Student's t-test was used to verify the average difference between both groups. The low-marbling grade groups had an average of 2.26 ± 0.05 and high-marbling grade groups had an average of 3.26 ± 0.12. In addition, the tender meat (LOW group) had an average of 4.11 ± 0.30 (kgf), and the tough meat groups (HIGH group) interactions between lncRNA and mRNAs, and of how these interactions may affect meat quality traits.

K E Y W O R D S
genic lncRNA, lincRNA, marbling, RNA-Seq, tenderness had an average of 8.93 ± 1.23 (kgf). The groups were statistically different by the Student's t-test (p-value < 0.001).
The RNA extraction was performed by methodologies described by Fonseca et al. (2017), Santos et al. (2020), and Muniz et al. (2020) that used the same dataset. The integrity of the RNA samples was assessed in the Agilent 2100 Bioanalyzer (Agilent, 2009), and the RNA concentration and genomic DNA contamination were determined using the Qubit ® 2.0 Fluorometer (Invitrogen, 2010) (Fonseca et al., 2017). RNA integrity number values were higher than 7.0 for all muscle samples, indicating good RNA quality. Paired end (2 × 100 bp) reads were sequenced using the HiSeq 2500 sequencer (Illumina).

| RNA-seq expression analyses
For the initial steps of the analysis, quality control of reads was performed based on Phred score, and GC content and over-represented sequence parameters, using the CLC Genomics Workbench software 20.0.4 (CLC Bio), following the parameters described in Cánovas et al. (2014).

| Large gap mapping and transcript discovery
The "Large Gap Read Mapping" tool, implemented in CLC Genomics Workbench software 20.0.4 (CLC Bio), was used to map the paired-end sequence reads according to the annotated reference genome ARS.UCD1.2 (ftp://ftp.ensem bl.org/pub/relea se-95/genba nk/ bos_tauru s/) following the parameters: length fraction and similarity fraction = 0.8; two mismatches; and three insertions and three deletions per read were allowed. The Large Gap Read Mapping tool maps reads to a reference, allowing for large gaps in the mapping. It is developed to support transcript discovery using RNA-seq data, since it is able to map RNA-seq reads that span introns without requiring prior transcript annotations, for more details see Muniz et al. (2021).
CLC Genomic workbench's Transcript Discovery plugin (CLC Bio; Release 20.4) was used for transcript discovery in bovine based on genome reference (ARS.UCD1.2; ftp.ensembl.org/pub/release-100/ fasta/bos_taurus/). This tool has an algorithm that takes large gap read mapped as input, here each trait was analyzed separately, tenderness and marbling phenotypes (i.e., one group of high and other group of low for each trait). For reads mapping, the following criteria were allowed: mismatch (2), insertion (3), and deletion (3) costs; were allowed a minimum similarity (0.8) and length fraction (0.9) between a mapped segment and the reference; and a gap with maximum 50 Kbp distance between mapped read segments to span the introns from RNA-Seq data was considered (Etebari et al., 2015). In addition, splice events were setting for: minimum length of ORF = 100; ignore duplicate or nonspecific read matches; minimum spliced and unspliced reads = 10 and minimum spliced and unspliced coverage ratio = 0.05 were allowed; chimeric or unknown splice signatures were ignored, gene merging distance = 50 bp; minimum reads in gene = 10; and minimum predicted gene length = 200 (Muniz et al., 2021). Then, a single assembly (GTF file) was generated for each trait.

| Long noncoding RNA (lncRNA) identification and reads mapping
To identify potential long noncoding RNAs from average of 887,304 predicted transcript (gtf file) models for each trait, Feelnc filter pipeline was used (Wucher et al., 2017). The transcripts with a length shorter than 200 bp, biotype-coding protein, single-exon transcripts, and biexonic transcripts having one exon size shorter than 25 bp were discarded (Etebari et al., 2015;Gupta et al., 2019). After filtering, the Feelnc coding potential tool, setting up with default parameters, was used to compute the coding potential score of the model transcripts, classifying them into putative lncRNAs and protein-coding RNAs. Then, the lncRNA transcript file, containing 5903 transcripts was merged with the bovine genome reference (ARS.UCD1.2; ftp.ensembl.org/pub/release-100/ fasta/bos_taurus/), which was used to create gene and RNA track making.
The RNA-Sequencing analysis tool (CLC Genomics Workbench environment) was used to align the sequences of each sample against the predicted gene and transcript tracks, using the bovine reference genome ARS.UCD1.2 (ftp://ftp.ensem bl.org/pub/relea se-100) as map, see more information about assembly parameters in Muniz et al. (2021).

| Differential lncRNA expression
Differential expression analyses were performed using the CLC Genomics Workbench software 12.0 (CLC Bio). A two-group experiments [tender vs. tough tenderness and high vs. low marbling groups] were carried out from each trait. Empirical analysis of DGE tool implemented in CLC workbench was performed for each set up experiment, using the original count values and two parameters related to the estimation of the dispersion were specified: (1) "total count filter cut-off" >5, that specifies which features should be considered when estimating the common dispersion component; (2) "Estimate tag-wise dispersions," which allows a weighted combination of the tag-wise and common dispersion for each transcript.
The Empirical analysis of DGE tool implements the "Exact Test" for two-group comparisons developed by Robinson and Smyth (2008), which is designed to deal with situations in which many features are studied simultaneously, and where a few biological replicates are available for each of the experimental groups studied. This test also accounts for overdispersion caused by biological variability; more details can be found in CLC manuals (CLC Manuals-clcsupport.com (qiagenbioinformatics.com). Then, transformed expression values were calculated by logarithm transformation (log10), and normalized (Reads Per Kilobase Million (RPKM) normalization method: by Total = 1,000,000) (Muniz et al., 2021). Additionally, the fold change (FC) absolute value >2 and p-value (<0.01) were used, aiming to having the most significant transcripts and filtering from them the lncRNAs.

| Enrichment analysis and QLT annotation
The FEElnc classifier tool (FEElnc software) was used to classify lncRNAs as following: type of interactions: Genic, when the lncRNA gene overlaps an RNA gene from the reference annotation file, and intergenic(lincRNA) otherwise; then it is classified into subtypes and locations, which are defined according to orientation of the interactions and localization of the transcription of proximal RNA transcripts, more details can be found in the workflow ( Figure S1) and in the FEElnc github database (https://github.com/tderr ien/ FEELnc). This classification was used to explore the DE lncRNA results separately by category. The GO enrichment analysis tool in Gene Ontology (http://geneo ntolo gy.org/; release 2020-10-09) was performed to identify significant enriched biological annotations and pathways associated with genes related to DE lncRNA, through PANTHER Overrepresentation Test with a cut-off of p-value < 0.05.

Genomic Annotation in Livestock for positional candidate Loci
(GALLO) is an R package developed for the accurate annotation of genes and quantitative trait loci (QTLs) located in regions identified in common genomic analyses performed in livestock (Fonseca, Suárez-Vega, et al., 2020). Thus, we used the "find_genes_qtls_ around_markers" function in the GALLO package implemented in R (https://CRAN.R-proje ct.org/packa ge=GALLO; Fonseca, Suárez-Vega, et al., 2020) to annotate QTL overlapping with DE lncRNAs coding regions according to the bovine QTL annotation database. For that, the position of the DE lncRNAs and the bovine QTL annotation database (https://www.anima lgeno me.org/ cgi-bin/QTLdb/ index) were used as inputs. Furthermore, the "find genes qtls around markers" tool is used to annotate of genes and QTLs around candidate regions, its output is a data frame composed of the columns present in the input file and the genes or QTLs mapped within or around (if interval provided) the candidate regions. In our study, we annotated QTLs only within DE lncRNAs coding regions.

| Differentially expressed long noncoding RNA (lncRNA) annotated in the genome reference (ARS. UCD 1.2)
For tenderness and marbling traits, three DE long noncoding RNA were identified ( Table 1)

| Novel differentially expressed long genic noncoding RNA (lncRNA)
For tenderness, seven long genic noncoding RNA were identified as DE in tough beef animals in relation to tender beef animals (Table 2). Three upregulated lncRNA (lncRNA_595.1; lncRNA_71.3; and lncRNA_3440.2), and four downregulated (lncRNA_3097.2; lncRNA_775.7; lncRNA_688.5 and lncRNA_3129.5), in tough beef animals in relation to tender beef animals, were identified. These lncRNAs were located in antisense direction of intronic and exonic regions of potential coding protein transcripts. The three upregulated DE lncRNAs were interacting with mRNAs associated with glutamate receptor subtypes, monoamine plasma membrane transporters, nonvoltage-gated sodium channels, cardiomyocyte proliferation, as well as ATP-binding cassette (ABC) transporters, respectively.
For marbling, eight genic lncRNA were differentially expressed in high marbling group in relation low marbling group ( Table 2). These lncRNAs were located in antisense direction of intronic and exonic

| Novel differentially expressed long intergenic noncoding RNA
For tenderness, 19 novel long intergenic noncoding RNA (lincRNA) were differentially expressed in tough beef animals in relation to tender beef animals ( Table 3), out of them, nine were upregulated and 10 downregulated, respectively. Most of these lincRNAs were located in antisense direction (81.25%) in relation to their gene interactions. The upregulated novel lincRNA (Table 3) identified as DE for tenderness trait in tough beef animals in relation to tender were associated with the coding proteins related to RNA binding protein, metalloprotease, ubiquitin-protein ligase, ribosomal protein, and phospholipase. In addition, downregulated novel lincRNAs (Table 3) were associated with the genes related to protein classes, such as transcription co-factor, dehydrogenase homeodomain transcription factor, DNA-directed DNA polymerase, and oxidase.
Regarding the marbling trait, 37 novels differentially expressed lincRNAs were found (

| Enrichment analysis and QTL annotation
Functional classification of GO terms (i.e., molecular function, biological process, and cellular component) was performed for genes associated with DE lncRNAs for tenderness ( Figure S2) and marbling ( Figure S3) traits.
For tenderness, five significant pathways (p-value < 0.05) were identified (Figure 1  e Interaction location = are defined according to the type of interactions (genic or intergenic), assuming that intergenic lncRNA are not overlapping any RNA transcript (mRNA/gene), then it can be further classified in: upstream (the lncRNA is upstream transcribed in head-to-head or tail-to-tail orientation with RNA partner or yet both same orientation) or downstream (the lncRNA is downstream transcribed in head-to-head or tail-to-tail orientation with RNA partner or yet both in same orientation). e Interaction location = are defined according to the type of interactions (genic or intergenic), assuming that intergenic lncRNA are not overlapping any RNA transcript (mRNA/gene), then it can be further classified in: upstream (the lncRNA is upstream transcribed in head-to-head or tail-to-tail orientation with RNA partner or yet both same orientation) or downstream (the lncRNA is downstream transcribed in head-to-head or tail-to-tail orientation with RNA partner or yet both in the same orientation).
functionality of them. However, due to their significantly different expression levels ( Table 1) in animals divergent for tenderness and marbling traits, we could hypothesize that these lncRNAs are important regulators of coding-protein genes, contributing to the expression of these phenotypes.
The vast majority of the genome produces numerous lncRNAs, which comprise various RNA species longer than 200 nucleotides (nt) that are not translated into proteins (Yao et al., 2019). Studies indicate that lncRNA structure is one of the most critical factors determining its interactions and thus its specific function (Duran et al., 2019;Gupta et al., 2019). Usually, lncRNAs are categorized according to their genomic locations in relation to the nearest protein-coding.
Here, we categorized lncRNAs according to their genic location and identified two general categories: the genic lncRNA, for tenderness and for marbling traits ( Table 2), and intergenic lncRNAs, for tenderness ( Table 3) and for marbling ( Table 4) traits.
Genic lncRNAs are long noncoding RNAs overlapping codingproteins, intersecting gene exons, introns, or an entire gene respectively, and affect the functionality of their associated mRNA (Duran et al., 2019), being further classified in exonic and intronic.
For tenderness, we observed four exonic lncRNA (lncRNA_595.1, lncRNA_71.3, lncRNA_775.7 and lncRNA_688.5) overlapping part of the genes such as: PICK1, DIPK2A, RANBP3, and SLC49A3; and three intronic lncRNAs (lncRNA_3440.2, lncRNA_3097.2, and ln-cRNA_3129.5) intersecting the genes: ABCC1, GADL1, and PSMD6 (Table 2). All these lncRNAs were located in antisense direction of their associated mRNA/genes. Some of the genes associated with our DE lncRNA, such as genes PICK1, GADL1, and PMD6 were enriched in the pathways associated with Ionotropic glutamate receptor, gamma-aminobutyric acid synthesis and ubiquitin proteasome, respectively ( Figure 1). The PCK1 is a main control point for the gluconeogenesis regulation and was found involved in Inotropic glutamate receptor pathways (Figure 1), which are receptors ligand-gated ion channels permissive to cation flux across the cell membrane (Kim et al., 2020). A study comparing muscle and intramuscular fat tissue transcriptome analysis  reveals that PCK1 gene was upregulated in adipose tissue compared to muscle tissue in both buffalo and cattle. However, it seems to be more abundant in muscle tissue than in adipose tissue in cattle and may be a potential gene affecting IMF deposition, and meat quality traits. In this study, the lncRNA, lncRNA_595.1, associated with PCK1, was highly upregulated in tough meat animals in relation to tender, which may be affecting the PCK1 expression levels, contributing to the regulation of the complex transcriptional network behind the diversity of meat quality phenotype expression profiles.
The lncRNA (lncRNA_3097.2) was downregulated in tough beef animals compared to tender ( Table 2) and was associated with a GADL1 gene. Some studies have shown that GADL1 overexpression suppressed cell migration, and are associated with morphological changes in these cells (Wu, Chen, Lee, et al., 2019;Wu, Chen, Liu, et al., 2019). This gene encodes a glutamate decarboxylase like 1 protein, which is part of a large family of protein responsible for catalyzing the production of gamma-aminobutyric acid from l-glutamic acid.
This gene was enriched to one of the pathways (gamma-aminobutyric acid synthesis) associated with meat tenderness (Figure 1). Gamma-Aminobutyric acid (GABA) synthesis (l-glutamic acid +H+ 4 GABA +COz) is rapidly stimulated by a variety of stress conditions (Crawford et al., 1994) (e.g., cell stress originated by oxygen lacking) to enhance cell's resistance. Chun et al. (2014), studying the effects of γ-aminobutylic acid (GABA) on the quality and sensory properties of meat products, revealed that a GABA percentual increase contributed to a greater tenderness index, and can be used to replace NaCl of a meat product. In beef cattle, Mashima et al. (2019) correlated muscle fiber type and free amino acids (e.g., glutamine) in Japanese Black steers, expecting to find more free amino acids in slow-twitch fibers than in fasttwitch fibers due to the slow-twitch fiber mitochondria content. They reported a correlation of 0.11% and 0.63% between the proportions of myosin (MyHC1), a main muscle protein and highly expressed in slowtwitch (type 1) muscle fibers, and glutamine and gamma-aminobutyric acid, free amino acid contents in muscle, suggesting that an increase in slow-twitch fiber content induces an increase in the total free amino acid content, possibly due to high oxidative metabolism of bovine muscles. The intronic DE lncRNA (lncRNA_3129.5) was downregulated in tough beef animals compared to tender and associated with the PMD6 gene (Table 2), that enriched to ubiquitin-proteasome pathway ( Figure 1). While ubiquitin-proteasome system genes can affect muscle mass including myofibrillar protein degradation, and myogenesis inhibition (Khalil, 2018). Regardless, there is still few information about how the lncRNAs identified in this study, have been acting together with genes, and how they may be contributing to meat tenderness.
F I G U R E 1 Panther pathways associated (p-value < 0.05) with tenderness in beef cattle Therefore, future studies need to be performed to establish the function of lncRNAs in muscle tissue of tender beef animals.  (Campbell, 2005). In addition, integrin signaling pathway, that was enriched by KRAS and ITGAL genes, is involved in intracellular signals in response to the extracellular matrix including cellular shape, mobility, and mediate the progression through the cell cycle (Harburger & Calderwood, 2009 These lincRNAs were related to several genes, among them the ECE1, UBE2K, and TLE1 (Table 3). These genes were enriched in the three pathways associated with endothelin, ubiquitin-proteasome, and Wnt signaling (Figure 1). The ECE1 (Endothelin converting enzyme 1) is involved in proteolytic processing of endothelin precursors for biologically active peptides. It is a key enzyme to generate Endothelin (ET-1), a potent vasoconstrictor peptide, and has been related to hypoxia resistance in pigs (Wang et al., 2015) and human (Kon et al., 2014). Hypoxia is an important modulator of endurance exercise-induced oxidative adaptations in skeletal muscle (Kon et al., 2014). The ubiquitin-proteasome has been related to muscle atrophy; thus, a lot of ubiquitin-proteasome system genes are involved in different processes affecting muscle mass, including myofibrillar protein degradation and myogenesis inhibition (Khalil, 2018). In some gene expression studies on meat tenderness in Nellore cattle, genes related to ubiquitin metabolism have been identified (Fonseca et al., 2017;Gonçalves et al., 2018;Muniz et al., 2020Muniz et al., , 2021. Wnt signaling is one of the most important developmental signaling pathways, and controls cell fate decisions and tissue patterning during early embryonic and later development (Buechling & Boutros, 2011).
The activation of Wnt signal transduction pathways can regulate diverse processes including cell proliferation, migration, polarity, and differentiation (Eisenmann, 2005). In pigs, Zhu et al. (2016) (Figure 2). The most part of these pathways is connected to each other, which triggers various biological processes contributing to multiple molecular functions ( Figures S2 and S3) that are involved in processes such as lipid metabolism (Jesudason & Wittert, 2008), neuromuscular junctions restoration during muscle repair (Daneshvar et al., 2020), skeletal muscle regeneration (Bryan et al., 2005), skeletal maturation delaying (Balasubramanian & Crowley, 2019), skeletal muscle hypertrophy (Boppart & Mahmassani, 2019), increased glucose metabolism and reduced inflammation in skeletal muscle (Dagdeviren et al., 2016), and collagen synthesis (Purslow et al., 2012). Therefore, the regulatory function of these lincRNAs into these pathways to expression of marbling and tenderness traits should be furthermore investigated.
The OMICs sciences have contributed greatly to the knowledge of complex traits, such as tenderness and marbling (Braz et al., 2019;Magalhães et al., 2016;Santos et al., 2019), revealing genes and their mechanisms in the genomic and transcriptional levels.
However, the integration of results from these multiple approaches is still challenging, magnifying the need for efficient and accurate integrative methods to puzzle out the relationship between transcriptional regulation and several phenotypes, allowing to identify genetic variants associated with changes in the expression of economically important traits. In this study, we have performed QTL annotation into DE lncRNA regions aiming to point out QTL regions that may be selected to further meat quality investigations. As shown in Figure  genes associated with reproductive events involved in metabolism, p53 signaling, Axon guidance, and ubiquitin pathways, some of these pathways had been associated with meat quality traits in this study. Irano et al. (2016) identified 10 genomic windows located on chromosomes 5, 8, and 14 associated with the occurrence of early pregnancy and scrotal circumference. Thus, these results may be indicating a possible genetic correlation between meat quality (e.g., tenderness and marbling) and reproductive traits (e.g., scrotal circumference and age at puberty), however, studies reporting genetic correlations between meat quality and reproductive traits are scarce, which does not allow us to further speculate about those associations and how these QTLs can be affecting both, meat quality and reproductive traits.

| CON CLUS ION
This study targeted specific genomic regions of lncRNA codification and pathways that showed a multifactorial interaction that resembles transcription factor recruitment at specific genomic sites.
Thus, the findings obtained further advance our understanding of the transcriptional regulation roles of lncRNAs and the growing importance of these molecules in the muscle biology system and for meat quality traits. Therefore, further investigation to understand the interaction between lncRNA and mRNAs affecting meat quality traits, and possible pleiotropic effects regarding reproductive traits are needed.

CO N FLI C T O F I NTE R E S T
The authors have no conflicts of interest to declare.

DATA AVA I L A B I L I T Y S TAT E M E N T
The dataset utilized in this study belongs to a Qualitas Nelore breeding program company and could be available on request. The authors do not have the authorization to share the data.