Muscle transcriptome provides the first insight into the dynamics of gene expression with progression of age in sheep

The dynamic synergy of genes and pathways in muscles in relation to age affects the muscle characteristics. Investigating the temporal changes in gene expression will help illustrate the molecular mechanisms underlying muscle development. Here we report the gene expression changes in skeletal muscles through successive age groups in Bandur, a meat type sheep of India. RNA sequencing data was generated from the longissimus thoracis muscles from four age groups, ranging from lamb to adult. Analysis of 20 highest expressed genes common across the groups revealed muscle protein, phosphorylation, acetylation, metal binding and transport as significant functions. Maximum differentiation was observed after 2.5–3 years on transition from lambs to adult. Transcriptional regulation by the TFAP2 transcription factors, IL-6 signaling and PI3K/AKT signaling pathways were enriched in younger animals. The gene-protein network demarcated key interactive genes involved in muscle development and proliferation that can be used as candidates for future research on improvement of muscle characteristics.

www.nature.com/scientificreports/ genes were muscle contraction, glucose metabolism, respiratory electron transport and creatine metabolism (Table S3). The expression profile of some genes previously associated with meat quality in livestock species 11 like HSPB1 and HSPB8 conspicuously increased after 3 years of age. Expression of DNAJB5 showed a mild increase whereas that of HSPA6 decreased with age. Transcription of CAPN3 also exhibited a very slight increase with advancement of age in investigated animals. Other genes like FABP3, CAPN1 and CAST did not express prominent difference across the groups. Similarly, collagen genes are also known to influence meat quality 15 . Although the expression of some collagen genes, namely, COL15A1, COL23A1, COL4A2, COL6A2, COL1A1, COL4A1, COL1A2 and COL3A1 was observed to be quite low (≤ 7.3 RPKM), their expression was highest in young animals and decreased considerably with age.
Differentially expressed genes and enriched pathways. Differential expression of genes was analyzed across all combinations of the four age groups. Considerably larger number of genes was observed to be differentially expressed between Group1-3 and Group1-4 only, which were therefore used for further analysis (Fig. 4). There were 301 significantly differentially expressed genes in Group 1 versus 3 (p adj < 0.05; FDR < 0.05), followed by 141 in Group 1 versus 4. The up-and down-regulated genes observed in Group 1 versus 3 were 222 and 79, respectively, whereas 98 up-regulated and 43 down-regulated genes were identified in Group 1 versus 4 (Table S4). Gene ontology analysis revealed negative regulation of cell proliferation, response to hypoxia, positive regulation of GTPase activity, negative regulation of cysteine-type endopeptidase activity involved in apoptotic process, calcium ion binding etc., as significant functions (p adj < 0.05) associated with differentially expressed genes in Group 1 versus 3 (Fig. 5a). Significant pathways associated with the up-regulated genes were transcriptional regulation by the TFAP2 (AP-2) family of transcription factors, IL-6 signaling, PI3K/AKT signaling, ECM proteoglycans, muscle contraction etc. (Table S5), while down-regulated genes were related with transcriptional regulation by RUNX1 and RUNX3. Similarly, enriched functions in Group 1 versus 4 included regulation of transcription from RNA polymerase II promoter, myoblast fusion, regulation of epithelial to mesenchymal transition and protein phosphatase regulator activity (Fig. 5b). The top canonical pathways identified again were transcriptional regulation by the TFAP2 (AP-2) family of transcription factors, IL-6 signaling, PI3K/ AKT signaling. Other important pathways included MAPK1 and MAPK3 activation, DAG and IP3 signaling and PKA activation (Table S6).

Discussion
Muscle traits relevant to meat quality are modulated by age, nutrition, genetics and environmental factors which in turn affect the structure and composition of muscle fibres as well as intramuscular fat. Differential expression of genes in muscles between different breeds in relation to feed, carcass and muscle properties has been examined by several studies [8][9][10]12 . It is well established that muscle characteristics are affected by age 16 . However, investigating the simultaneous change in gene expression with age in skeletal muscles of sheep has not garnered much interest.  www.nature.com/scientificreports/ The present study therefore, explored the pattern of gene expression with progression of age in skeletal muscles of Bandur sheep. The sarcomere proteins, mainly alpha-actinins, myosin heavy chains and Z disc form the basic units of muscle fibres and are coded by ACTA1, ACTN2, ACTN3, MYH4, MYH7 and MYOZ genes 17,18 . The troponin complex modulates the interaction between actin and myosin during muscle contraction 19 . TNNI2 and TNNC2 encode for their respective troponin isoforms and are known to play a key role in muscle composition and aging 20 . The genes mainly involved in muscle contraction and muscle fibre composition (ACTA1, ACTN2,  ACTN3, ANKRD1, TCAP, NRAP, MYH4, MYH8, TNNI2, TNNC2, MYOZ1, LDB3, TPM1) were abundantly expressed irrespective of age group in our study. Earlier investigations have revealed no significant effect of sex on muscle fibre type, however, the fibre diameter was affected 21 . Hence, further research on gene expression in longissimus thoracis muscles, across both rams and ewes is warranted. Besides structural components of muscles the highly expressed genes also represented glycogen metabolism and storage (LDHA, PYGM, PGAM2) as well as energy metabolism (COX1, COX3, CKM, ND4). The calpains are another family of genes associated with muscle tenderization 22 . The calpain-calpastatin activity was reported to increase with age in beef cattle 23 . It is well established that the activity of calpains is enhanced during aging 24 . A slight increase in the expression of CAPN3 with age was also detected in Bandur sheep. The domination of expression of these genes in our study reinforces their fundamental role in muscle physiology. Major pathways identified by differential analyses also had relevance to myogenesis. During muscle development, the transition from epithelial to mesenchymal cells is known to be regulated by TFAP2 25 . Studies on mice have demonstrated that cytokines of the IL-6 family play a crucial role in regulating myogenesis 26 . It is interesting to note that this activity is in turn dependent on MAPK and NFκβ signaling 27 . Recent studies on Bandur sheep have also reported the enrichment of the PI3K-Akt and MAPK pathways that are associated with hypertrophy as well muscle differentiation 12,13 . These pathways were enriched in the young animals.
As expected, lesser differentiation of gene expression was observed between successive age groups. Differential expression was prominent after 2.5-3 years in Bandur sheep reflecting transition from lamb to adult sheep. Major differences observed between lambs and adult sheep muscles were evident in the expression of genes associated with muscle structure, growth and lipid metabolism. The gene-protein network demarcated important interactive hub genes that may be involved in age related muscle development in sheep. The hub genes enriched in younger animals included COL1A2, COL1A1, TWIST1, NR4A1, VIPR1, CCNE1, ESR1, IL17A and ISG15.It is worth mentioning that these highly connected genes identified in our study are known to be involved in muscle growth and regeneration. Collagens form the major connective tissue that imparts structural strength to muscles 28 . The crosslinking between the collagen molecules increases with age resulting in increased toughness of meat in adult animals 18 . The expression of COL1A1 and COL1A2 was observed to be decreased in bovine fetuses that were myostatin deficient 29 . A decrease of collagen expression in older animals was also observed in mice 30 . Several studies have associated the significance of expression of collagen type 1 gene in muscle composition and regulation of intra muscular fat in pork 31,32 . Age affects the structure and distribution of collagen in the extracellular matrix of skeletal muscle 15 , which is also reflected in the change in expression of collagen precursor genes in skeletal muscles of Bandur sheep. TWIST1 is another gene encoding a transcription factor that is involved in the early growth of skeletal muscles but its function in adults is less explored 33 . It is also known to inhibit DNA binding of MYOD1 during myotube formation 34 . Other genes with higher expression in lambs that are associated with muscle growth include NR4A1 and VIPR1. Experimental evidence using NR4A1 knockout mice suggests its association with growth and proliferation of muscle cells 35 . Studies on mouse myoblasts revealed NR4A1 as an important factor regulating myoblast differentiation 36 . VIPR1 is a G-protein coupled receptor expressed in the skeletal muscles that may be involved in muscle atrophy/hypertrophy 37 . Further, experimental evidence supports the role of the estrogen receptor (ESR1) in muscle strength by modulation of myosin regulatory light chain in mice 38 . The ESR1 gene regulates expression via the estrogen signaling pathway. The IL17A gene, on the other hand has been reported to activate PPARγ which is involved in increased adipogenesis 39 . In addition, NR4A receptor family in conjunction with CRTC2 and PPARγ assists in myogenesis and adipogenesis via the CREB pathway 40 . Thus, most of the highly connected genes identified in Bandur lambs were associated with growth, proliferation and adipogenesis.
Notable hub genes in older or mature animals included MYOD1, RUNX1, TLE3. SOCS3, CBLB, HSPB1, IFIT1, CDLK2, NR4A2 and PER1. Myoblast differentiation is orchestrated by an assortment of muscle-specific regulatory factors. MYOD1, a myogenic regulatory factor and an important marker for myoblast development, facilitates myogenesis during embryogenesis. MYOD1 along with RUNX1 and TFAP1 transcription factors, has been associated with myoblast regeneration or proliferation 41 . Although the function of MYOD1 in adult muscle is less investigated, there are some reports that have established its role in growth and proliferation in adult muscles 42,43 . RUNX1 has also been suggested to contribute to myotube generation in humans 44 and its expression is activated by muscle damage 41 . TLE is a family of transcriptional co-factors that are involved in cell differentiation. The transcriptional activity of MyoD has also been reported to be regulated by TLE3 45 . Another important gene SOCS3, identified in this study has been linked to myogenic differentiation 46 . Previous reports suggest the role of CBLB gene in skeletal muscle atrophy in mice by negative regulation of growth factors during cell differentiation and development 47 .
The small heat shock proteins (Hsps) long associated with heat stress, are being increasingly implicated in cellular senescence and apoptosis 48 . The HSPB1 is highly expressed in muscle tissues and codes for Hsps 49 . The Hsps are anti-apoptotic and prevent the muscles from denaturation, which has implications in meat tenderness 50 . Several studies have reported the involvement of HSPB1 in meat tenderness in cattle 51,52 . Age related expression of genes associated with meat quality detected considerable increase in the expression of HSPB1 and HSPB8 in older age (3.5-4 years) in Bandur sheep. A similar increase in hspb1 protein was observed in the skeletal muscles of chicken 53 . Although much information on HSPB8 and aging are not available, it is known to be involved in protein damage through translational arrest and autophagy 54  www.nature.com/scientificreports/ significant effect of sex on the tenderness of longissimus thoracis muscles 55 . Therefore, it may be speculated that a similar expression of hsps could be expected across both sexes. The role of some of the identified hub genes is not very well defined in skeletal muscles. These include the interferon stimulated genes ISG15 and IFIT1 that are suspected to be involved in muscular diseases 56 . Although inadequate information is available on CCNE1 (Cyclin E1) and Cyclin Dependent Kinase (CDK) gene in muscle tissues, both participate in the cell growth during the cell cycle 57,58 . CDKL2 is known to be involved in cell proliferation in human cancer cells 59,60 , while PER1 gene of the circadian clock is known to regulate cell proliferation and apoptosis in association with human cancer cells 61 . Further investigations are required to determine the precise role of these genes in sheep muscle development.

Conclusion
This study is an attempt to gain an insight into the dynamics of gene expression with progression of age in sheep skeletal muscles. Our results highlight the temporal changes in gene expression in skeletal muscles from lambs to adult sheep. Noteworthy demarcation between lambs and adult sheep was observed in the expression of genes associated with muscle growth, adipogenesis and apoptosis. The enrichment of PI3K-Akt and MAPK pathways in lambs suggests their relevance in muscle growth. This differential expression analysis will contribute towards understanding the genetic basis of physiological changes in muscles with age. Several highly connected genes identified in our study are known to be involved in muscle development and proliferation and may serve as candidates for future research in myogenesis.

Materials and methods
Samples. Healthy animals were selected in coordination with the sheep rearers and butchers from Mandya district, Karnataka. As the animals were selected from the field, they were considered healthy if they did not show any visible signs of infection or ailment and were inspected by a veterinary officer. Care was taken to select unrelated male animals from flocks reared under the same management and feed. Since no records of the animals were maintained, the age of the animals was determined by information provided by the rearers as well as the dentine pattern. A two-tooth stage confirmed the age about 12 months (Group 1), 4-tooth between 18 and 24 months (Group 2), 6-tooth between 30 and 36 months (Group 3) and 8-tooth between 42 and 48 months (Group 4). Four animals in each group were selected. The animals were slaughtered according to standard commercial 'halal' procedures.
RNA sequencing. The longissimus thoracis muscles were collected in RNALater (Thermo Fisher Scientific Baltics, UAB, Vilnius, Lithuania) solution. RNeasy kit (Qiagen, Hilden, Germany) was used for extraction and purification of total RNA. Samples with RIN value ≥ 7.0 (Agilent Bioanalyzer) were selected for further processing. Libraries of four biological replicates from each age group were prepared by TruSeq RNA Library Prep Kit v2 (Illumina, San Diego, CA, USA). Paired end sequencing (150 bp) was performed on Illumina HiSeq-2000 Platform.
Data analysis. FastQC (v 0.11.5) was used for assessment of the quality of samples 62 . The raw reads were trimmed and filtered using FastXToolKit. CLS Genomics Workbench 6.5.1 (CLC Bio, Aarhus, Denmark) was used for data analysis. The reads were mapped against the ovine genome assembly (Oar_v4.0). The RNA sequencing data have been deposited in the NCBI Bio Project PRJNA416678 with accession numbers SAMN16191735-750. Normalization of reads was done as reads per kilobase million (RPKM) and reads with RPKM values > 0.01 were included in the study. The CLC transcriptomics analysis tool was used for differential expression analysis across the different age groups.
Validation by real time quantitative PCR (RT qPCR). Ten differentially expressed genes were selected for RT qPCR analysis (Table S7). Primers for selected genes were designed by Primer 3 software 63 or taken from published sequences [64][65][66] . cDNA was synthesized from 2 µg of purified RNA, using Super Script III Reverse Transcriptase (Thermo Fisher Scientific, Carlsbad, CA), as per manufacturer's protocol. Each sample was analyzed in triplicate qPCR reactions. The final reaction volume of 10 µl contained 2 µl of cDNA, 8 µl of qPCR master mix (5 µl of SYBR Green Real-Time master mix (Applied Biosystems,Vilnius, Lithuania), 0.3 µl of each primer, 2.4 µl of DNA/RNA-free water). The samples were run on QuantStudio 5 Real-Time PCR System (Applied Biosystems). Standard curve calculation using four points of cDNA serial dilutions was used to estimate the PCR efficiency. Statistical analysis. Differentially expressed genes with log 2 fold change ≥ 2.0 and p value (p adj ) < 0.05 were used for further analyses. DAVID 67 , Consensus Pathway Data Base 68 and Reactome were utilized for the functional and pathways analysis 69 . Cytoscapever 3.6.0 70 along with cytoHubba app was utilized for gene-protein network analysis 71 . For the RT qPCR the mean cycle threshold (Ct) values of the genes were normalized to geometric mean of the reference genes PPIB and β-ACTN 72 . The 2 −ΔΔCT method was used for data analysis 73 . Ethics approval. Animal  www.nature.com/scientificreports/