Circulating skeletal muscle related microRNAs profile in Piedmontese cattle during different age

Piedmontese cattle is known for double-muscle phenotype. MicroRNAs (miRNAs) play important role as regulators in skeletal muscle physiological processes, and we hypothesize that plasma miRNAs expression profiles could be affected by skeletal muscle growth status related to age. Plasma samples of cattle were collected during four different ages from first week of life until the time of commercial end of the fattening period before slaughter. Small-RNA sequencing data analysis revealed the presence of 40% of muscle-related miRNAs among the top 25 highly expressed miRNAs and, 19 miRNAs showed differential expression too. Using qRT-PCR, we validated in a larger bovine population, miRNAs involved in skeletal muscle physiology pathways. Comparing new-born with the other age groups, miR-10b, miR-126-5p, miR-143 and miR-146b were significantly up-regulated, whereas miR-21-5p, miR-221, miR-223 and miR-30b-5p were significantly down-regulated. High expression levels of miR-23a in all the groups were found. Myostatin, a negative regulator of skeletal muscle hypertrophy, was predicted as the target gene for miR-23a and miR-126-5p and we demonstrated their direct binding. Correlation analysis revealed association between miRNAs expression profiles and animals’ weights along the age. Circulating miRNAs could be promising for future studies on their biomarker potentialities to beef cattle selection.


Experimental design
A pilot study to determine the presence and, the differential expression of bovine plasma miRNAs was strategized as depicted in Fig. 1. In brief, the target animals were divided into four age groups: new-born (NB), 4-6 months old (4-6M), 10-12 months old (10-12M) and 15-17 months old (15-17M) corresponding to different phases of skeletal muscle growth. Plasma miRNAs of three animals for each age group were sequenced to identify a panel of detectable and/or differentially expressed (DE-miRNAs) ci-miRNAs at four time points (Fig. 1a-c). The expression of candidate miRNAs using quantitative real-time PCR (qRT-PCR) was conducted using a larger bovine population, N = 15, N = 18, N = 18, and N = 16 animals of each age group, respectively (Fig. 1d,e). Body weight in kilograms during the four age periods beginning from birth until getting slaughtered was collected.
Results miRNAs sequencing data. Plasma miRNA were sequenced to explore the panel of expression of ci-miR-NAs in Piedmontese cattle across the four ages of growth and to identify DE-miRNA among the different groups. In total, 237 Bos taurus annotated miRNAs (bta-miRNAs) were detected. The number of reads per million (RPM) obtained for each sample after removing reads with low quality (reads without adapter, short reads and reads with multiple undetermined base calls are provided in Supplementary Table S1). Whereas Supplementary  Fig. S1 displays the size distribution plot and read length distributions that showed a peak at 20-23 nucleotides, which is the main feature of mature miRNAs.
According to the GLM analysis of sequencing data, 19 ci-miRNAs resulted differentially expressed (DE) ( www.nature.com/scientificreports/ to predict which specific molecular pathways showed possible participation in our biological setting. KEGG analysis showed statistically significant 65 pathways (Fig. 2). The significantly enriched pathways were mainly involved in energy metabolism including phosphatidylinositol 3-kinase (PI3K)/protein kinase B (AKT) signaling, signaling pathways regulating pluripotency of stem cells and skeletal muscle growth, including FoxO, transforming growth factor (TGF)-β and mTOR pathways.
miR-23a and miR-126-5p directly target the 3′-UTR of bovine MSTN. Among the 7 ci-miRNAs selected on the base of expression characteristics described above, TargetScan version 7.2 predicted MSTN as one of the target genes for miR-23a and miR-126-5p. To investigate the mechanism of miR-23a and miR-126-5p, transfection of miRNAs mimic with luciferase reporter gene cloned to the wild-type 3′-UTR of MSTN and mutated MSTN 3′-UTR were performed in both 293T cells and myoblasts (BoSC in GM). miR-27b mimic and scramble were used as positive 15 and negative controls, respectively (Fig. 5a,b). Both the miRNAs confirming the direct binding with the wild-type MSTN 3′-UTR. miR-23a significantly inhibited the dual luciferase activity around 34% and 57% in 293T and myoblasts respectively, and no inhibition of MSTN 3′-UTR-mutated was identified ( Fig. 5a). Although miR-126-5p significantly inhibited in 293T cells and in myoblasts, the binding activity was lower in both the cell systems, 11% and 43% respectively, probably ascribable to the different context score attributed to TargetScan (Fig. 5b).

miR-23a and miR-143 overexpression modulates skeletal muscle hypertrophy-related target genes. Several published results suggested that MSTN as well as Insulin Like Growth Factors (IGFs)
were crucial mediators in physiological skeletal muscle processes, including development, hypertrophy, and regeneration [33][34][35] . Based on the above results related to luciferase activity inhibition, we decided to proceed To test the ability of miR-23a and miR-143 to alter the expression of their targets in myoblasts under GM and DM conditions, the overexpression of both miRNA mimics was induced and the post-transfection expression of miR-23a and miR-143 and their target genes MSTN and IGF1R respectively, at mRNA levels were defined through qRT-PCR (Fig. 5c,d).
In proliferating myoblasts (GM), comparable results were obtained with the overexpression and the inhibition of miR-143 with significant different levels of IGF1R expression when compared with control after 24 h of transfection (Fig. 5d). Overexpression of miR-143 (log2fold change = 32, p < 0.05, Fig. 5d) significantly reduced the expression of IGF1R (Fig. 5d). In differentiated condition, no significant differences in IGF1R expression emerged as well as the inhibition of miR-143 showed no different significant expression in both the myoblasts culture conditions (Fig. 5d).

Relationships between miRNAs expression and animals' body weight.
To further deepen the relationship between the animal weight and the expression of all the 14 ci-miRNAs validated in qRT-PCR, we analysed the Spearman's rank correlation coefficients. Considering all the time points, the statistical test high- www.nature.com/scientificreports/ lighted that miR-126-5p and miR-146b were positively correlated to animals' weight, rho = 0.26 (p = 0.02) and rho = 0.36 (p = 0.002) respectively ( Fig. 6a,b), while a negative correlation was observed for miR-223 expression and beef cattle weights, rho = − 0.39 (p = 0.001) (Fig. 6c).

Discussion
In the recent years, tissue and circulating miRNAs have gained a lot of attention due to their legit fine-tune role in orchestrating several physiological processes, metabolic changes and adaptive response pathways underlying skeletal muscle growth, differentiation, and hypertrophy 37 . However, mostly the studies have been carried out on tissue miRNAs in humans and rodent models. The knowledge about ci-miRNAs in livestock and particularly in bovine is limited. To characterize the changes in skeletal-muscle related ci-miRNAs during the growth of Piedmontese beef cattle, we examined the plasma samples of animals categorized into four groups based on the age and corresponding weight. The goal of this study was to characterize the plasma miRNA signatures of cattle from first week of life until the time of commercial end of the fattening period before slaughter to integrate information between ci-miRNAs expression, muscle physiological age-growth, and animals' weight. We adopted a multistep approach, using firstly a miRNAomics approach to identify DE-miRNAs in plasma, followed by validation of miRNAs profiling results through qRT-PCR. Furthermore, the investigation of ci-miRNAs along with gene targets involved in skeletal-muscle physiology mechanisms including hypertrophy was carried out demonstrating MSTN as new direct target of miR-23a (Fig. 4a) and miR-126-5p (Fig. 4b). At last, we found correlation between ci-miRNAs expression and animals body weight.
In our experiments, the levels of myomiRNAs including miR-1, miR-133a/b, and miR-206 were found negligible in all the age groups on sequencing. Muroya and colleagues reported a similar result obtained by microarray analysis done in the plasma of 22 months aged Japanese beef cattle 38 . Circulating levels of myomiRNAs are described increased in patients of muscular pathologies 39,40 , and their higher abundance in blood stream of patients with skeletal muscle pathologies suggest a possible relation with the injury of the fibres.
Among the DE-miRNAs identified by qRT-PCR, miR-10b, miR-126-5p, miR-143 (Fig. 3a), and miR-146b were found to be increased in expression along the age-groups when compared with NB and reached the highest expression at 10-12 months for the first three miRNAs, and at 15-17 months of life for miR-146b (Fig. 4a). www.nature.com/scientificreports/ Previously, miR-10b, miR-143 and mir-146b were described as involved in PI3K/AKT pathway in the skeletal muscle tissue and myoblast cells 18,41,42 . It seems that all these miRNAs negatively regulated the myoblasts proliferation. A study reported that miR-10b significantly suppressed PIK3CA expression and decreased PI3K/Akt/ mTOR pathway activity, promoting expression of TGF-β 43 . Furthermore, miR-10b expression steadily decreased during myoblasts proliferation, but significantly increased during myoblasts differentiation in C2C12 mouse cell line by binding with NFAT5 gene (isoform of Nuclear Factor of Activated T cells) which is recognized as a key regulator in myoblast migration and differentiation 44 .
Additionally, miR-10b was previously identified upregulated in grazing cattle plasma compared with grainfed animals coincidently with the expression change of miR-10b in the Longissimus dorsi. This observation was correlated with possible alterations of skeletal muscle regulation through PTEN targeting (the phosphatase that converts phosphatidylinositol 3,4,5-triphosphat), involved in muscle cell differentiation and hypertrophy 45 .
Likewise, miR-143 was demonstrated to have a role in skeletal muscle and adipose tissues through the regulation of target genes involved in PI3K, IGF and FoxO signaling pathways such as IGF-I, IGF-II, IGF-IR, and IGFBPs 18,46 . The pathway analysis prediction (KEGG) of our sequencing data revealed a relation between the same pathways and ci-miRNAs identified in bovine plasma (Fig. 2).
In in vitro bovine cell models, miR-143 exerted important function in bovine myogenesis by targeting IGFBP5 (insulin like growth factor binding protein 5) 18 , which is an important component of IGF signalling pathway, but also its expression increased in differentiating bovine satellite cells 47 . At the end of pigs' fattening period, miR-143 expression was found upregulated in the psoas muscle, assuming its involvement in processes for the transformation of muscle fibre type 19 , but an important positive role of miR-143 was already identified in adipogenesis mechanisms 48 and recently demonstrated highly expressed in the adipose tissue of cows 49 .
In our findings, the increased levels of expression of circulating miR-10b and miR-143 during 10-12 months of age could be suggestive of their involvement in differentiation processes when the animal substantially gains the body weight, but the increase in miR-143 expression could be also linked with adipocyte differentiation. So, it can be hypothesized that this trend might mediate the maintenance of proliferating status of myoblasts In addition, we found miR-146b differentially expressed in the plasma when NB subjects and the elder groups were compared (Fig. 4a) and it was positively correlated with the body weight gain during the four growth points (Fig. 6b). In mouse and chicken models, miR-146b was described as a positive regulator of myogenic differentiation, acting through multiple targets (e.g., Smad4, Notch1, and AKT1) 20 but this function was demonstrated in tissue. At present, no information is available about the role of plasma miR-143 and miR-146b in cattle, but their differential expression among ages in the bloodstream might attribute to an important role in the regulation of growth mechanisms probably related also to skeletal muscle 42 .
MSTN is a secreted growth and differentiation factor that belongs to TGF-β superfamily. MSTN gene mutation or targeted degradation of its mRNA has shown increased muscle mass attributing to muscle hypertrophy phenotype 15 .
In our qRT-PCR assay, miR-23a showed high expression in terms of Cq values in all the age-groups. Even if not differentially expressed, this result combined with the targeting of MSTN mRNA by miR-23a are suggestive of an important role of this miRNA in Piedmontese breed affecting the double-muscled phenotype. Our previous study has shown functional effect of miR-27b on MSTN and as one of possible contributing factors of hypertrophy in Piedmontese cattle 15 . Here we have demonstrated for the first time the direct binding between miR-23a and MSTN gene. Indeed, both in cardiac 50 and in skeletal muscle 51 it was previously demonstrated that miR-23a can mediate the hypertrophic signals, through the binding of an anti-hypertrophic protein (NFATc3) and two muscle-specific ubiquitin ligases (MAFbx/atrogin 1 and MuRF1) that promote hypertrophy by protecting muscle atrophy-associated protein degradation. These findings suggest that the high level of miR-23a expression in Piedmontese breed may have role in modulation of skeletal muscle mass. Future studies are required to specifically elucidate the role of miR-23a and if it can exert a synergic effect with other miRNAs such as miR-126-5p and miR-27b in the hypertrophic program.
Interestingly, miR-223 showed a decreased expression in the elder animals when compared with NB and its expression was negatively correlated with their weights. TargetScan prediction indicated IGF1R as a putative target of miR-223. Recent studies have described miR-223 as involved in various myocardial disorders related with cardiomyocytes hypertrophy also through the IGF1R binding 52,53 . It seems that the NB subjects with higher levels of miR-223 show higher body weight at 10-12M, further analysis will confirm with a larger sample number if miR-223 could be defined as predictive marker in differential growth.

Conclusion
For the first time, our findings provide evidence of high numbers of skeletal-muscle related miRNAs in beef cattle plasma of different age-groups. Among these ci-miRNAs, miR23a and miR-126-5p were demonstrated to directly bind MSTN mRNA in a bovine satellite cells in vitro model. Three of DE ci-miRNAs (miR-146b, miR-126-5p and miR-223) have shown correlation between miRNAs expression value and body weight during different age. All these results are the basis for future investigations about a possible role of ci-miRNAs as biomarkers of muscle post-natal growth and development. Remarkably, the information about plasma miRNAs associated with production traits in beef cattle could prove to be promising biomarkers for the genetic selection of meat-purpose www.nature.com/scientificreports/ animals for future breeding. The pattern of miRNAs expression and its correlation with age growth status are, obviously, intriguing to predict at calf birth and the trend of weight during the first months of life.

Material and methods
Ethics statement. Blood  Blood collection and plasma separation. Blood samples were collected in 10 ml K2 EDTA Vacutainer tubes (Becton Dickinson, USA) by jugular vein puncture, using 18G needles (Becton Dickinson) and instantly stored in ice. Within 2 h of collection, samples were centrifuged at 1900g for 10 min at 4 °C to remove blood cells, followed by second centrifugation at 16,000g for 10 min at 4 °C to remove cellular debris and platelets 54  For the RNA extraction from cells, TRI Reagent (Sigma-Aldrich, St. Louis, MO, USA) was used following the manufacturer's protocol. Cells were homogenized in 1 ml of TRI reagent and RNA pellets were resuspended in variable amount of RNAse free water corresponding to the size of pellet. RNA quantity and quality were determined using Nanodrop ND1000 Spectrophotometer (Thermo Fisher Scientific, USA). The ratio of the optical densities measured at 260 and 280 nm was > 1.9 for all the RNA samples. cDNA was synthesized using iScript Reverse Transcription Supermix for qRT-PCR kit (Bio-Rad, USA) for each sample taking the volume of eluted miRNAs sample equivalent to 500 ng. . Libraries were then prepared for sequencing and sequenced on single-end 75 bp mode on NextSeq500 (Illumina, San Diego, CA). The Bcl2Fastq 2.0.2 version of the Illumina pipeline was used to process raw data for both format conversion and de-multiplexing. Adapter sequences were masked with Cutadapt v1.11 from raw Fastq data using the following parameters: -anywhere (on both adapter sequences) -overlap 5 -times 2 -minimum-length 35 -mask-adapter. Lower quality bases and adapters were trimmed by ERNE v1.4.6 software. Above mentioned small RNA-Seq procedures were commercially performed by IGA Technology Services (IGAtech), Udine, Italy (Fig. 1b). GLM analysis was performed cumulatively for all the four groups using DESeq software taking into consideration log2 fold change > ± 0.58 and p < 0.05 55 . miRNAs pathway analysis and target gene prediction. Bovine miRNAs were related to their miRNA family through TargetScan 7.2 association table 56 . miRNA-target genes were associated to our miRNAs of interest by selecting from TargetScan 7.2 predicted target genes with at least one conserved binding site and a cumulative weighted context++ scores (CWCS) > − 0.05 57 KEGG (Kyoto encyclopaedia of genes and genomes) enrichments and functional analysis were conducted through Bioconductor 58 package ClusterProfiler, version 3.12.0. A 0.05 cut off value was chosen for both p value and q value and BH (Benjamini-Hochberg) and False Discovery Rate control method for multiple testing was considered. Human functional annotations were based on org.Hs.eg.db packages. All analyses were run in R, a free software environment for statistical computing and graphics, release 3.6.3 59,60 .
Quantitative assessment of miRNAs and mRNA target gene expression. qRT-PCR with total cDNA was performed using SYBR green II PCR Kit (Qiagen, Germany). PCR amplifications were performed on Bio-Rad CFX Connect Real-Time System (Bio-Rad, Hercules, CA, USA). The qRT-PCR parameters specific for miRNAs and gene expression are detailed in Supplementary Table S4. Differential expression among miRNAs was carried out by comparing the normalized Cq values (ΔCq) for all biological replicates between the two group of samples. Based on stable number of sequencing reads and stable value of qRT-PCR expression among the groups, miR-378 was selected as the reference gene. The fold change of expression of transcript/miRNA was calculated by 2 −ΔΔCq method where ΔCq = Cq of the target gene/miRNA − Cq of the reference gene/miRNA 12  www.nature.com/scientificreports/ For the quantitative expression analysis of target genes, cDNA (5 μg) was prepared in a single run to perform qRT-PCR experiments for all the selected genes. To determine the relative amount of specific IGF1R and MSTN transcripts (Supplementary Table S3), the primers for target and reference genes were designed on Bos taurus GenBank mRNA sequences using Primer 3 Software (version 4.0). Oligonucleotides were designed to cross the exon/exon boundaries to minimize the amplification of contaminant genomic DNA and were analysed with the IDT tool 61 for hairpin structure and dimers formation. Primer specificity was verified with BLAST analysis against the genomic NCBI database. To establish primers efficiency, the dilution method was used. Hypoxanthine phosphoribosyl transferase 1 (HPRT-1) gene was used as a reference gene for RNA concentration and reverse transcription efficiency 12 . Tables 1 and 2 summarize primer assays information including sequences and target gene assays with gene accession number and amplicon sizes. Primers were pre-designed and commercially synthesized by LNA technology (Qiagen, USA) (see Supplementary Table S2).
Transfections and dual-luciferase reporter assay. 293T cell line and BoSC were used to perform transfections for dual-luciferase assay. Cells were seeded at the count of 100,000 cells per well of a 24-well plate. At 60-70% confluence cells were transfected. A total of 300 ng of the appropriate psicheck-2 luciferase reporter construct cloned with 3′-UTR of the MSTN transcript 15 was co-transfected with 30 nM final concentration of double-stranded RNA oligonucleotides designed to mimic miR-23a, miR-126-5p (Qiagen, Germany) and miR-27b (Exiqon, USA) molecules using Lipofectamine 2000 (Invitrogen, USA) according to the manufacturer's protocol. Twenty-four hours after transfection, the medium was changed, and cells were grown for additional 24 h before assay. Firefly and Renilla luminescent signals were quantified according to the manufacturer's instructions by Dual-Luciferase Reporter Assay System (Promega, USA) with a VICTOR Multilabel Counter luminometer (PerkinElmer, Waltham, MA). All values are given relative to transfections with the appropriate negative (Scramble) and positive (miR-27b) controls.
To study the effect of selected miRNAs on their target gene expression, BoSC were cultured in 6-well plates in GM and were transiently transfected with 30 nM of specific mimics and 100 nM of inhibitors of miR-23a and miR-143 (Qiagen, Germany) along with a negative control scramble (Exiqon, USA) using Lipofectamine 2000 (Invitrogen, USA). After 24 h, BoSC were harvested for the expression analysis of MSTN and IGF1R by qRT-PCR. Statistical analysis. DESeq software package was employed to perform the GLM analysis and pair-wise analysis on sequencing data. Statistical significance concerning qRT-PCR data was assessed on GraphPad Prism (9.0.2 version) by one-way ANOVA test and Student's t test. Distribution normality assumption was assessed through Shapiro Wilk test. Data were expressed as mean ± SEM. Differences were considered as significant at a level of p < 0.05. Spearman's rank correlation coefficient was evaluated by cor.test function in R and used to test possible relationship between animal weight (expressed in Kg) and miRNA expression (https:// www.r-proje ct. org/). The level of significance has been set at p < 0.05 and trends were considered at 0.05 ≤ p < 0.1.
Ethical approval and consent to participate. The cattle blood samples used in this scientific work were collected after the authorization of Ethical Animal Welfare Committee of Department of Veterinary Science, University of Turin (Prot. No. 663). All methods were carried out in accordance with relevant guidelines and regulations. Methods are reported in the manuscript following the recommendations in the ARRIVE guidelines.

Data availability
All data generated or analysed during this study are included in this published article [and its supplementary information files].