影响藏猪与大约克夏猪肌肉品质的lncRNA的筛选及功能分析

宋明坤，薛明明，张力戈，颜 铎，夏 宁，商 鹏，李新建，韩雪蕾，乔瑞敏，李秀领，李 明，王克君*

(1. 河南农业大学动物科技学院，郑州 450002； 2. 西藏农牧学院动物科学学院，林芝 860000)

在生猪产业中，猪的肉质性状是重要的经济性状，主要包括肉色、滴水损失、大理石花纹、肌内脂肪含量等指标。猪肉品质与肌肉生长发育和脂肪沉积密切相关，在很大程度上直接影响消费者对猪肉的选择[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个。

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)。

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)。

根据皮尔逊相关性系数|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。

A. 6个样本的reads分布统计图:TA代表藏猪，YA代表大约克夏猪；B. 藏猪与大约克夏猪中鉴定出的mRNAs数量；C. 藏猪与大约克夏猪鉴定出的lncRNAs数量；D. lncRNAs的分类统计图:lincRNAs代表基因间lncRNAs，lncNATs代表反义lncRNAs，ilncRNAs代表内含子lncRNAsA. 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

红色代表上调基因；灰色代表表达无显著变化的基因；绿色代表下调基因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.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级联的调控的生物学过程中。

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)；POPDC2、PER2和P-选择素(selectin P，SELP)共3个基因与肌内脂肪含量显著相关(P<0.05)；转化生长因子β2(transforming growth factor beta 2，TGFB2)和Janus激酶2基因(Janus kinase 2，JAK2)与眼肌面积(LMA)显著相关;POPDC2、PI16、FOS和PER2共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)。

3 讨 论

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

左下图所富集的聚类根据聚类的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.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信号通路下游的重要靶基因可通过抑制MYOD和MYOG的表达控制激活的肌肉干细胞增殖和分化之间的平衡[41]。与藏猪相比，MSTRG.151.1、MSTRG.21406.1、MSTRG.26709.1和HES1在大约克夏猪的表达显著下调，推测MSTRG.151.1、MSTRG.21406.1和MSTRG.26709.1可能通过正向调控HES1的表达进而参与肌肉的生长发育过程。

表3 差异表达靶基因显著富集的与肌肉生长发育和脂肪沉积相关的GO生物学过程

差异表达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，为深入解析肌肉品质的形成提供了重要基础依据，其与靶基因的调控作用和具体的作用机制还有待进一步探究。

图中的肉质性状均使用原始数据除以对应个体胴体重的方法进行了归一化校正。圆颜色和圆大小代表基因表达水平和肉质性状的相关系数，即圆颜色越深和圆越大代表相关系数越高，蓝色代表正相关，红色代表负相关。采用FDR校正法进行显著性检验。*. P<0.05The 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

4 结 论

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