Differential microRNA Expression in Longissimus Dorsi Muscle of Anhui White Goat

1. Introduction

Anhui white goat, commonly known as Huaihe goat, is a local breed of goat valued for its skin and meat. It has been included in the provincial catalogue of genetic resources for livestock and poultry. The Anhui white goat is characterised by its strong adaptability and high reproductive rate. However, it faces challenges in the form of a slower growth rate and lower meat yield compared to the Boer goat. The Boer goat is the first type of meat goat from South Africa to be introduced to our country, having come from Germany. It has a good cross effect, a fast growth rate, and a high meat yield. The report mentions that the Anhui white goat and the Boer goat have large allele differences and a rich genetic structure. Eight microsatellite loci are highly polymorphic and the genetic differentiation between the two populations is significant. The Anhui white goat has a higher genetic diversity and stronger environmental adaptability than the Boer goat. This was also confirmed in previous study, in which seven microsatellite markers were selected to analyse the genetic diversity of Anhui white goats, resulting in the detection of 45 alleles and 54 genotypes [1]. The PIC and heterozygosity of seven microsatellite markers were high. These results suggest that the Anhui white goat retains a wide range of genetic traits and is well-suited to breeding.

The quality of meat is an important economic trait of domestic animals. It is mainly determined by factors such as breed, age, sex, feeding method, feed ratio, muscle colour, pH level, water-holding capacity, tenderness, intermuscular fat content and flavour substance content [2]. In addition, from a genetic perspective, the type and composition of skeletal muscle fibers not only affect skeletal muscle growth and development, but also affect meat quality. Skeletal muscle growth and development is a complex process. In addition to mRNA, a large number of ncRNAs (non-coding RNAs) have been shown to be involved in skeletal muscle growth and development in a variety of species, including pig, sheep, chicken, and cattle [3,4,5,6].

MicroRNAs (miRNAs) are small, non-coding RNAs that are evolutionarily conserved and are usually located in intergenic or intronic regions. They repress target gene expression at the post-transcriptional level by degrading or preventing the translation of specific messenger RNAs (mRNAs). It also plays a key role in cell proliferation and maintaining metabolic stability. Many microRNAs (miRNAs) have now been confirmed to play a role in regulating myoblast proliferation and differentiation. Key specific regulators involved in myogenesis include miR-1, miR-133, miR-206 and miR-486. In addition, many microRNAs (miRNAs) that are widely expressed in skeletal muscle, such as miR-489, miR-195, miR-128, miR-29 and miR-24, are also necessary for regulating muscle development. Myoblasts are one of the most important components of embryonic muscle development and endogenous repair. Studies have shown that microRNA-27b-3p (miRNA-27b-3p) can inhibit the proliferation of skeletal muscle satellite cells in Great Black sheep [7]. Targeting ANGPT1, miR-452 can promote the proliferation of C2C12 cells and inhibit their differentiation [8]. The expression of miR-193b-3p is high in goat skeletal muscle. It activates IGF2BP1 by binding to the 3’UTR region, thereby promoting myoblast proliferation [9]. This study has confirmed that the expression of microRNAs (miRNAs) 133 and 299-5p in the skeletal muscle of Urank sheep is higher than in other tissues [10]. Previous study found that miR-487b-3p targets IRS1 through the IRS1/PI3K/Akt signalling pathway to inhibit the proliferation and differentiation of goat myoblasts. miR-378a-3p, on the other hand, can promote differentiation and inhibit proliferation by targeting HDA4 [11]. Several studies have demonstrated the important role of microRNA (miRNA) in muscle growth and development across a variety of organisms.

However, the regulatory mechanism of microRNA (miRNA) in the muscle growth and development of Anhui White goats has not been thoroughly explored. Using bioinformatics methods based on whole-transcriptome sequencing technology, we predicted and analysed the functions of microRNA (miRNA) expression profiles in the longissimus dorsi muscle of Anhui white goats and Boer goats. We also screened new molecular targets that regulate muscle growth and development in Anhui white goats. This provided a scientific theoretical basis for the breeding selection of Anhui white goats and the improvement of meat quality.

2. Materials and Methods

2.1. Sample Preparation

Three healthy Boer goats and three Anhui white goats of a similar age and body condition were selected from a sheep farm in Funan County. They were all fed the same level of nutrition. On the same day, the experimental goats were slaughtered and their longissimus dorsi muscle removed between the 12th and 13th ribs. This muscle was then placed in a cryopreserved tube and taken back to the laboratory in liquid nitrogen. A portion of muscle tissue from each goat was sent to Shanghai Yun sequence Biotechnology Co., Ltd. for whole transcriptome sequencing, while the remaining samples were stored in a refrigerator at -80 °C for subsequent use.

2.2. Total RNA Isolation, Library Construction, and Sequencing

Total RNA was extracted from the longissimus dorsi muscle of goats using the Trizol method. After extraction, the concentration of RNA in each sample was measured using a NanoDrop ND-1000 spectrophotometer, and the integrity of the RNA was assessed using denaturing agarose gel electrophoresis. The GenSeq® rRNA Removal Kit was then used to remove ribosomal RNA (rRNA) from the samples. The samples were then used to construct sequencing libraries using the GenSeq^® Low Input RNA Library Prep Kit, and small RNA libraries were sequenced on an Illumina HiSeq sequencer after quality control.

2.3. Data Quality Control and Read Mapping

Quality control and quantification of the constructed sequencing library were performed using a BioAnalyzer 2100 system. Raw reads were generated from the sequencing data following image analysis, base identification and quality control. The raw reads contain sequence and quality information for each read, as well as some reads with connectors and low quality. Using a Q30 quality control standard of >80%, the cutadapt software (v1.9.3) was used to remove sequences with adapters at the 3′ end, as well as low-quality reads with a base quality of less than 20, in order to obtain clean, high-quality reads.

2.4. Identification and Screening of Differentially Expressed miRNA

Align the trimmed reads of each sample to the combined human pre-miRNAs database (miRNA database (V22) + newly predicted pre-miRNAs). Count the number of tags aligned to each mature microRNA (miRNA) as its original expression level, and then normalise this using the transcripts per million (TPM) method. Merge the trimmed reads of all samples and use the miRDeep2 software (version 2.0.0.5) to predict new miRNAs, labelling them as ‘novel’. Calculate the fold change between the two samples and screen for differentially expressed miRNAs using the criteria of |log2FC| ≥ 1.5 and P < 0.05.

2.5. Functional Analysis of Differentially Expressed miRNAs

The target genes of differentially expressed microRNAs (DE-miRNAs) were subjected to Gene Ontology (GO) functional enrichment analysis and Kyoto Encyclopedia of Genes and Genomes (KEGG) pathway analysis. The focus was on the biological processes and signalling pathways involved in muscle growth and development that were significantly enriched by DE-miRNAs.

2.6. Construction of the ceRNA(lncRNA–miRNA–mRNA) Network

In this study, the miRNA and TargetScan software packages were used to analyse the interaction between long non-coding RNAs (lncRNAs) and microRNAs (miRNAs) and between miRNAs and 3’UTR mRNA separately. The miRNAs identified in the lncRNA-miRNA and miRNA-mRNA analyses were then intersected, and the corresponding lncRNAs and mRNAs were extracted for negative correlation analysis. A ceRNA network was constructed using Cytoscape software. 2.8 Validation of differentially expressed RNA using RT-qPCR.

2.7. Validation of DE RNA Using RT-qPCR

To verify the accuracy and reliability of the differentially expressed RNAs (DE RNAs) identified by RNA-seq, four DE mRNAs and three DE miRNAs were randomly selected for RT-qPCR. The RNA was reverse transcribed into cDNA using the PrimeScript™ RT reagent Kit with gDNA Eraser (Takara). The mRNA primers were designed using Primer Premier (version 5.0) software, and their specificity was tested at the NCBI. Due to the short miRNA sequence, reverse transcription primers were designed using the tail-adding method for qPCR primers, and qPCR reactions were performed using the TB Green® Fast qPCR Mix kit. β-Actin was used as the endogenous control for mRNA and U6 for miRNA. To verify the accuracy of the sequencing data, four differentially expressed genes (DE-mRNAs) and three differentially expressed miRNAs were randomly selected for quantitative detection by qPCR, and the relative expression was calculated using the 2^(−ΔΔ) CT method. GraphPad Prism (version 8.0) software was used to analyse the data. A T-test was used to compare the data between the two groups, and P < 0.05 was considered statistically significant.

3. Results

3.1. Total RNA Quality Test Results

Six samples of longissimus dorsi muscle tissue from Boer and Anhui White goats were collected for sequencing. Total RNA was extracted using the Trizol method, and the quality of the total RNA was assessed in terms of purity and integrity. In this study, all the OD260/OD280 values were above 1.8, indicating that the RNA was pure, the RNA bands were clear, there was no tailing phenomenon, the integrity of the RNA was good, and there was no DNA contamination. After sequencing, FastQC software was used to assess the quality of the data in terms of base and read length distributions and GC content.

3.2. Screening of Differentially Expressed mRNAs

A total of 15,954 genes were identified and screened for differentially expressed genes (DE-mRNAs) by analysing and comparing their expression profiles. A total of 1,552 differentially expressed genes were identified following statistical analysis (Figure 1A). Of these, 733 were found to be up-regulated, while 819 were found to be down-regulated. The standardised count number was used as the original data for the cluster analysis of differentially expressed genes, in order to obtain a cluster heat map. (Figure 1B). The volcano map was obtained using the difference fold and p-value (Figure 1C), with the standardised average expression serving as the original data. The log₂ transformation was performed to obtain scatter plots (Figure 1D). Further analysis of the three plots revealed significant differences between groups, strong correlations within groups, and a greater number of up-regulated genes (red) than down-regulated genes (blue). The top 10 differentially expressed genes with significantly up-regulated expression (see Table 1) and the top 10 differentially expressed genes with significantly down-regulated expression (see Table 2) are listed here.

GO functional enrichment and KEGG pathway analyses were performed on these differentially expressed genes to further speculate and annotate their possible biological functions. GO functional enrichment analysis revealed that these differentially expressed genes were enriched in 577 functional groups (P < 0.05). Of these, 555 were enriched in the biological process (BP) category (Figure 2A,B). The term with the most enriched genes was ‘ATP metabolic process’. KEGG enrichment analysis (Figure 3A,B) revealed that the differentially expressed genes were enriched in 53 pathways involved in regulating muscle growth and development. This article lists GO entries and genes (Table 3) and signalling pathways and genes (Table 4) that are related to muscle growth and development.

To more accurately explore the key genes affecting goat muscle growth and development, the GO entries in Table 2.10 were cross-referenced with the signalling pathways in Table 2.11. Twenty-eight genes were found to be present in three or more cross-references. They were MYOG, BCL2, ACTN3, NRG1, MAPK10, ADCY5, ADCY8, JUN, FOS, FGF6, CREB5, CASP3, MYLK3, GNGT1, GNG4, HGF, and PLA2G4A, SOX8, LOX, ATP2A2, DNER, CASQ2, WFIKKN1, PROX1, EGR1, TBX3, TNNI3, TNNI1. Among them, MYOG, BCL2, ACTN3, NRG1 and MAPK10 were 5 up-regulated genes, and the remaining 23 down-regulated genes. These genes were all present in the longissimus dorsi muscle of Anhui White goats and Bohr goats (P < 0.05). Four up-regulated genes (BCL2, ACTN3, NRG1, MAPK10) were significantly different (P < 0.01). There were 17 down-regulated genes (ADCY5, JUN, SOX8, FOS, FGF6, CREB5, MYLK3, HGF, PLA2G4A, ATP2A2, CASQ2, WFIKKN1, PROX1, EGR1, TBX3, TNNI3, TNNI1) were highly significant, suggesting that these 28 genes might be key candidates for muscle growth and development in goats.

3.3. Screening of Differentially Expressed miRNAs

A total of 490 microRNAs (miRNAs) were obtained in this sequencing, of which 232 were novel (Figure 4A). These 490 miRNAs were considered to be the known miRNAs of goats. Analysis of these 490 known miRNAs revealed that highly conserved, muscle-specific miRNA families such as miR-208, miR-499 and miR-206 were differentially expressed in the longissimus dorsi muscle of both the Anhui white goat and the Boer goat.

Similarly to the screening and analysis of DE-mRNAs, the heat map (Figure 4B), volcano map (Figure 4C), and scatter plot (Figure 4D) were used to visualize the DE-miRNAs. Using |log2FC| ≥ 1.5 and P < 0.05 as the criteria, differentially expressed miRNAs (DE-miRNAs) between Anhui White goat and Boer goat were screened. A total of 19 differentially expressed miRNAs (Figure 4C) were obtained, including 7 significantly up-regulated and 12 significantly down-regulated DE-miRNAs. Among the significantly down-regulated DE-miRNAs, two novel miRNAs were identified: chi-miR-novel-7_9719 and chi-miR-novel-LWLT01000027.1_26273.

In this study, target gene prediction was performed for all DE-miRNAs, and then the miRNA-mRNA network was plotted using Cytoscape software (v3.9.0). In order to explore the function of miRNAs more comprehensively, GO functional enrichment analysis and KEGG signaling pathway analysis were performed on the target genes of the top 10 DE-miRNAs.

The predicted target genes of microRNA, combined with the results of GO and KEGG enrichment, were used to screen out genes that may affect muscle growth and development (Table 5). The genes in the GO and KEGG pathways related to muscle growth and development were analysed together, and the key candidate genes that may affect goat muscle development and are targeted by DE-miRNAs were found to be BDNF, MAPT, MYB, STC1 and SMAD7. Then, a correlation analysis was performed between these key candidate genes and the expression profile of DE-miRNAs. Finally, five key candidate miRNAs that may regulate the growth and development of goat muscles were obtained. These genes were identified as chi-miR-10b-5p, chi-miR-122, chi-miR-novel-7_9719, chi-miR-143-3p and chi-miR-424-3p.

To verify the accuracy of the sequencing data, three differentially expressed miRNAs were randomly selected for quantitative detection by qPCR, and the 2-ΔΔCT theory was used to calculate their relative expression. We found that the expression trend of these three miRNAs was essentially consistent with the sequencing results. (Figure 5).

4. Discussion

In this study, the microRNA target gene prediction software programs miRanda and TargetScan were used to predict microRNA target genes. The overlapping results from these two programs were used as the final microRNA target gene predictions. GO and KEGG function analysis of the predicted target genes showed that they were mainly enriched in cell metabolism, nitrogen compound metabolism, cell biological system process, organic matter biosynthesis process, carbon metabolism, Toll-like receptor, PI3K/Akt, cAMP, MAPK and other signaling pathways. Cross analysis of genes in these GO items and KEGG pathways, combined with the network map of miRNA target genes, The key candidate mirnas (chi-miR-10b-5p, chi-miR-122, chi-miR-novel-7_9719, chi-miR-143-3p, chi-miR-424-3p) with significant differences were screened out.

Previous studies have shown that microRNA-10b-5p is widely expressed in the heart, vertical body, brain, lungs and leg muscles of chickens, and that it has evolved relatively conservatively. It may regulate the proliferation and differentiation of muscle cells via the P53 and Wnt pathways [12]. It can also encourage the growth of C2C12 myoblasts by preventing the production of NFAT5 [13] . The function of miR-122 is to downregulate the TGF-β/Smad pathway. It can inhibit the TGF-β/Smad signalling pathway by directly targeting TGFBR2, thereby restoring myogenesis [14]. Upregulation of microRNA-143 (miRNA-143) can increase the number of skeletal muscle satellite cells, but inhibit the expression of insulin-like growth factor binding protein 5 (IGFBP5), thereby inhibiting the differentiation of bovine skeletal muscle satellite cells and reducing the number of myotubes [15]. Over-expression of miR-424-5p can target HSP90AA1 and promote the proliferation of myoblasts while inhibiting their differentiation [16].

Although the regulatory effects of these microRNAs on muscle growth and development in other animals have been preliminarily explored, their potential to regulate muscle growth and development in the Anhui white goat has not yet been reported. These miRNAs are significantly differentially expressed in the muscles of Anhui White goats and Boer goats. Among them, chi-miR-10b-5p, chi-miR-novel-7_9719, and chi-miR-143-3p were significantly different, and their target genes were all found in signaling pathways related to muscle growth and development. Therefore, it is speculated that these microRNAs (miRNAs) may play a role in regulating muscle growth and development in Anhui white goats.

5. Conclusions

In this study, the whole transcriptome sequencing technology was used to perform differential expression analysis, functional enrichment analysis and ceRNA network analysis of the longissimus dorsi muscle of Anhui white goat and Boer goat. A total of 490 miRNAs were identified in this sequencing, and 19 miRNAs were differentially expressed (|log2FC|≥1.5, P < 0.05), including 7 up-regulated mirnas and 12 down-regulated mirnas. Muscle-derived mirnas such as miR-206, miR-208 and miR-499 were also found. The target genes of the differentially expressed mirnas were predicted, and the results showed that the differentially expressed mirnas targeted BDNF, MAPT, MYB, STC1, SMAD7 and other genes. GO and KEGG enrichment analysis of these target genes was performed. It was found that these target genes were mainly involved in MAPK, GABAergic synapse, HIF-1, Toll-like receptor, PI3K/Akt, cAMP signal pathway, cell metabolism, organic matter biosynthesis and other biological processes Five key candidate mirnas that might be involved in muscle growth and development were chi-miR-10b-5p, chi-miR-122, chi-miR-novel-7_9719, chi-miR-143-3p and chi-miR-424-3p.