Transcriptomics analysis reveals the effect of Broussonetia papyrifera L. fermented feed on meat quality traits in fattening lamb

Overview of experimental program

This paper mainly explores the molecular regulation mechanism of different ratios of B. papyrifera fermented feed on sheep meat quality. We randomly selected 12 crossbred weaned male lambs divided into four groups (three lambs per group) fed that ratio of the B. papyrifera fermented feed in the total mixed diet (0%, 6%, 18%, and 100%) for 60 days. HE staining and analysis of fatty acid content of LD muscle, a suitable ratio of fermented feed for B. papyrifera was initially selected. Then conduct transcriptomics analysis to screen the key genes that affect meat quality and predict the biological functions of these genes in fatty acid synthesis and metabolism. Finally, combining the phenotypic data and the transcriptome results, explore the effect of B. papyrifera fermented feed on the quality of lamb meat, and select the appropriate addition ratio (Fig. 1).

Chemicals and reagents

NaOH (Jinhuitaiya, Tianjin), CH₃OH (jingke, Wuxi), TRIpure Reagent (DiNing Biotech, China), CHCl₃ (Beyotime), NaCl (Jingke, Wu Xi), 14%BF3-CH3OH (Xiya, Shandong), Eosin Staining Solution (Beyotime), Hematoxylin Staining Solution (Beyotime), Evo M-MLV RT Kit with gDNA Clean for qPCR (Accurate Biotechnology (Hunan) Co., Ltd), YBR Green Premix Pro Taq HS qPCR Kit (Accurate Biotechnology (Hunan) Co., Ltd).

Ethics approval

Experimental animals were reviewed and approved by the Animal Committee of Gansu Agricultural University (GSAU-2019-76).

Preparation of experimental animals

12 crossbred weaned male lambs Dorper ♂ × Hu ♀ were selected, which with same weight (20 kg/lamb), similar ages and good health. All fed in Gansu Zhongtian Sheep Industry Co., Ltd. (Longxi, China). Randomly divided into four groups (three lambs per group) based on the feed that ratio of the B. papyrifera fermented feed in the total mixed diet (0%, 6%, 18%, and 100%). There is transition period of 7 days, a pre-test period of 7 days and a formal period of 60 days. Refer to NRC (2007) standard for 20 kg lamb with a daily weight gain of 200g/d. The diet formula was shown in Table 1. Feed the same amount twice a day at 08:00 and 18:00, respectively, and drink freely. 12 lambs were euthanized after 60 days of feeding each lamb was given intravenous injection of 360 mg sodium pentobarbital without heartbeat, continuous non-spontaneous breathing for 2–3 min, and no blinking reflex, and then dissected. Collecting the LDs: one sample was immediately placed in liquid nitrogen and subsequently cryopreserved at −080 °C for the extraction of total RNA, the second was placed in an ice box and subsequently cryopreserved at −20 °C to extract fatty acids, and the last was fixed with 4% paraformaldehyde for approximately 48 h, dehydrated in a gradient series of ethanol, cleared inxylene and embedded in paraffin. The tissues were sectioned at five mm thickness and used for HE staining.

Fatty acid analysis

Frozen samples were thawed 12 h prior to analyse at 4 °C, and weighed 5g LD muscle on an electronic balance. In order to accurately and quickly extract the lipids in the LD, total lipids were extracted with CHCl₃/CH₃OH (2:1, V/V) as described by Zlatkis, Zak & Boyle (1953) with some modifications. Fatty acid methyl esters (FAME) were prepared as following from the methods as described previously (Liu , 2015) with some modification. Extracted lipids (approximately 10 mg) were combined with BHT and 2.0 mL NaOH-MeOH solution (0.1 mol/L), saponify at 60 °C for 60 min. After the solution was clear and transparent, test the tube under 4 °C running water cool to room temperature; then add 2.0 mL of 14% BF₃-CH₃OH solution, methylate at 60 °C for 15 min, and cool the test tube to RT under 4 °C of tap water; add three mL of n-hexane, one mL of 25% NaCl solution and six mL of distilled water was shaken vigorously, centrifuged at 6,000 r/min for 10 min, the upper layer of FAME was collected, the lower layer solution was repeatedly extracted once with n-hexane, combined FAME, added an appropriate amount of anhydrous sodium sulfate, and all were sucked out with a disposable needle, and the organic phase was filtered. The membrane is placed at the front of the needle tube for filtration, concentrated to dryness under the protection of nitrogen, the residue is dissolved in n-hexane and fixed to a 2.0 mL volumetric flask, stored at −20 °C for future use, until GC detection.

Derivatized methyl esters of fatty acids were separated and quantified by gas chromatograph (7820A; Aglilent Technologies, California, USA). The capillary column was a 100 m ×0.25 mm ×0.2 µm (SPTM-2560; Sigma Louis, MO, USA). The carrier gas (Nitrogen) flow rate was 1.5 mL min⁻¹, The flow rate of gas (Hydrogen) was 40 mL min⁻¹, The flow rate of air was 400mL min⁻¹. Split injection, split ratio 60:1, injection volume 1.0 μL. Injector temperature was set as 260 °C. The initial oven temperature was programmed at 140 °C and maintained for 5 min, then increased to 200 °C at 2 °C min⁻¹, increased to 230 °C at 6  °C/min and keep 20 min.

Haematoxylin and eosin staining

Section of LD of was stained with H&E as described previously (Asuka et al., 2017). Then was observed under a microscope, and photographed using ImageView.

Muscle fiber measurement

The HE-stained section was photographed under the 20 × and 40 ×microscope using the ImageView microphase system, and the diameter and area of the muscle fibers under the 40 ×microscope were determined using ImageJ software.

Total RNA extraction, cDNA library preparation

Total RNA was extracted from each sample using TRIpure Reagent according to instructions. The concentrations and quality of RNA samples were examined by NanoDrop2000 (Thermo Fisher, Waltham, MA, USA) and Agilent 2100 (Agilent). The total RNA of the aforementioned samples was used to construct an RNA-Seq library, which were sequenced in parallel on Illumina HiSeq 2000 system by Biomarker Technologies (Beijing, China).

Identification and quantification of DEGs

Identification of differential genes used FPKM (fragments per kilobase per million reads) (Mortazavi et al., 2008). Differential expression analyses of genes among four groups were implemented using the Cuffdiff (Trapnell et al., 2013). Genes with indexes fold change (FC) > 1.5 and FDR < 0.05 were considered as DEGs.

Functional annotation and pathway analysis of DEGs

The functional of DEGs were identified through GO and KEGG database. GO annotation and KEGG pathways analyses by the GOseq R package and KOBAS software, the p value ≤ 0.05 were defined as significant enriched by DEGs.

Validation of RNA-seq results

To verify reliability of the transcriptomic profiling data, quantitative real-time PCR (qRT-PCR) was performed for nine randomly selected DEGs (FABP3, ECI1, ACADVL, GOT1, LPIN1, MLYCD, PLIN5, RASD1 and SELENOW). GAPDH was used as a reference gene in quantitative analysis. Using Evo M-MLV RT Kit with gDNA Clean for qPCR to reverse transcribe RNA into cDNA, and using SYBR Green Premix Pro Taq HS qPCR Kit and run on the Roche LightCycler96. The relative expression was calculated using the 2^–ΔΔCt method (VanSlyke & Musil, 2000). The qRT-PCR primers are listed in the electronic supplementary material, Table S1 .

Statistical analysis

All statistical analyses were conducted with SPSS 22.0 software (SPSS Inc., Chicago, IL, USA), the least significant difference (LSD with Fischer’s) method in one-way analysis of variance (ANOVA). The results were expressed as mean ±standard error. P value < 0.05 were considered as statistically significant.

Comparison of differences in fatty acid composition and content among different groups

Comparison of morphological differences of LD in different groups

HE staining of the LD of lamb in different groups, the cross-section staining results showed that the muscle fibers were polygonal, all arranged tightly, the muscle bundle gap was large, and the intramuscular membrane was obvious (Figs. 1A–1D). The longitudinal section shows that the muscle fibers are long fusiform, and the nuclei are blue elliptical or rod-shaped (Figs. 2E–2H). From Figs. 1I–1J and Table S3, the difference between the three groups of muscle fiber’s area and diameter 0%, 18% and 100% was significant (P < 0.05), while the difference between 6% and 18 was not significant (P > 0.05), but 18% was less than 6%, indicating that with the increased in the amount of B. papyrifera fermented feed the muscle fiber’s area and diameter was decreased.

Data analysis from RNA-seq

In view of the influence of different proportions of B. papyrifera fermented feed on the meat quality of fattening lamb, transcription level changes in the LD were analyzed using the Illumina HiSeq 2000 platform, and prepared to construct cDNA libraries for RNA-seq. As shown in Table 5, After filtering, the clean data of each sample reached 6.07 GB, the percentage of Q30 bases was greater than 93.57%. The percentage of clean reads that mapped into the sheep reference genome ranged from 96.79 to 97.57%. About 80% of clean reads were uniquely mapped and used for subsequent analysis.

Gene expression analysis

Comparing the three groups to yield a total of 443 DEGs, among them, 48 (19 up-regulated and 29 down-regulated), 104 (38 up-regulated and 66 down-regulated) and 198 (84 up-regulated and 141 down-regulated) DEGs belonged to 6%,18% and 100%, respectively. The 6 DEGs (2 up-regulated and 4 down-regulated) were commonly regulated by effecting of fermented feed from B. papyrifera in the three groups (Fig. 3A). A total of 89, 157 and 269 DEGs were identified in the groups of 0% vs 6%, 0% vs 18%, and 0% vs 100%, respectively (FC > 1.5, FDR < 0.05). After adding 6%, the number of down-regulated DEGs (52) were more than the number of up-regulated DEGs (37). When fed all B. papyrifera fermented feed, the number of up-regulated DEGs increased from 37 (6%) to 98 (100%), and the down-regulated DEGs increased from 52 (6%) to 171 (100%) (Fig. 3B). Hierarchical clustering represents the difference and similarity of 443 DEGs (Fig. 3C and Table S4). The results showed that there was a significant difference in gene expression profiles differences among the four groups. With the increase of the added amount, the genes with high expression at 0% were gradually down-regulated, while those with low FPKM at 0% were gradually up-regulated. It indicated that the addition of B. papyrifera fermented feed a set of gene expression or induces another set of gene expression, although these changes were to regulate the meat quality performance.

Expression pattern and functional analysis of the DEGs in LD of sheep

According to the expression profiles, 443 DEGs were classified into 8 clusters by co-expression clustering (Fig. 4 and Table S5). A total of 16 and 67 genes were classified as clusters 4 and 6, which were rapidly down-regulated within 0% to 18% and then gradually up-regulated within 100%. Most of them participated in pathways such as “Rap1 signaling pathway”, and “protein digestion and absorption”, it indicates that those genes were transiently inhibited by B. fermenta. While the genes in clusters 1 and 5 were up-regulated from 0% to 18%, and gradually down-regulated within 100%. This result indicates that those genes were transiently promoted by B. fermenta. Most of them participated in pathways such as “glucagon signaling pathway”, “PPAR signaling pathway” and “fat digestion and absorption”. The genes belonging to clusters 2, 7 and 8 were continuously down-regulated. Transcription is inhibited by B. fermented feed, most of which are rich in “MAPK signaling pathway” and “cAMP signaling pathway” pathways. Cluster 3 contains fewer genes (31), and the expression levels of these genes are significantly up-regulated with the addition of fermented feed from B. papyrifera, most of which are rich in the “AMPK signaling pathway” and “Circadian rhythm” pathways.

GO classification analyses

Annotated genes were divided into three major functional categories: biological processes (BP), cellular components (CC) and molecular functions(MF). It is obvious that there were more functional terms for BP and relatively few transcripts for CC and MF. Compared with down-regulated transcripts, detoxification (GO: 0098754), synapse (0045202), synapse part (0044456) antioxidant activity (GO:0016209) terms were peculiar in up-regulated DEGs; nucleoid (0009295) was unique in upregulation GO terms (Figs. 5A, 5B and Table S6).

KEGG annotation analyses

KEGG pathway analysis was conducted to investigate whether the genes in LD of the sheep participate in some special pathways with the addition of different B. fermenta. fermented feed. Top20 of KEGG enrichment showed that up-regulated DEGs were highlighted in “fatty acid degradation (ko00071)”, “circadian rhythm(ko04710)”, “tryptophan metabolism (ko00380)”, “arginine and proline metabolism (ko00330)”, “regulation of lipolysis in adipocytes (ko04923)”, “glycerolipid metabolism (ko00561)” and “cGMP-PKG signaling pathway (ko04022)”. Moreover, “pertussis (ko05133)”, “proteoglycans in cancer (ko05205)”, “complement and coagulation cascades (ko04610)”, “ECM-receptor interaction (ko04512)”, “pathways in cancer (ko05200)”and “Rap1 signaling pathway (ko04015)” annotated the most down-regulated DEGs (Figs. 6A, 6B, and Table S7).

RNA-Seq expression validation by qRT-PCR

As shown in Fig. 7 and Table S8, we selected 9 DEGs that were closely associated with the meat quality. The result showed that qRT-PCR expression patterns were consistent with the changing trends from RNA-seq data.

PPI network

The PPI network of DEGs related to meat quality in the LD muscle of lamb were presented in Fig. 8. There were 43 edges and 17 nodes in the network. Based on the PPI natwork, ACSL1 (degree = 9), ACADVL (degree = 9), ECI1 (degree = 9), ECHS1 (degree = 9) were the top four hub genes.

Effect of B. papyrifera fermented feed on meat quality traits

The internal quality of meat is dependent on certain characteristics of muscle fibers. Besides, there are many factors that affect the differences of muscle fibers, such as growth and development stage, gender, environment, and nutrition. The size and number of muscle fibers are important factors that affect meat quality. Therefore, the characteristics of muscle fibers are of great significance to the meat quality characteristics and growth of animals (Shen et al., 2014; Kim et al., 2018). Muscle bundles were composed of muscle fibers, and the size of muscle fiber’s diameter is also related to the shearing force. The smaller the muscle fiber’s diameter, the tenderer the meat quality. At the same time, the cross-sectional area of the muscle fiber was also one of the important indicators for judging the quality of the meat (Bidanel et al., 1991).In this study, the diameter and area of muscle fibers decreased with the increase of the ratio of B. papyrifera fermented feed (Refer to Fig. 2), indicating that the addition of B. papyrifera fermented feed was directly proportional to muscle tenderness. Therefore, this feed can be added as a protein feed in ruminant diets for improving muscle tenderness.

Lauric acid (12:0), myristic acid (14:0) and palmitic acid (16:0) in long-chain SFA are closely related to cholesterol and can harden human blood vessels, which was the reason to cause cardiovascular disease and arteriosclerosis, there are some epidemiology showing that stearic acid (C18:0) is related to coronary heart disease (Hu et al., 1999; Hunter, Zhang & Kris-Etherton, 2010). In this study, the content of lauric acid (12:0) and myristic acid (14:0) in the group100% was higher than other groups (Refer to Table 2), indicating that excessive addition of B. papyrifera fermented feed in the ruminant feed would cause long chains SFA increase. The intake of unsaturated fats can effectively reduce the level of cholesterol in the body and prevent atherosclerosis (Cameron et al., 2000; Wood et al., 2008). PUFA have been proven to prevent and treat cardiovascular diseases, and also have important physiological effects such as delaying aging and anti-cancer. MUFA plays an important role in lowering cholesterol (Harris, Poston & Haddock, 2007). The n3 series of PUFA play a variety of functions in the human body. They can reduce triglycerides in blood lipids, improve blood vessels and eliminate inflammation, and reduce the morbidity and mortality of cardiovascular disease patients (Philip & Calder, 2013). In the present study, PUFA were higher in other three groups than that in control group, and it might due to the antimicrobial activities in the diet supplemented with B. papyrifera fermented feed (Refer to Table 4), as it is reported that prenylated flavonoids isolated from B. papyrifera have antimicrobial activity (Sohn et al., 2004; Sohn, Kwon & Son, 2010), which affects the biohydrogenation of unsaturated fatty acids by ruminant rumen microorganisms, and also dietary antioxidant activity might increase PUFA content in meat (Sanchez-Muniz et al., 2012; Bellés et al., 2019 ). Combining the content of beneficial fatty acids such as DHA(C22:6n3), linoleic acid(C18:2n6c) and α-linolenic acid(C18:3n3) and harmful fatty acids such as lauric acid(C12:0), myristic acid(C14:0) and stearic acid(C18:0) in the 4 groups, it showed that the 18% group has a higher content of beneficial fatty acids and a lower count of harmful fatty acid, which would be more beneficial for the consumer’s health.

Effects of B. papyrifera fermented feed on meat quality-related DEGs and fatty acid metabolic pathways

In this study, we screened 443 genes whose mRNA abundance changes not less than 1.5 times based on FDR (Refer to Fig. 3). Intramuscular fat is one of the indicators reflecting meat quality and flavor, and is closely related to carcass and meat quality. B. papyrifera fermented feed mainly improves the intramuscular fat content of fattening lamb through nutritional regulation, thereby affecting meat quality performance. By analyzing the KEGG pathway of differentially expressed genes, it can be seen that 4 DEGs (ACSL1, AQP7, FABP3 and PLIN5) are enriched in the PPAR signaling pathway (Refer to Fig. 8). FABP3 mainly regulates the body’s fat and glucose balance in the PPAR signaling pathway, and plays an important role in the transportation and metabolism of fatty acids in cells. It can transport long-chain fatty acids from the cell membrane to fatty acid oxidation, triglyceride and phospholipid synthesis (Veerkamp & Maatman, 1995). FABP3 gene mutations and mRNA expression levels have a significant impact on intramuscular fat, which in turn affects muscle tenderness (Wang et al., 2015). ACSL1 is the most important synthetase for acyl-CoA to synthesize triglycerides, which helps fatty acid transport and triglyceride deposition (Richards et al., 2006). AQP7, a water/glycerol transporting protein, regulates adipocyte glycerol efflux and influences lipid and glucose homeostasis (Oikonomou et al., 2020). PLIN5 is expressed on both lipid droplets and mitochondria, and may participate in the interaction between lipid droplets and mitochondria (Wang et al., 2013). Muscle-specific overexpression of PLIN5 increases the storage of lipid droplets in muscle cells and also increases the rate of oxidative gene expression and metabolism (Harris-Ann et al., 2015). The formation of larger lipid droplets with a higher degree of esterification, increased methylene content and more saturated lipids (Nils et al., 2015). The above results indicate that these four genes mainly promote the deposition of intramuscular fat by participating in the synthesis of triglycerides and polyunsaturated fatty acids, thereby changing muscle tenderness. In this study, the expression of ACSL1, AQP7, FABP3 and PLIN5 genes increased with the addition of B. papyrifera fermented feed, and the expression of all genes decreased in 100% group (Refer to Fig. 7), indicating that 18% amount of added B. papyrifera fermented feed can improve meat quality.

By analyzing the KEGG pathway of differentially expressed genes, it could be seen that there were four genes (ECHS1, ACSL1, ACADVL and ECI1) that were significantly enriched in the fatty acid metabolism signaling pathway (Refer to Fig. 8). ECHS1 in cells could activate mammalian target proteins, and animal target proteins may enhance LPL expression through the PPAR pathway. LPL was involved in the metabolism and transport of lipids, and mainly hydrolyzes the chylomicrons and very low-density lipoproteins present on the surface of capillary endothelial cells to produce fatty acids (Li et al., 2014). ACADVL mainly catalyzes the oxidation of fatty acids to form C 2-C 3 double bonds, and usually catalyzes C 16-acetyl COA or even longer chain fatty acids (Aoyama et al., 1995). The long-chain fatty acyl COA is transported to the mitochondria, firstly binds to VLCAD, and interacts with the mitochondrial trifunctional protein. After 2-3 cycles of oxidation, the medium-chain fatty acyl CoA is formed (Houten & Wanders, 2010). The protein encoded by ECI1 was an important mitochondrial enzyme involved in the β oxidation of unsaturated fatty acids, which was mainly involved in the metabolic process of fatty acids (Williams & Tjian, 1991). In this study, the expression level of ECI1 increased with the addition of B. papyrifera fermented feed, reaching the highest in the 18% group, but decreased in the 100% group. The expression level of ACAVDL and ECHS1 rose to the highest after 6% B. papyrifera fermented feed, and then decreased with the increase of the added amount (Refer to Fig. 7). It showed that adding the appropriate amount of B. papyrifera fermented feed would increase the synthesis of long-chain fatty acids, while excessive addition of B. papyrifera fermented feed would increase fatty acid metabolism and reduce fatty acid content in the body, thereby affecting the fat deposition on meat quality in the LD of lamb.