畜牧兽医学报  2022, Vol. 53 Issue (1): 53-65. DOI: 10.11843/j.issn.0366-6964.2022.01.006    PDF    
影响藏猪与大约克夏猪肌肉品质的lncRNA的筛选及功能分析
宋明坤1, 薛明明1, 张力戈1, 颜铎1, 夏宁1, 商鹏2, 李新建1, 韩雪蕾1, 乔瑞敏1, 李秀领1, 李明1, 王克君1     
1. 河南农业大学动物科技学院,郑州 450002;
2. 西藏农牧学院动物科学学院,林芝 860000
摘要:旨在分析藏猪和大约克夏猪(180日龄)背最长肌组织中与猪肉质性状相关的长链非编码RNA(lncRNAs)和基因,探讨lncRNAs在猪肉品质中的分子调控作用,为猪肉质性状的改良提供理论依据。本研究根据先前公布的数据,以饲养在相同条件下的180日龄健康藏公猪(n=3)和大约克夏公猪(n=3)为研究对象,收集背最长肌肉品质表型数据,利用相同个体背最长肌转录组数据筛选差异表达lncRNAs和差异表达基因(DEGs)。根据差异表达lncRNAs和DEGs位置关系以及表达的相关性预测lncRNAs顺式作用(cis)和反式作用(trans)的靶DEGs,并对trans靶DEGs进行GO生物学过程和KEGG通路富集分析。进一步筛选潜在影响肌肉生长发育和脂肪沉积的候选差异表达lncRNAs和靶DEGs,构建lncRNAs-DEGs-生物学过程调控网络。最后将筛选的候选靶基因与藏猪和大约克夏猪的肉质性状(滴水损失、pH、肉色、大理石花纹、肌内脂肪含量、肌纤维横截面积、眼肌面积和维生素B1含量)进行关联分析。结果显示,在藏猪和大约克夏猪背最长肌组织中共鉴定出1 546个lncRNAs和20 219个mRNAs,其中包括49个差异表达lncRNAs和154个DEGs,差异表达lncRNAs分别顺式和反式调控7个和138个靶DEGs。功能富集分析发现,大量靶DEGs显著富集到与肌肉生长发育和脂肪沉积相关的生物学过程,包括磷脂酰肌醇3-激酶信号的正调控、肌肉结构发育、脂肪细胞分化、JAK-STAT级联的正调控、MAPK级联的正调控等。调控网络显示,MSTRG.34836.10、MSTRG.151.1调控的靶DEGs最多,其次是MSTRG.21406.1、MSTRG.26709.1。关联分析发现,这些lncRNAs调控的候选靶DEGs有11个与肉质性状显著相关。综上表明,本研究鉴定的部分差异表达lncRNAs与肌肉生长发育和脂肪沉积有关,分析发现这些差异表达lncRNAs可能通过调控HES1、LEPTGFB2和PER2等基因的表达影响猪肉品质,为进一步解析猪肉质性状的潜在分子机制奠定了基础。
关键词藏猪    大约克夏猪    RNA-seq    lncRNAs    肉质性状    
Screening and Functional Analysis of lncRNA Affecting Muscle Quality in Tibetan and Yorkshire Pigs
SONG Mingkun1, XUE Mingming1, ZHANG Lige1, YAN Duo1, XIA Ning1, SHANG Peng2, LI Xinjian1, HAN Xuelei1, QIAO Ruimin1, LI Xiuling1, LI Ming1, WANG Kejun1     
1. College of Animal Science and Technology, Henan Agricultural University, Zhengzhou 450002, China;
2. Animal Science College, Tibet Agricultural and Animal Husbandry University, Linzhi 860000, China
Abstract: This study aimed to identify the long non-coding RNAs (lncRNAs) and genes related to pork quality traits of the longissimus dorsi in Tibetan pigs and Yorkshire pigs (180 days old), and to further explore their molecular regulatory role, which could provide theoretical basis for the improvement breeding of pork quality traits. According to previously published data, 180-day-old healthy Tibetan boars (n=3) and Yorkshire boars (n=3) raised under the same conditions were used as the research objects to collect phenotype data of the quality of the longissimus dorsi and the transcriptome data of the longissimus dorsi of the same individual were used to screen differentially expressed lncRNAs and differentially expressed genes (DEGs). Cis- and trans-target DEGs were predicted based on the positional relationship and expression correlation of differentially expressed lncRNAs and DEGs, then trans-target DEGs were enriched by GO biological processes and KEGG pathway analysis. Subsequently, lncRNAs-DEGs-biological process regulatory network related to muscle growth and fat deposition was constructed. Then, association events between candidate target genes and meat quality traits were analyzed, which comprised of drip loss, pH, meat color, marbling, intramuscular fat content, muscle fiber cross-sectional area, loin muscle area and vitamin B1 content. A total of 1 546 lncRNAs and 20 219 mRNAs were identified in the longissimus dorsi of Tibetan pigs and Yorkshire pigs, there were 49 differentially expressed lncRNAs and 154 DEGs of them based on statistical analysis. Of these, 7 cis and 138 trans regulatory DEGs targeted by differentially expressed lncRNAs were predicted. Functional enrichment analysis showed that a large number of candidate DEGs were significantly enriched in biological processes of muscle growth and fat deposition, including positive regulation of phosphatidylinositol 3-kinase signaling, muscle structure development, fat cell differentiation, positive regulation of JAK-STAT cascade, and positive regulation of MAPK cascade, etc. Evidence from regulatory network showed that MSTRG.34836.10 and MSTRG.151.1 had the most target DEGs, followed by MSTRG.21406.1, MSTRG.26709.1. Association analysis indicated that 11 candidate target DEGs were significantly related to the meat quality traits. Taken together, some differentially expressed lncRNAs identified in this study were related to muscle growth and fat deposition. It was also suggested that these differentially expressed lncRNAs had an effect on pork quality by regulating the expression of HES1, LEP, TGFB2, and PER2. This study laid the foundation for further analysis of the potential molecular mechanism of pork quality traits.
Key words: Tibetan pig    Yorkshire pig    RNA-seq    lncRNAs    meat quality traits    

在生猪产业中,猪的肉质性状是重要的经济性状,主要包括肉色、滴水损失、大理石花纹、肌内脂肪含量等指标。猪肉品质与肌肉生长发育和脂肪沉积密切相关,在很大程度上直接影响消费者对猪肉的选择[1-2]。藏猪是我国宝贵的地方种质资源,主要起源于西藏高原。由于长期在恶劣的环境中经历自然选择[3-5],藏猪形成了耐粗饲、耐低氧、抗病力较强、脂肪沉积能力强等品种特性[6-7]。藏猪肉因其肉质鲜美,富含蛋白质和氨基酸等特点深受广大消费者的喜爱[8]。大约克夏猪在世界范围内广泛分布,是世界著名的瘦肉猪品种,具有生长速度快、瘦肉率高等优良特点。已有研究表明,藏猪和大约克夏猪在肌肉品质上存在较大差异[9],因此两者是比较猪肌肉品质的理想试验模型。

长链非编码RNA(long non-coding RNA,lncRNA)是一类长度大于200 bp,几乎不具有编码蛋白能力的非编码RNA。大量研究表明,lncRNAs在细胞增殖、分化、发育、疾病等生命活动中发挥着重要的调节功能[10-11]。随着高通量测序技术的不断发展与应用,发现lncRNAs在家养动物肌肉生长发育和脂肪沉积过程中发挥着重要作用[12-13]。目前已有一些关于不同猪种背最长肌转录组的比较分析。Wang等[14]通过对大约克夏猪和淮南猪胚胎期背最长肌的组织进行比较分析,发现大量差异表达的lncRNAs可能在肌肉生长发育和脂肪沉积过程中发挥着重要作用,并对候选差异表达lncRNA进行功能验证,发现IMFlnc1可抑制脂肪的生成。Li等[15]通过分析圩猪和大约克夏猪背最长肌的转录组数据,发现许多差异表达lncRNAs可能通过调节其潜在靶基因从而影响肌内脂肪的发育过程。然而,关于藏猪和大约克夏猪背最长肌组织转录组lnc- RNAs的比较分析还鲜有报道。

本研究基于藏猪和大约克夏猪的背最长肌组织转录组测序,对不同猪种中的lncRNAs和基因进行生物信息学分析,以确定影响藏猪和大约克夏猪肉品质的候选lncRNAs和基因。这些研究结果将有助于从肌肉生长发育和脂肪沉积的角度解析藏猪和大约克夏猪肉质性状形成的分子机制,为猪的遗传改良提供理论参考。

1 材料与方法 1.1 肌肉品质表型和RNA-seq数据的收集

本研究试验动物为180日龄健康藏公猪(n=3)和大约克夏公猪(n=3),均饲养于相同饲养环境,且保证饲养条件一致[16]。背最长肌肉品质测定数据来源于四川省畜牧科学院[16]。肉质指标主要包括滴水损失、pH45 min、pH24 h、肉色、大理石花纹、肌内脂肪含量、眼肌面积、肌纤维横截面积和维生素B1(VB1)含量。RNA-seq数据来源于同一样本藏猪和大约克夏猪的背最长肌组织,均收集于NCBI数据库(登录号:SRP090525)[16]

1.2 lncRNAs的鉴定流程

测序数据包含带接头和低质量的reads,故使用FastQC软件对测序数据进行质量控制。过滤后的高质量序列(clean reads)用于后续分析。使用HISAT2(v2.0.0-beta)软件将FASTQ格式的clean reads比对到猪的参考基因组Sus scrofa 11.1。

使用StringTie(v1.3.0)软件将比对后的数据用于转录本的组装。利用gffcompare(v0.10.1)软件将组装好的转录本与Ensemble基因集(Sscrofa 11.1.91)相比较,将具有“i”,“u”,“o”和“x”这类代码的转录本进行后续lncRNA的鉴定。为了避免不完全组装和假阳性率,对转录本以长度>200 bp和外显子数>2的标准进行过滤。使用CNCI[17]和CPC[18]软件对过滤后的转录本进行蛋白编码潜力的预测,CNCI和CPC得分均低于0的转录本被认为是lncRNAs。

1.3 lncRNAs和基因的差异表达分析

利用HTseq软件[19]计算比对到每个基因的reads数。使用DEseq2软件[20]鉴定差异表达lncRNAs和基因,筛选标准为差异表达倍数|log2(fold change)|≥1且显著性P-adj<0.05。

1.4 lncRNAs的靶DEGs预测及功能富集分析

lncRNAs的顺式调控作用被定义为对邻近基因的影响[21]。Bedtools软件[22]用于鉴定位于差异表达lncRNAs上游或下游100 kb的顺式调控靶DEGs。利用皮尔逊相关性检验计算lncRNAs和DEGs之间的表达相关性来预测差异表达lncRNAs反式调控的靶DEGs。

基于Metascape数据库[23]对差异表达lncRNAs的靶DEGs进行GO生物学过程(biological process)和KEGG信号通路功能富集分析。Cytoscape(v3.8.2)软件[24]用于构建与肌肉生长发育和脂肪沉积相关的候选差异表达lncRNAs、靶DEGs和生物学过程三者间的调控网络。

1.5 lncRNAs的靶DEGs与肉质性状的关联分析

基于R语言环境下使用“Corrplot”和“Psych”软件包对差异表达lncRNAs的靶DEGs的表达水平和肉质性状的表型数据进行关联分析。皮尔逊相关系数用于度量基因与表型间的相关程度。FDR校正法用于确定每个相关事件的显著性水平。

1.6 数据统计分析

使用EXCEL 2019对试验表型数据进行初步整理,利用SPSS 22.0统计软件进行独立样本t检验。结果使用“平均值±标准差”表示,P<0.05为差异显著判断标准,P<0.01为差异极显著判断标准。

2 结果 2.1 藏猪和大约克夏猪肌肉品质的比较

180日龄藏猪和大约克夏猪的背最长肌肌肉品质测定结果显示,大约克夏猪组的肌纤维横截面积和眼肌面积均极显著高于藏猪组(P<0.01);藏猪组的肌肉pH45 min和肉色均显著高于大约克夏猪组(P<0.05),且肌内脂肪含量、VB1含量极显著高于大约克夏猪组(P<0.01);而肌肉滴水损失、pH24 h和大理石花纹在两组间不存在显著性差异(P>0.05,图 1)。

pH45 min和pH24 h分别代表屠宰后45 min和24 h肌肉的pH。下图同 pH45 min and pH24 h represent the pH of the muscle at 45 min and 24 h after slaughtering, respectively. The same as below 图 1 藏猪和大约克夏猪肌肉品质的比较分析 Fig. 1 Comparative analysis of meat quality between Tibetan and Yorkshire pigs
2.2 转录组数据评估和lncRNAs的鉴定

本研究测序数据质量控制后共获得227 087 843个clean reads,包括藏猪组的114 520 669个和大约克夏猪组的112 567 174个。clean reads比对到猪参考基因组Sscrofa 11.1的比对率超过82%,其中超过85.5%的reads比对到CDS、5′UTR和3′UTR区域(图 2A)。Reads进行转录本组装后共获得21 765个转录本。对这些转录本进行鉴定分析,共筛选出20 219个mRNAs和1 546个lncRNAs。藏猪组和大约克夏猪组共有的mRNAs有15 329个,在藏猪组和大约克夏猪组中特异性表达的mRNAs分别有2 839和2 051个(图 2B);藏猪组和大约克夏猪组中共有的lncRNAs有1 314个,在藏猪组和大约克夏猪组中特异性表达的lncRNAs分别有145和87个(图 2C)。对lncRNAs进行分类,如图 2D所示,1 546个lncRNAs中基因间lncRNAs(lincRNAs)有654个,反义lncRNAs(lncNATs)有556个,内含子lncRNAs(ilncRNAs)有336个。

A. 6个样本的reads分布统计图: TA代表藏猪,YA代表大约克夏猪;B. 藏猪与大约克夏猪中鉴定出的mRNAs数量;C. 藏猪与大约克夏猪鉴定出的lncRNAs数量;D. lncRNAs的分类统计图: lincRNAs代表基因间lncRNAs,lncNATs代表反义lncRNAs,ilncRNAs代表内含子lncRNAs A. Reads distribution of 6 samples: TA is Tibetan pigs, and YA is Yorkshire pigs; B. The number of mRNAs identified in Tibetan and Yorkshire pigs; C. The number of lncRNAs identified in Tibetan and Yorkshire pigs; D. The classification and statistics of lncRNAs: lincRNAs represent intergenic lncRNAs, lncNATs represent antisense lncRNAs, ilncRNAs represent introns lncRNAs 图 2 lncRNAs和mRNAs的鉴定 Fig. 2 Identification of lncRNAs and mRNAs
2.3 差异表达lncRNAs和DEGs的筛选

根据|log2(fold change)|≥1且P-adj<0.05的标准筛选差异表达基因,结果表明在两组中共筛选出49个差异表达lncRNAs和154个DEGs。与藏猪组相比,大约克夏猪组中有28个上调lncRNAs和21个下调lncRNAs(图 3A),74个上调DEGs和80个下调DEGs(图 3B)。

红色代表上调基因;灰色代表表达无显著变化的基因;绿色代表下调基因 Red represents up-regulated genes; gray represents genes with no significant change in expression; green represents down-regulated genes 图 3 差异表达lncRNAs和DEGs的火山图 Fig. 3 Volcano plots of differentially expressed lncRNAs and DEGs
2.4 差异表达lncRNAs靶基因的预测

对差异表达lncRNAs进行顺式作用靶基因预测,共获得了8对具有潜在顺式调控作用的差异表达lncRNAs和靶DEGs(表 1)。其中MSTRG.42387.1上下游100 kb共存在3个潜在顺式作用靶DEGs,包括ATP合酶亚基8(ATP synthase F0 subunit 8,ATP8)、NADH脱氢酶亚单位1(NADH dehydrogenase subunit 1,ND1)和NADH脱氢酶亚单位2(NADH dehydrogenase subunit 1,ND2)。MSTRG.27311.5、MSTRG.42191.1、MSTRG.43545.1分别靶向ENSSSCG00000031583、β-防御素-1(beta defensin 1,DEFB1)和diaphanous相关成蛋白2(diaphanous related formin 2,DIAPH2)基因。MSTRG.8729.4和MSTRG.8729.8同时靶向肌球蛋白重链13(myosin heavy chain 13,MYH13)。

表 1 差异表达lncRNAs顺式作用靶DEGs Table 1 Differentially expressed target genes of 100 kb upstream and downstream of lncRNAs

根据皮尔逊相关性系数|correlation|≥0.95且P<0.05的筛选标准发现48个差异表达lncRNAs与138个DEGs存在493个反式作用关系对。在大约克夏猪组和藏猪组中,差异表达倍数前5的lncRNAs及其反式作用的靶DEGs如表 2所示。值得注意的是,不同的lncRNAs其潜在的反式作用靶DEGs的数量差别很大,MSTRG.42621.1反式作用调控25个DEGs,MSTRG.151.1反式作用调控18个DEGs,而MSTRG.8729.3和MSTRG.34836.20分别只有2个反式作用靶DEGs。

表 2 差异倍数最大的前5 lncRNAs及其反式差异表达靶基因 Table 2 Top 5 up-regulated and top 5 down-regulated lncRNAs with the greatest fold changes and their trans-differentially expressed target genes
2.5 差异表达lncRNAs反式差异表达靶基因的功能富集分析

为了揭示差异表达lncRNAs潜在的生物学功能,本研究对差异表达lncRNAs的反式作用DEGs进行GO生物学过程和KEGG通路功能富集分析,结果显示,靶DEGs显著富集到了循环系统过程(circulatory system process)、骨骼系统发育(skeletal system development)、心脏发育(heart development)、磷脂酰肌醇3-激酶信号的正调控(positive regulation of phosphatidylinositol 3-kinase signaling)、脂肪细胞分化(fat cell differentiation)等生物学过程和一个糖尿病并发症的AGE-RACE信号通路(AGE-RAGE signaling pathway in diabetic complications)(图 4)。为进一步挖掘与猪肌肉生长发育和脂肪沉积密切相关的候选差异表达lncRNAs和靶DEGs,本研究构建了差异表达lncRNAs及其反式作用DEGs共表达调控网络,其中共有39个差异表达lncRNAs与23个差异表达基因存在较强的共表达关系且与肌肉生长发育和脂肪沉积相关(图 5),这些靶DEGs显著富集在磷脂酰肌醇3-激酶信号的正调控、肌肉结构发育(muscle structure develop- ment)、脂肪细胞分化(fat cell differentiation)、JAK-STAT级联的正调控(positive regulation of JAK-STAT cascade)、MAPK级联的正调控(positive regulation of MAPK cascade)等生物学过程(表 3)。lncRNAs-DEGs-生物学过程调控网络显示, MSTRG.34836.10、MSTRG.151.1调控的靶DEGs最多,其次是MSTRG.21406.1和MSTRG.26709.1。与lncRNAs存在靶向关系最多的基因为TGFB2,其次是JAK2,这些靶基因共同被MSTRG.12095.1、MSTRG.27311.5、MSTRG.33150.6、MSTRG.42387.11、MSTRG.42387.18和MSTRG.42387.2所调控,并共同参与到磷脂酰肌醇3-激酶信号的正调控、MAPK级联的正调控和MAPK级联的调控的生物学过程中。

左下图所富集的聚类根据聚类的ID添加颜色。右上图代表聚类的P The enriched clusters in the left bottom figure were colored by the cluster ID. The color of clusters in the upper right figure represents the P-value of the clusters 图 4 差异表达lncRNAs反式作用靶DEGs的功能富集分析图 Fig. 4 Functional enrichment analysis of trans-acting target DEGs of differentially expressed lncRNAs
菱形代表差异表达lncRNAs,圆代表靶DEGs,正方形代表与肌肉生长和脂肪沉积相关的生物学过程,红色代表基因表达上调,绿色代表基因表达下调 The diamonds represent the differentially expressed lncRNAs, the circles represent the target DEGs, and the squares represent the biological processes associated with muscle development and fat deposition. The red edge represents up-regulated gene expression, and the green edge represents down-regulated gene expression 图 5 差异表达lncRNAs和与肌肉生长和脂肪沉积相关的生物学过程相关靶DEGs的调控网络 Fig. 5 Regulation network of differentially expressed lncRNAs and target DEGs enriched in biological processes related to muscle growth and fat deposition
表 3 差异表达靶基因显著富集的与肌肉生长发育和脂肪沉积相关的GO生物学过程 Table 3 Significantly enriched GO biological processes related to muscle growth and fat deposition in DE target genes
2.6 靶DEGs与肉质性状的关联分析

为了更好地解析差异表达lncRNAs反式作用靶DEGs的潜在功能,本研究将上述筛选的23个与肌肉生长发育和脂肪沉积生物学过程相关的靶DEGs与藏猪和大约克夏猪的肉质性状进行关联分析。结果表明,有11个差异表达靶基因与肉质性状显著相关(图 6)。CCN1(cellular communication network factor 1)与滴水损失显著相关(P<0.05);肽酶抑制蛋白16(peptidase inhibitor 16,PI16)和原癌基因FOS(Fos proto-oncogene,FOS)与肌肉pH24 h显著相关(P<0.05);popeye结构域2(popeye domain containing 2,POPDC2)、PI16、FOS和生物钟基因PER2(period circadian regulator 2)共4个基因与肉色(CLR)显著相关(P<0.05);PI16、FOS、钙调磷酸酶B同源蛋白2(RB1 inducible coiled-coil 1,RB1CC1)、腱蛋白(chordin,CHRND)和受体活性修饰蛋白3(receptor activity modifying protein 3,RAMP3) 共5个基因与大理石花纹(MBL)显著相关(P<0.05);POPDC2PER2和P-选择素(selectin P,SELP)共3个基因与肌内脂肪含量显著相关(P<0.05);转化生长因子β2(transforming growth factor beta 2,TGFB2)和Janus激酶2基因(Janus kinase 2,JAK2)与眼肌面积(LMA)显著相关; POPDC2、PI16、FOSPER2共4个基因与维生素B1含量(VB1)显著相关(P<0.05);关联分析没有发现任何基因与pH45 min和肌纤维横截面积(CSAF)显著相关(P>0.05)。有趣的是,PI16和FOS基因与pH24 h、CLR、MBL和VB1含量均显著相关(P<0.05);POPDC2和PER2基因与CLR、IMF和VB1含量均显著相关(P<0.05)。

图中的肉质性状均使用原始数据除以对应个体胴体重的方法进行了归一化校正。圆颜色和圆大小代表基因表达水平和肉质性状的相关系数,即圆颜色越深和圆越大代表相关系数越高,蓝色代表正相关,红色代表负相关。采用FDR校正法进行显著性检验。*. P<0.05 The meat quality traits in the figure are normalized and corrected by dividing the original data by the corresponding individual's carcass weight. The color and size of the circles represent the correlation coefficients between gene expression levels and meat quality traits, that is, the darker the circle color and the larger the circle represent the higher the correlation coefficient, blue represents positive correlation, and red represents negative correlation. FDR correction method was used for significance test. *. P < 0.05 图 6 lncRNAs的反式差异表达靶基因的表达水平与肉质性状的关联分析 Fig. 6 Correlation analysis between expression of trans differentially expressed target genes of lncRNAs and meat quality traits
3 讨论

肉质性状属于微效多基因控制的数量性状,了解影响肉质性状的潜在分子机制对提高猪肉产量和肉质改良具有重要意义。先前的研究表明,哺乳动物基因组中含有大量的lncRNAs,一些lncRNAs已被发现在动物生长发育等生命过程中发挥着非常重要的调节作用[12-13]。目前,lncRNAs的研究主要集中在人和小鼠等模式动物[25-26],在猪上的研究相对较少,对猪肌肉品质相关lncRNAs的认识更是不足。藏猪和大约克夏猪的肌肉品质存在较大差异, 是研究猪肌肉品质理想的试验模型[27]。与藏猪组相比,大约克夏猪组具有较高的眼肌面积和较低的肌内脂肪含量。本研究通过分析藏猪和大约克夏猪的背最长肌组织的转录组测序数据,筛选潜在影响肌肉生长发育和脂肪沉积相关的候选lncRNAs和候选基因,以探究lncRNAs在藏猪和大约克夏猪肌肉品质调控中的作用。

3.1 差异表达lncRNA可能通过顺式作用参与猪肌肉生长发育过程

lncRNAs一般不具有编码蛋白的功能,可通过顺式调控模式直接结合到临近位置的DNA调节基因表达或以反式调控的模式调控lncRNA远端位置的基因表达[28]。Chen等[29]通过识别肌肉特异性表达lncRNAs临近位置的蛋白质编码基因(<100 kb)以预测lncRNA的功能。本研究通过预测差异表达lncRNAs的顺式调控靶基因,发现共有8对具有潜在顺式调控作用的差异表达lncRNAs-靶DEGs对,其中lncRNA MSTRG.43545.1与顺式调控靶基因DIAPH2存在较强的相关性。与藏猪相比,大约克夏猪的MSTRG.43545.1和DIAPH2的表达水平均显著提高。有研究报道,DIAPH2在果蝇成肌细胞融合过程中起着至关重要的作用,可诱导成肌细胞融合成核并促进肌动蛋白丝的延长[30]。MSTRG.43545.1可能通过上调DIAPH2的表达促进肌肉的生长发育过程。

3.2 差异表达lncRNAs可能通过反式作用调控猪肌肉生长和肌内脂肪沉积过程

Li等[31]通过lncRNA-mRNA共表达关系分析了lncRNAs在骨骼肌纤维形成的反式调控作用。在本研究中,共确定了48个差异表达lncRNAs与138个DEGs存在493个反式作用关系对。为了明确这些差异表达lncRNAs的生物学功能,本研究对其反式作用的靶DEGs进行富集分析,结果表明大量靶基因显著富集到与肌肉生长发育和脂肪沉积相关的生物学过程。其中Notch信号通路已被证明在骨骼肌的发育和再生过程中起着关键作用,在骨骼肌卫星细胞再生期间激活Notch信号通路可使肌肉干细胞保持静息状态[32]。此外,Notch信号还参与了调节白细胞和棕色脂肪细胞祖细胞的增殖和脂肪分化的过程[33]。磷脂酰肌醇3-激酶信号通路作为经典的胰岛素信号通路,在调节脂肪代谢中起着关键作用[34],还可以通过激活AKT、mTOR等下游靶蛋白参与到骨骼肌生长发育的过程中[35]。MAPK级联[36-37]和JAK-STAT级联[38-39]也均已被证明在肌肉生长发育和脂肪沉积过程中发挥着重要作用。这些生物学过程可能是导致藏猪和大约克夏猪肌肉生长发育和脂肪沉积性状差异显著的重要因素。

为进一步挖掘潜在影响肌肉生长发育和脂肪沉积的基因的互作关系,本研究构建了lncRNAs-DEGs-生物学过程调控网络,并将这些DEGs与藏猪和大约克夏猪的肉质性状进行关联分析。调控网络显示,MSTRG.12095.1、MSTRG.14538.2、MSTRG.25815.1、MSTRG.26232.1等共11个lncRNAs对TGFB2存在反式调控作用。TGFB2作为TGF-β/SMAD信号通路的成员之一,具有抑制骨骼肌卫星细胞增殖和分化,促进凋亡的作用[40]。本研究中,TGFB2与眼肌面积显著相关,且显著富集到Notch信号通路、肌肉结构发育等肌肉发育相关生物学过程。由此推测,在大量lncRNAs的调控作用下TGFB2可能通过多种途径共同参与到肌肉的生长发育过程中。在肌肉发育和再生过程中,HES1作为Notch信号通路下游的重要靶基因可通过抑制MYODMYOG的表达控制激活的肌肉干细胞增殖和分化之间的平衡[41]。与藏猪相比,MSTRG.151.1、MSTRG.21406.1、MSTRG.26709.1和HES1在大约克夏猪的表达显著下调,推测MSTRG.151.1、MSTRG.21406.1和MSTRG.26709.1可能通过正向调控HES1的表达进而参与肌肉的生长发育过程。

差异表达lncRNAs调节的靶DEGs还通过参与脂肪沉积的相关生物学过程导致藏猪与大约克夏猪肌肉品质差异。与藏猪相比,MSTRG.42191.1、MSTRG.11899.1和 LEP在大约克夏猪中显著下调,且MSTRG.42191.1是差异倍数第二大的lncRNAs。LEP基因多态性与猪[42]和绵羊[43]的脂肪特征相关。LEP基因是白色脂肪的标志基因,调节脂肪沉积的重要因子,可促进脂肪的生成和脂滴的形成[44]LEP基因被显著富集到了磷脂酰肌醇3-激酶信号的正调控、JAK-STAT级联和MAPK级联的调控和脂肪细胞分化等GO生物学过程中。据此推测,MSTRG.42191.1和MSTRG.11899.1可能通过靶向调控LEP基因的表达促进藏猪的肌内脂肪沉积,MSTRG.42191.1可能在肌内脂肪沉积中发挥着更为重要的作用。PER2是一个节律基因,可通过直接调节PPARγ基因影响脂肪细胞的分化,PER2基因缺陷的小鼠表现出脂质代谢的改变以及总三酰甘油和非酯化脂肪酸急剧减少[45]。本研究中,PER2与肌内脂肪含量显著相关,并显著富集到了脂肪细胞分化生物学过程。与藏猪相比,MSTRG.3334.1、MSTRG.8729.8在大约克夏猪中显著上调,而PER2在大约克夏猪中显著下调,推测MSTRG.3334.1、MSTRG.8729.8可能通过下调PER2的表达进而阻碍脂肪细胞分化过程,从而导致了藏猪和大约克夏猪的肌内脂肪含量的显著差异。本研究筛选到许多与肌肉生长发育和脂肪沉积相关的新lncRNAs,为深入解析肌肉品质的形成提供了重要基础依据,其与靶基因的调控作用和具体的作用机制还有待进一步探究。

4 结论

藏猪和大约克夏猪背最长肌中部分lncRNAs和基因存在显著性差异,推测这些lncRNAs在肌肉品质的调节中发挥重要作用。MSTRG.12095.1、MSTRG.151.1、MSTRG.42191.1、MSTRG.3334.1等一些lncRNAs可能通过调节TGFB2、HES1、LEPPER2等基因的表达来影响肌肉生长发育和肌内脂肪沉积的过程,从而导致藏猪和大约克夏猪在肌肉生长发育和肌内脂肪含量方面的差异。本研究从lncRNA的角度为猪肉质改良工作中候选功能基因的研究提供了理论依据。

参考文献
[1]
ZAPPATERRA M, GIOIOSA S, CHILLEMI G, et al. Dissecting the gene expression networks associated with variations in the major components of the fatty acid Semimembranosus muscle profile in Large White heavy pigs[J]. Animals (Basel), 2021, 11(3): 628.
[2]
ZHANG J, CUI L, MA J, et al. Transcriptome analyses reveal genes and pathways associated with fatty acid composition traits in pigs[J]. Anim Genet, 2017, 48(6): 645-652. DOI:10.1111/age.12597
[3]
LI M Z, TIAN S L, JIN L, et al. Genomic analyses identify distinct patterns of selection in domesticated pigs and Tibetan wild boars[J]. Nat Genet, 2013, 45(12): 1431-1438. DOI:10.1038/ng.2811
[4]
SHANG P, LI W T, LIU G, et al. Identification of lncRNAs and genes responsible for fatness and fatty acid composition traits between the Tibetan and Yorkshire pigs[J]. Int J Genomics, 2019, 2019: 5070975.
[5]
SHANG P, LI W T, TAN Z K, et al. Population genetic analysis of ten geographically isolated Tibetan pig populations[J]. Animals (Basel), 2020, 10(8): 1297.
[6]
MA Y F, HAN X M, HUANG C P, et al. Population genomics analysis revealed origin and high-altitude adaptation of Tibetan pigs[J]. Sci Rep, 2019, 9(1): 11463. DOI:10.1038/s41598-019-47711-6
[7]
李梦柔, 韦明邦, 张健, 等. 基于iTRAQ技术挖掘猪肾组织高原低氧适应性关键蛋白[J]. 畜牧兽医学报, 2021, 52(5): 1238-1246.
LI M R, WEI M B, ZHANG J, et al. Mining key proteins of high altitude hypoxic adaptability in porcine kidney tissue based on iTRAQ technology[J]. Acta Veterinaria et Zootechnica Sinica, 2021, 52(5): 1238-1246. (in Chinese)
[8]
MI S, LI X, ZHANG C H, et al. Characterization and discrimination of Tibetan and Duroc×(Landrace×Yorkshire) pork using label-free quantitative proteomics analysis[J]. Food Res Int, 2019, 119: 426-435. DOI:10.1016/j.foodres.2019.02.016
[9]
GAN M L, SHEN L Y, FAN Y, et al. High altitude adaptability and meat quality in Tibetan pigs: a reference for local pork processing and genetic improvement[J]. Animals (Basel), 2019, 9(12): 1080.
[10]
STATELLO L, GUO C J, CHEN L L, et al. Gene regulation by long non-coding RNAs and its biological functions[J]. Nat Rev Mol Cell Biol, 2021, 22(2): 96-118. DOI:10.1038/s41580-020-00315-9
[11]
YANG S M, LIM K H, KIM S H, et al. Molecular landscape of long noncoding RNAs in brain disorders[J]. Mol Psychiatry, 2021, 26(4): 1060-1074. DOI:10.1038/s41380-020-00947-5
[12]
WANG L J, YANG X, ZHU Y H, et al. Genome-wide identification and characterization of long noncoding RNAs of brown to white adipose tissue transformation in goats[J]. Cells, 2019, 8(8): 904. DOI:10.3390/cells8080904
[13]
CHEN R, LEI S, JIANG T, et al. Roles of lncRNAs and circRNAs in regulating skeletal muscle development[J]. Acta Physiol (Oxf), 2020, 228(2): e13356.
[14]
WANG J, CHEN M Y, CHEN J F, et al. LncRNA IMFlnc1 promotes porcine intramuscular adipocyte adipogenesis by sponging miR-199a-5p to up-regulate CAV-1[J]. BMC Mol Cell Biol, 2020, 21(1): 77. DOI:10.1186/s12860-020-00324-8
[15]
LI Q Q, HUANG Z Y, ZHAO W J, et al. Transcriptome analysis reveals long intergenic non-coding RNAs contributed to intramuscular fat content differences between Yorkshire and Wei Pigs[J]. Int J Mol Sci, 2020, 21(5): 1732. DOI:10.3390/ijms21051732
[16]
TAO X, LIANG Y, YANG X M, et al. Transcriptomic profiling in muscle and adipose tissue identifies genes related to growth and lipid deposition[J]. PLoS One, 2017, 12(9): e0184120. DOI:10.1371/journal.pone.0184120
[17]
SUN L, LUO H T, BU D C, et al. Utilizing sequence intrinsic composition to classify protein-coding and long non-coding transcripts[J]. Nucleic Acids Res, 2013, 41(17): e166. DOI:10.1093/nar/gkt646
[18]
KONG L, ZHANG Y, YE Z Q, et al. CPC: assess the protein-coding potential of transcripts using sequence features and support vector machine[J]. Nucleic Acids Res, 2007, 35(Web Server issue): W345-W349.
[19]
ANDERS S, PYL P T, HUBER W. HTSeq—a Python framework to work with high-throughput sequencing data[J]. Bioinformatics, 2015, 31(2): 166-169. DOI:10.1093/bioinformatics/btu638
[20]
LOVE M I, HUBER W, ANDERS S. Moderated estimation of fold change and dispersion for RNA-seq data with DESeq2[J]. Genome Biol, 2014, 15(12): 550. DOI:10.1186/s13059-014-0550-8
[21]
GUIL S, ESTELLER M. Cis-acting noncoding RNAs: friends and foes[J]. Nat Struct Mol Biol, 2012, 19(11): 1068-1075. DOI:10.1038/nsmb.2428
[22]
QUINLAN A R. BEDTools: the Swiss-army tool for genome feature analysis[J]. Curr Protoc Bioinformatics, 2014, 47: 11.12.1-11.12.34.
[23]
ZHOU Y Y, ZHOU B, PACHE L, et al. Metascape provides a biologist-oriented resource for the analysis of systems-level datasets[J]. Nat Commun, 2019, 10(1): 1523. DOI:10.1038/s41467-019-09234-6
[24]
SAITO R, SMOOT M E, ONO K, et al. A travel guide to Cytoscape plugins[J]. Nat Methods, 2012, 9(11): 1069-1076. DOI:10.1038/nmeth.2212
[25]
KAZIMIERCZYK M, KASPROWICZ M K, KASPRZYK M E, et al. Human long noncoding RNA interactome: detection, characterization and function[J]. Int J Mol Sci, 2020, 21(3): 1027. DOI:10.3390/ijms21031027
[26]
JOSHI M, RAJENDER S. Long non-coding RNAs (lncRNAs) in spermatogenesis and male infertility[J]. Reprod Biol Endocrinol, 2020, 18(1): 103. DOI:10.1186/s12958-020-00660-6
[27]
SHANG P, WANG Z X, CHAMBA Y, et al. A comparison of prenatal muscle transcriptome and proteome profiles between pigs with divergent growth phenotypes[J]. J Cell Biochem, 2019, 120(4): 5277-5286. DOI:10.1002/jcb.27802
[28]
TSAI M C, MANOR O, WAN Y, et al. Long noncoding RNA as modular scaffold of histone modification complexes[J]. Science, 2010, 329(5992): 689-693. DOI:10.1126/science.1192002
[29]
CHEN G T, CHENG X F, SHI G L, et al. Transcriptome analysis reveals the effect of long intergenic noncoding RNAs on pig muscle growth and fat deposition[J]. Biomed Res Int, 2019, 2019: 2951427.
[30]
DENG S, BOTHE I, BAYLIES M K. The formin diaphanous regulates myoblast fusion through actin polymerization and Arp2/3 regulation[J]. PLoS Genet, 2015, 11(8): e1005381. DOI:10.1371/journal.pgen.1005381
[31]
LI R Y, LI B J, JIANG A W, et al. Exploring the lncRNAs related to skeletal muscle fiber types and meat quality traits in pigs[J]. Genes (Basel), 2020, 11(8): 883. DOI:10.3390/genes11080883
[32]
GOPINATH S D, WEBB A E, BRUNET A, et al. FOXO3 promotes quiescence in adult muscle stem cells during the process of self-renewal[J]. Stem Cell Reports, 2014, 2(4): 414-426. DOI:10.1016/j.stemcr.2014.02.002
[33]
SHAN T Z, LIU J Q, WU W C, et al. Roles of Notch signaling in adipocyte progenitor cells and mature adipocytes[J]. J Cell Physiol, 2017, 232(6): 1258-1261. DOI:10.1002/jcp.25697
[34]
BERETTA M, BAUER M, HIRSCH E. PI3K signaling in the pathogenesis of obesity: the cause and the cure[J]. Adv Biol Regul, 2015, 58: 1-15. DOI:10.1016/j.jbior.2014.11.004
[35]
YOSHIDA T, DELAFONTAINE P. Mechanisms of IGF-1-mediated regulation of skeletal muscle hypertrophy and atrophy[J]. Cells, 2020, 9(9): 1970. DOI:10.3390/cells9091970
[36]
SIN T K, ZHANG G H, ZHANG Z C, et al. Cancer-induced muscle wasting requires p38β MAPK activation of p300[J]. Cancer Res, 2021, 81(4): 885-897. DOI:10.1158/0008-5472.CAN-19-3219
[37]
DA SILVA C, DURANDT C, KALLMEYER K, et al. The role of Pref-1 during adipogenic differentiation: an overview of suggested mechanisms[J]. Int J Mol Sci, 2020, 21(11): 4104. DOI:10.3390/ijms21114104
[38]
COLLOTTA D, HULL W, MASTROCOLA R, et al. Baricitinib counteracts metaflammation, thus protecting against diet-induced metabolic abnormalities in mice[J]. Mol Metab, 2020, 39: 101009. DOI:10.1016/j.molmet.2020.101009
[39]
ADDINSALL A B, CACCIANI N, AKKAD H, et al. JAK/STAT inhibition augments soleus muscle function in a rat model of critical illness myopathy via regulation of complement C3/3R[J]. J Physiol, 2021, 599(11): 2869-2886. DOI:10.1113/JP281220
[40]
YIN H D, HE H R, SHEN X X, et al. MicroRNA profiling reveals an abundant miR-200a-3p promotes skeletal muscle satellite cell development by targeting TGF-β2 and regulating the TGF-β2/SMAD signaling pathway[J]. Int J Mol Sci, 2020, 21(9): 3274. DOI:10.3390/ijms21093274
[41]
LAHMANN I, BRÖHL D, ZYRIANOVA T, et al. Oscillations of MyoD and Hes1 proteins regulate the maintenance of activated muscle stem cells[J]. Genes Dev, 2019, 33(9-10): 524-535. DOI:10.1101/gad.322818.118
[42]
PÉREZ-MONTARELO D, RODRÍGUEZ M C, FERNÁNDEZ A, et al. Haplotypic diversity of porcine LEP and LEPR genes involved in growth and fatness regulation[J]. J Appl Genet, 2015, 56(4): 525-533. DOI:10.1007/s13353-015-0284-7
[43]
HAJIHOSSEINLO A, JAFARI S, AJDARY M. The relationship of GH and LEP gene polymorphisms with fat-tail measurements (fat-tail dimensions) in fat-tailed Makooei breed of Iranian sheep[J]. Adv Biomed Res, 2015, 4: 172.
[44]
DARDENO T A, CHOU S H, MOON H S, et al. Leptin in human physiology and therapeutics[J]. Front Neuroendocrinol, 2010, 31(3): 377-393. DOI:10.1016/j.yfrne.2010.06.002
[45]
GRIMALDI B, BELLET M M, KATADA S, et al. PER2 controls lipid metabolism by direct regulation of PPARγ[J]. Cell Metab, 2010, 12(5): 509-520. DOI:10.1016/j.cmet.2010.10.005

(编辑   郭云雁)