家猪长链非编码RNA的组织特异性表达及其对肉质的影响

王浩杰， 张珍华， 李琰， 张霞， 彭锐

（四川大学生命科学学院，生物资源与生态环境教育部重点实验室，成都610065）

哺乳动物的基因组中约2/3可转录，而仅约1.9%的部分可合成蛋白（Mattick，2001；Djebaliet al.，2012）。多数非编码基因的转录产物称为非编码RNA（non-coding RNA，ncRNA），其中长度超过200 nt的RNA被定义为长链非编码RNA（long non-coding RNA，lncRNA）。与 mRNA相同的是，lncRNA也是由RNA聚合酶Ⅱ转录，具有5’-和3’-poly A尾部的结构（Mercer & Mattick，2013；Marcheseet al.，2017）。除了无法编码蛋白外，与mRNA相比，lncRNA有表达量更低、更高的组织特异性表达模式等的特点（Mongelliet al.，2019）。LncRNA曾被视为转录组的“噪音”，但它可通过与蛋白或蛋白编码基因（protein-coding gene，PCG）等生物大分子作用，从而参与一系列重要的生物过程，如转录调控、细胞分化和细胞自噬等（Gil & Ulitsky，2020；Sunet al.，2020）。

在包括人类和家畜在内的多个物种中，lncRNA都可呈现出组织特异性的表达特点（Tsoiet al.，2015；Kernet al.，2018；Cavaet al.，2019）。一些组织特异性lncRNA（tissue-specific lncRNA，TS lncRNA）已被证实对其所在的组织发育有着重要的影响，Grote等（2013）在小鼠Mus musculus的胚胎细胞中发现了特异性表达于侧中胚层的TS lncRNA Fendrr，并证实对Fendrr的沉默会导致小鼠的心脏和体壁发育异常。Morán等（2012）筛选了人胰岛内的TS lncRNA，并发现它们是胰岛β细胞分化过程中的重要组成，还发现了lncRNA HILNC25对糖尿病相关基因GLIS3的正调控作用。

猪肉不仅在肉类消费市场上占有重要地位，也是人类摄入蛋白的重要来源。因此，提升猪肉的品质对生猪养殖业十分重要。肌肉内脂肪（intramuscular fat，IMF）是影响猪肉品质的重要因素，较高的IMF含量可以提升猪肉的多汁度、柔软度和风味（Daszkiewiczet al.，2005）。影响IMF含量的因素包括品种、年龄、饲养环境、饲料中蛋白质含量等（Parket al.，2018；Malgwiet al.，2022），然 而lncRNA-mRNA的相互作用对IMF的影响还不清楚。

本研究利用公共数据库中的RNA-seq数据，重新组装了家猪Sus scrofa domestica的转录组并鉴定了lncRNA，筛选了在肌肉、脂肪、脾脏、肾脏、肝脏和卵巢中的TS lncRNA并预测了其功能。此外，还鉴定了影响其IMF含量的lncRNA并建立了相关的调控网络。这些数据为生猪养殖业和未来lncRNA的研究提供了参考。

1 材料与方法

1.1 数据来源

参考基因组（ftp：//ftp.ensembl.org/pub/release-98/fasta/sus_scrofa/）和注释文件（ftp：//ftp.ensembl.org/pub/release-98/gtf/sus_scrofa/）均 来 自 Ensembl数据库。RNA-seq数据均下载于基因表达（Gene Expression Omnibus，GEO）数据库，包括6个组织或器官（肌肉、脂肪、脾脏、肝脏、肾脏和卵巢）共36组用于 lncRNA 鉴定（Fushanet al.，2015；Draget al.，2017；Liet al.，2017；Limet al.，2017；TangQZet al.，2017；TangZLet al.，2017；Cheet al.，2018；Kimet al.，2018；Oczkowiczet al.，2018；Liuet al.，2019；Wanget al.，2019）。3组高IMF含量和3组低IMF含量的转录组测序数据（Limet al.，2017）用于差异表达分析。

表1 数据集信息汇总Table 1 The summary of datasets

1.2 转录组组装

RNA-seq 由 Trim Galore 0.6.4（Martin，2011）去除接头和低质量的测序片段后，使用Hisat2 2.1.0（Kimet al.，2019）将数据比对到参考基因组上并得到sam格式的比对文件。最后，使用StringTie 2.0（Perteaet al.，2015）将比对的结果进行组装并定量，合并后得到gtf格式的转录组注释文件。

1.3 LncRNA的鉴定

首先将注释文件中被标记为“lncRNA”的转录本视作已知的lncRNA。然后使用Cuffcompare 2.2.1工具将所有新鉴定的转录本进行分类，获取其中被标记为“i”“j”“u”“x”和“o”5 种类型之一的转录本，并过滤掉长度小于200 nt以及外显子数量＜2的转录本。随后，分别使用CNCI（Sunet al.，2013）、PLEK 1.2（Liet al.，2014）和 Pfam（Batemanet al.，2004）对上一步所得转录本进行编码能力预测，获取在3种方法中都显示无法编码蛋白的转录本。

1.4 TS lncRNA的鉴定

使用 Tspex 0.6.1（https：//tspex.lge.ibi.unicamp.br/）计算每个lncRNA的组织特异性指数（tissue specific index，TSI），TSIi=，式中，N为组织的总数，xi为转录本x在第i个组织中经归一化的表达量（Julienet al.，2012）。TSI的值在0～1之间，越接近1表明组织特异性程度越高。将TSI＞0.9且在所有样本中的平均表达量＞0.1 FPKM的lncRNA视作具有组织特异性。

1.5 加权基因共表达网络分析

使 用 R 包 WGCNA（Langfelder & Horvath，2008）构建加权基因共表达网络分析（weighted gene co-expression network analysis，WGCNA）。首先，根据表达量计算每对转录本之间的Pearson系数，选择5为软阈值并构建邻接矩阵。其次，根据邻接矩阵计算拓扑重叠度并得到拓扑重叠矩阵，再将其转换为相异矩阵并据其构建层次聚类树，通过动态树切割将表达模式相似的转录本聚集到同一模块。随后，计算出每个组织与模块之间的相关性，值越接近1或-1表示该组织与模块之间的相关性越高。最后，将与每个TS lncRNA相关性最高的10个mRNA视作其潜在的反式作用位点。

1.6 功能富集分析

首先使用R包的BiomaRt 2.42.1（Smedleyet al.，2009）将每个转录本转换为与之对应的Entrez id，随后使用 R 包的 ClusterProfiler 3.14.0（Yuet al.，2012）进行基因本体GO和KEGG富集分析，结果中P-value≤0.05的注释条目视为有效结果。

1.7 影响IMF的lncRNA-mRNA筛选

首先利用 R 包的 DESeq2 1.26.0（Loveet al.，2014）对不同IMF水平的RNA-seq进行差异分析。结果中log2（Fold Change）≥1且FDR≤0.05的mRNA被视作显著差异，而log2（Fold Change）≥1且P-value≤0.05的lncRNA被视作显著差异。其次，计算所有差异表达的lncRNA和mRNA之间的Pearson相关系数，选择Pearson系数绝对值＞0.9的转录本作为共表达的lncRNA-mRNA作用对。随后，利用Lnc-Tar（Liet al.，2015）估算这些lncRNA-mRNA之间的归一化结合自由能（ndG），以-0.13作为阈值筛选出可能相互作用的lncRNA-mRNA。最后进行功能富集分析，将与脂质代谢相关的lncRNA-mRNA视作潜在影响IMF的lncRNA-mRNA作用对。

1.8 LncRNA的二级结构预测

利用 ViennaRNA 2.4.18（Lorenzet al.，2011）基于成环结构以及动态规划算法预测目标lncRNA的二级结构：最小自由能结构和质心二级结构。

1.9 RNA干扰与过表达实验

首先将猪PK-15细胞培养于DEME（高糖）培养基，置于37 ℃，5%CO2环境。在PLKO.1质粒中插入shRNA构建目标lncRNA的干扰载体，对照组质粒中插入scramble shRNA。过表达实验中，将目标lncRNA的全长cDNA插入pcDNA-3.1质粒构建过表达载体，对照组使用空质粒。均使用High-Gene转染试剂将载体转染入细胞。

1.10 实时荧光定量PCR（RT-qPCR）

总RNA由TRIzol试剂提取并反转录为cDNA，随后使用 2×Taq SYBR®Green qPCR Premix对MSTRG.7256.5、PPARγ、SCD、AQP7和SORBS1进行RT-qPCR实验，选用管家基因GADPH为内参，每个样本3组重复，采用2-ΔΔCt法计算每个基因的相对表达水平（表2）。

表2 实验所用引物Table 2 Primers used in this study

1.11 油红O染色实验

首先将收集的细胞用PBS缓冲液清洗3次，再用10%甲醛溶液固定5 min，使用60%异丙醇溶液清洗5 min后，将细胞用油红O染液染色10 min，最后将细胞用蒸馏水洗涤并置于显微镜下观察。

2 结果

2.1 LncRNA的鉴定与分析

共获得23 678条lncRNA，其中，新鉴定的lncRNA有14 309条，包括3 061条反义lncRNA，3 051条基因间lncRNA，2 212条内含子lncRNA和5 985条正义lncRNA。与mRNA相比，lncRNA的长度更短，且这一点在新鉴定的lncRNA中更显著（图1：A）。在每条染色体上，mRNA的数量都明显多于lncRNA，仅在第10条、第11条、第16条和Y染色体上，二者的数量差距不大（图1：B）。此外，lncRNA还有着更少的外显子数量和更低的表达量（图1：C，D）。

图1 LncRNA的特征Fig. 1 Characterization of lncRNAs

2.2 组织特异性lncRNA的鉴定与分析

通过计算每条lncRNA的TSI指数，共有1 218条被鉴定为具有组织特异性，其中，卵巢的TS lncRNA数量最多（429条），而肌肉的最少（85条）（图2：A）。将TS lncRNA的位置映射到染色体上，第1条上的数量最多，Y染色体上为0，TS lncRNA的数量基本与其所在染色体的长度呈正相关关系（r=4.25×10-7，P=4.62×10-6）（图2：B，C）。TS lncRNA的表达量经归一化后的中位数均＜1（图2：D）。尽管卵巢中TS lncRNA的数量最多，但表达水平最低（图2：A，D）。

图2 TS lncRNA的特征Fig. 2 Characteristics of TS lncRNAs

2.3 TS lncRNA的功能预测

首先，查找TS lncRNA上、下游50 kb内的mRNA，将其视作顺式作用的潜在位点，并进行GO和KEGG富集分析：仅有肾脏、肝脏和脾脏中TS lncRNA的功能与对应的组织相关。在肝脏中，显著的GO条目有“甘油三酸酯稳态”“化学内稳态”“有机酸代谢过程”等（图3：A）。在脂肪中，显著的GO条目与mRNA剪接、mRNA处理等相关（图3：B）。

图3 通过顺式作用的TS lncRNA功能预测Fig. 3 Functional prediction of TS lncRNAs by cis-acting

为获取TS lncRNA潜在的反式作用位点，选取所有lncRNA和mRNA构建加权基因共表达分析。最终表达模式相近的转录本被分到相同的模块，每个模块被标记为单独的颜色，共获得42个模块（图4：A）。所有TS lncRNA都位于与其所在组织高度相关的模块（P＜0.05）。选取与每个TS lncRNA相关程度最高的10个mRNA，并对其进行功能富集分析。结果显示，每个组织中TS lncRNA的功能都与该组织相关，肌肉中TS lncRNA功能预测中显著的GO条目有“肌肉组织发育”“横纹肌细胞分化”等（图4：B），而脾脏中对应的GO条目有“淋巴细胞激活”“免疫系统过程的调节”等（图4：C）。

图4 通过反式作用的TS lncRNA功能预测Fig. 4 Functional prediction of TS lncRNAs by trans-action

2.4 调控IMF的lncRNA-mRNA筛选

从GEO数据库中下载有关IMF差异的家猪RNA-seq数据，通过差异表达分析，共获得239条显著差异表达的mRNA和245条lncRNA（图 5：A，B）。共筛选出 120个 lncRNA-mRNA作用对，功能富集分析发现23对lncRNA-mRNA与脂质代谢相关的通路有关，包括“PPAR信号通路”“脂肪酸代谢”“甘油酯代谢”等（图5：C）。其中，新鉴定的lncRNA——MSTRG.7256.5与脂质相关基因的连接最多，其中部分参与PPAR通路（图5：D）。因此，推测MSTRG.7256.5可能是调控IMF含量的因素。

MSTRG.7256.5位于第12条染色体，具有2个外显子，全长300 nt，其对应的最小自由能结构和质心二级结构如图所示（图5：E，F）。

图5 与IMF相关的lncRNA筛选Fig. 5 Identification of lncRNAs responsible for IMF

2.5 MSTRG.7256.5功能的验证

RT-qPCR的结果显示，构建的shRNA载体有效降低MSTRG.7256.5的水平后（图6：A），AQP7、SORBS和PPARγ的表达量增加，其中，PPARγ的增加显著，而SCD的表达量显著下降（图6：B）。随后构建并转染了MSTRG.7256.5的过表达载体，MSTRG.7256.5的表达水平升高后（图6：C），AQP7、SORBS和PPARγ的表达量显著下降，而SCD的表达量显著上升，整体上与干扰实验相反（图6：D）。

图6 MSTRG.7256.5作用的RT-qPCR验证Fig. 6 Correlation validation of MSTRG.7256.5 by RT-qPCR

油红O染色实验结果显示，MSTRG.7256.5干扰后的细胞内脂滴数量高于对照组（图7：A）。过表达了MSTRG.7256.5的细胞脂滴数量则更低（图7：B）。这些结果表明，MSTRG.7256.5可能通过PPAR途径调控家猪的IMF含量。

图7 MSTRG.7256.5干扰（A）与过表达后（B）的细胞内脂滴变化（油红O染色实验）Fig. 7 Change of lipid droplets after RNA interference（A） and overexpression （B）（oil red O staining experiment）

3 讨论

由于其营养价值和经济价值，家猪已经成为许多研究中的重要生物模型。lncRNA在许多生物过程中发挥着关键作用（Zhuet al.，2019），然而lncRNA-mRNA的相互作用，及其对组织特异性表达和IMF沉积影响的研究还有待深入。在本研究中，利用公共数据库中的RNA-seq数据，综合分析了lncRNAs的组织特异性表达模式以及对猪IMF含量的影响，为未来相关养殖业的研究提供了数据和参考。

在组装了猪的转录组后，共鉴定出23 678个lncRNA，其中14 309个为新鉴定的。经过特征分析，发现lncRNA与编码蛋白的mRNA相比，具有长度更短、表达更低、外显子更少等特征，这与之前的研究一致（Zouet al.，2018）。新鉴定的与已注释的lncRNA也存在明显的差异，在其他研究中也有类似的现象（Liet al.，2020），这可能是由目前lncRNA鉴定流程不够成熟造成的。有研究报道，lncRNA与蛋白编码基因相比具有更高的组织特异性（Cabiliet al.，2011；Kornienkoet al.，2016）。因此，本文通过计算TSI指数筛选了所有的TS lncRNA，其中，卵巢组织中的TS lncRNA数量明显高于其他组织。曾有研究报道lncRNA在哺乳动物卵巢和睾丸中特异性表达的数量高于其他组织（Iwakiriet al.，2017；Cavaet al.，2019），因此，推测lncRNA可能在生殖系统中有更高的组织特异性。

lncRNA可通过顺式或反式作用调控其靶基因（Gil & Ulitsky，2020），因此通常使用2种方式预测lncRNA的功能：一是根据其基因座附近的编码基因预测；二是根据与其共表达的基因预测。因此，本文通过寻找TS lncRNA的邻近基因（＜50 kb）和与其共表达的mRNA来预测lncRNA的功能。富集分析的结果显示，与顺式作用相比，lncRNA的反式作用靶点与其所在组织功能的相关性更高。卵巢TS ln-cRNA的邻近基因在功能富集分析中未获得显著结果，但与其共表达的mRNA在生殖系统相关的GO条目中显著富集。推测TS lncRNA倾向于通过反式作用发挥功能。

最后通过一系列流程筛选可能调控IMF含量的lncRNA-mRNA作用对，并构建了一个lncRNA-mRNA相互作用网络。发现1条新鉴定的lncRNA——MSTRG.7256.5与脂质代谢相关的PCG连接数最多，其中一些PCG富集在PPAR通路中，已有报道称，PPAR通路与脂质调节相关（Walczak & Tontonoz，2002；Li & Glass，2004）。 因此 ，本 文 对MSTRG.7256.5分别进行了RNA干扰和过表达实验，以验证其与潜在靶基因的关系。结果表明，MSTRG.7256.5表达趋势与AQP7、SORBS和PPARγ的表达呈负相关，与SCD的呈正相关，但SORBS的表达量在干扰实验中的变化不显著，AQP7的表达量在干扰实验与过表达实验中的变化均不显著。此前已有研究表明，AQP7和SORBS与脂质代谢相关（Yanget al.，2003；Skowronskiet al.，2007），而SCD可以调节猪的脂肪酸沉积和IMF含量（Yokotaet al.，2012；Ros-Freixedeset al.，2016）。此外，油红O染色实验的结果显示，MSTRG.7256.5的表达水平与细胞内的脂滴数量负相关，但过表达实验的结果不如RNA干扰实验的结果明显。

综上所述，认为MSTRG.7256.5可能通过PPAR途径负向调节家猪的脂肪沉积。然而，由于此步骤目的在于验证MSTRG.7256.5能否对其靶基因和脂肪沉积产生影响，因此仅选用了PK-15细胞用于实验，而PK-15细胞较低的脂肪含量可能导致实验结果不显著。总之，MSTRG.7256.5对其靶基因具体的作用机制和对肌肉间IMF含量的影响还须通过实验进一步探究。