用高通量测序技术研究松辽黑猪与长白猪背最长肌mRNA 和lnc...

奥鹏网院作业

用高通量测序技术研究松辽黑猪与长白猪背最长肌mRNA 和lncRNA 的差异表达
吴垚群1，陈少康3，盛熙晖1，齐晓龙1，王相国1，倪和民1，郭勇1，王楚端2，邢凯1

（1 北京农学院动物科学技术学院，北京 102206；2 中国农业大学动物科技学院，北京 100193；3 北京市畜牧总站，北京 100101）

摘要：【背景】肌肉生长和脂肪沉积是猪生产中非常重要的经济性状，同时也是复杂的数量性状，受到众多因素的影响。长白猪原产于丹麦， 是目前世界上使用最为广泛的瘦肉型猪种之一， 具有生长速度快， 瘦肉率高等特点，但肉质欠佳。松辽黑猪是以民猪、长白猪和杜洛克猪为素材经过三元杂交育成的黑色母系地方培育品种，具有繁殖率高、肉质优良、适应性强、耐粗饲等特性。二者在肉质方面存在较大差异。【目的】采用高通量测序技术对松辽黑猪和长白猪的背最长肌转录组进行分析，筛选影响猪肌肉生长、肉质和脂肪沉积的关键mRNA 和长链非编码RNA（long non-coding RNA， lncRNA），为猪肉品质研究提供新的参考信息。【方法】 采集6 头松辽黑猪和6头长白猪的背最长肌作为样本，提取其RNA，将每个品种内的个体取等量RNA 混池，采用Illumina HiSeq 2500高通量测序技术对混池后的样本RNA 进行测序，然后对所获得的Reads 进行比对、注释和差异表达等分析，筛选出差异表达mRNA 和lncRNA 并进行相关生物学功能富集分析。而后在筛选出的差异表达mRNA 和lncRNA 中各挑选10 个基因，包括5 个上调基因和5 个下调基因，采用荧光定量PCR（Real-time PCR， RT-PCR）方法鉴定所选mRNA 和lncRNA 的表达水平，并验证测序结果的可靠性。【结果】将得到的原始数据进行质量控制后，每个样本分别得到103 M 和150 M reads，其中98%以上reads 质量较高，满足进行下一步分析的要求。将reads比对到参考基因组上，共得到26 774 个转录本，其中包括8 428 个新的转录本。采用edgeR 软件包进行差异表达分析，参考标准是错误发现率（false discovery rate， FDR）小于0.05 和差异倍数大于2（fold change， FC）。在2 个猪种中共得到4 239 个显著差异表达的mRNA，其中在松辽黑猪肌肉中2 023 个上调，2 216 个下调。其中HLCS、BTD、DGKα、LPIN1、FGF1、FGF10、FGFR1 和ZNF7 等基因参与脂类代谢和肌肉发育相关的调控。差异表达的mRNA 的生物学功能富集分析发现，其主要富集在细胞发育、生物代谢调控、肌肉发育、脂类代谢等相关的信号通路和生物代谢途径上。同时，在2 个猪种中还检测到的178 个差异表达显著的lncRNAs，其中在松辽黑猪肌肉组织中有84 个上调，94 个下调。联合差异表达lncRNAs 的靶基因预测和差异表达基因结果分析，显示有4 个差异表达的lncRNA 与差异表达mRNA 存在潜在的调控关系，分别是TCONS-00041713、TCONS-00041712、ENSSSCT00000034982、ENSSSCT00000034269。同时10 个差异表达mRNA 和10 个差异表达lncRNA 的RT-PCR 结果显示：其表达水平在两组中同样差异显著，且与转录组测序中表达的趋势相似，表明高通量测序结果具有可靠性。【结论】 获得了影响猪肌肉生长、肉质和脂肪沉积的关键mRNA 和lncRNA，同时发现了几对差异表达lncRNA和mRNA 之间可能存在着潜在的调控关系。这些潜在调控关系的发现，可能对于研究调控肌肉组织代谢活动的分子机制的具有一定意义，并为其提供新的参考信息。

关键词：松辽黑猪；长白猪；背最长肌；mRNA；lncRNAs

0 引言
【研究意义】 长链非编码RNA（long non-coding RNA, lncRNA）是一类转录长度大于200nt 的非编码RNAs。目前许多研究已证实，lncRNA 参与了广泛的基因表达的调控过程，在细胞增殖、分化、凋亡以及肿瘤的发生和发展中发挥重要作用[1-2]。随着高通量测序技术的广泛应用，已经有众多的lncRNA 在肌肉发育和脂肪沉积中的作用被发现[3-5]。由此可见lncRNA是基因表达调控的重要因素,本研究的意义在于通过对不同品种猪背最长肌组织中lncRNA 的研究，进一步解析猪骨骼肌生长发育和脂肪沉积的分子机理。【前人研究进展】骨骼肌占哺乳动物总体重的50%左右[6]。多个研究发现lncRNAs 在骨骼肌发育过程中起着关键作用。在小鼠成肌细胞中，一个lncRNA 位于发育多能性相关2（developmental pluripotency associated 2, Dppa2）基因上游，其转录产物可以招募多种DNA 甲基化转移酶到Dppa2 基因的启动子区域，使Dppa2 基因沉默，进而促进成肌细胞的分化和肌肉再生[7]。lncRNA H19 由H19/胰岛素样生长因子2（insulin like growth factor 2, IGF2）基因座母体等位基因转录，在胚胎发育和成年机体的肌肉组织中表达较高，H19 通过结合某些复合物，起到抑制IGF2 转录的作用[8]。LU 等人发现 lncRNA YY1 相关肌发生 RNA 1（YY1-associated myogenesis RNA 1, Yam-1），受转录因子YY1（YY1 transcription factor, YY1）的正向调控，Yam-1 可以激活miR-715、抑制Wnt 家庭成员7B（Wnt family member 7B, Wnt7b）的表达，可以作为肌细胞抑制因子发挥作用[9]。在牛骨骼肌组织中，lncRNA-133a 可以促进牛骨骼肌卫星细胞的增殖，对牛骨骼肌卫星细胞的增殖有着正向调控作用[10]。lncRNA 作为基因表达的重要调控分子，可以通过多种作用机制，从转录水平、转录后水平以及染色质水平调节靶基因的表达，进而调控生命活动[11]。【本研究切入点】猪骨骼肌的生长发育与其生长速度和瘦肉率息息相关，直接关系着生猪养殖的生产效益。猪肌肉生长发育的差异是导致不同品种猪肉品质差异的重要原因。松辽黑猪是我国培育的第一个瘦肉型黑色母系品种，是由杜洛克猪、长白猪和东北民猪3 个品种经过长时间杂交培育而成的一个培育品种。松辽黑猪具有繁殖率高、肉质优良、适应性强、耐粗饲等特性，同时松辽黑猪的肌内脂肪含量可达3.19%，猪肉品质极佳[12]。长白猪原产于丹麦, 是目前世界上使用最为广泛的瘦肉型猪种之一，具有生长速度快，瘦肉率高等特点。研究发现，二者在肉质方面存在较大差异[12-13]。因此，以松辽黑猪和长白猪作为试验对象来研究lncRNA 与不同猪品种骨骼肌发育差异的关系，有利于研究lncRNA 在猪肌肉发育过程中的分子机制。目前已有多项研究证实lncRNA 在猪肌肉发育过程中发挥作用，但lncRNA 在中外不同品种猪肌肉生长发育过程中的差异性以及参与的生物学信息还有待于深入研究。【拟解决的关键问题】采用高通量测序技术对其背最长肌转录组进行测序分析，筛选差异表达的mRNA 和lncRNAs，并进行相关生物信息学的分析。最终筛选出关于猪肌肉发育和脂肪沉积的关键mRNA 和lncRNAs，以加深对肌肉生长发育和脂肪沉积分子调控机制的了解，为国内外猪品种在肉质和生长性状方面形成差异的分子机理研究提供新的参考信息。

1 材料与方法
1.1 试验材料
本研究所用的长白猪和松辽黑猪来自天津宁河原种猪场。所有动物均在统一条件下饲养管理，自由采食和饮水。饲养至体重为100kg 左右时进行屠宰。屠宰在华都阳光食品公司进行。屠宰标准按照GB/T 17236-2008 生猪屠宰操作规程进行[14]。屠宰后，立即将背最长肌组织移至液氮中进行冷冻保存，待提取RNA 使用。

1.2 总RNA 提取和建库
采用Trizol 法按照产品说明提取每个背最长肌组织的总RNA。用1%的琼脂糖凝胶电泳初步检测总 RNA 的完整性以及有无 DNA 污染。采用Nanodrop、Qubit 2.0 和Aglient 2100 方法检测各样品的浓度、纯度和完整性等。保证样品浓度≥500 ng·µL-1，28S:18S ＞1.0，RIN≥8。将每个品种的6 个样品的总RNA 等量混合后，进行混池建库，随后利用Illumina HiSeq 2500 (Illumina, America)平台进行双末端测序。文库的构建及测序在广州基迪奥生物科技有限公司进行。

1.3 数据的质控、比对和组装
FastQC 软件（http://www.bioinformatics.babraham. ac.uk/projects/fastqc/）用于对原始数据的质量控制，质控的标准为：（1）去除含adapter 的reads；（2）去除含N 比例大于10% 的reads；（3）去除低质量reads （质量值Q≤20的碱基数占整条read的50%以上）[15]。TopHat2 将高质量的序列数据（clean reads）比对到猪的参考基因组Sus scrofa 11.1 中[16]。Cufflinks 流程用于数据的转录本组装、注释与合并[17]。HT-seq 软件用于每个样本中基因表达量的计算[18]。

1.4 差异表达基因的功能富集分析
基因本体(gene ontology, GO )数据库 (http://www. geneontology.org/) 被用于差异表达基因的GO 功能富集分析，包括细胞组分，分子功能和生物学过程。超几何检验被用于P-value 的计算其公式为

pagenumber_ebook=193,pagenumber_book=839
式中：N 为所有Unigene 中具有GO 注释的转录本数目；n 为N 中差异表达转录本的数目；M 为所有Unigene 中注释为某特定GO 条目的转录本数目；m为注释为某特定GO 条目的差异表达转录本数目。FDR 被用于校正原始的P value。FDR＜0.05 的GO terms 被鉴定为显著富集的通路。京都基因和基因组百科全书(kyoto encyclopedia of genes and genomes, KEGG) 数据库 (https://www.genome.jp/kegg/) 被用于通路研究。同样，采用上述的方法和阈值鉴定差异表达基因的显著富集通路。

1.5 lncRNA 的鉴别
首先，筛选cuffcompare 转录本类别符号为“I”、“j”、“o”、“u”、“x”的转录本，其可能是新的lncRNA。进一步，为了剔除其中的小的非编码RNA，过滤掉其中外显子总长度小于200nt 及只有1个或无外显子的转录本。将剩下的转录本通过CPAT（Computerized Progressive Attention Training ）软件[19]进行lncRNA 的初步预测。为了得到可靠的lncRNA，继而利用CPC（Cetylpyridinium Chloride）软件[20]对CPAT 软件预测的结果进行进一步lncRNA 的识别。Cufflinks 用于对筛选出来的lncRNA 进行定量分析[17]。

1.6 基因及lncRNA 的差异表达分析
采用FPKM（fragments per kb per million fragments）的方法校正基因的表达水平，其同时考虑到了测序深度和基因长度对基因计数的影响[21]。edgeR 软件包被用于筛选差异表达基因[22]，其标准为FDR＜0.05 且|log2FC|＞1。

1.7 差异表达lncRNAs 的靶基因的预测
基于lncRNA 的顺式作用和反式作用方式预测差异表达lncRNA 的靶基因。在lncRNA 基因组位置上下游100kb 的基因认为是该lncRNA 潜在的靶基因。利用RNAplex 软件来预测lncRNA 通过碱基互补配对的反式作用模式的潜在靶基因[23]。

1.8 测序结果可靠性分析
为了验证测序结果的可靠性，对差异表达的mRNA 和lncRNA 进行验证。挑选10 个差异表达的mRNA 和10 个差异表达的lncRNA，包括5 个上调基因和5 个下调基因，进行RT-PCR 验证。RT-PCR所使用的cDNA 由高通量测序使用的剩余的RNA样品反转录产生。设计所挑选基因的引物，使用甘油醛-3-磷酸脱氢酶（glyceraldehyde-3-phosphate dehydrogenase, GAPDH）基因作为内部对照，引物均由上海生工合成（表1）。进行RT-PCR，并使用2-ΔΔCt方法计算样品之间基因的相对表达水平。

2 结果
2.1 测序数据的质量及比对信息
本研究通过Illumina HiSeqTM 2000 测序平台对松辽黑猪和长白猪的背最长肌组织mRNA 进行测序，每个样本平均得到约10G 数据量。测序数据的Q20在96%以上，Q30 在91%以上，测序质量良好。对原始数据质控后，平均每个样本得到18 224 475 795 条150bp 的clean reads（表2）。除去比对到核糖体RNA上的reads 后，平均81.22% 的clean reads 比对到猪的参考基因组中，其中唯一比对到基因组位置的为79.49%。各项数据显示测序数据的质量较好，可以进行后续的分析。

2.2 基因表达水平和差异分析
对转录本进行统计分析，松辽黑猪和长白猪得到的所有转录本数目分别是21 866 个和21 972 个，其中检测到新的转录本分别有7 846 个和7 838 个。对转录本的表达量统计发现，多数转录本的RPKM值在1 000以下，其中FPKM 值在1 以下转录本最多，为11 047个；FPKM 值1—10 之间的转录本为9 905 个；FPKM值在10—100 的转录本为2 915 个；而FPKM 值大于100 的转录本只有363 个，占总数的1.65%。同时，还发现两个样本的转录本表达丰度模式基本相同（图1）。

利用edgeR 软件包来筛选差异表达基因，显著差异表达的条件为 FDR＜0.05 且| log2FC|＞1 。结果显示两组之间一共有4 239 个基因显著差异，相对于长白猪，松辽黑猪组中2 023 个显著上调，2 216 个显著下调（图2）。

2.3 差异表达基因功能分析
为了解差异表达基因涉及的生物学功能，使用超几何检验的方法对差异表达基因进行GO 条目和KEGG 通路的功能富集分析。富集GO 条目的结果显示，差异表达基因主要分布在生物学过程、分子功能和细胞组成的三大类中的3 780 个GO 条目中，其中 1 314 个是显著富集的（Q-value≥0.05）。主要参与的 细胞膜的组成、跨膜运输、蛋白质结合、转运蛋白活性、能量代谢和氧化磷酸化等过程（表3）。对差异表达基因进行KEGG 通路功能富集分析，共鉴定出25个显著的通路，其中包括生物素代谢、泛酸盐辅酶A生物合成、多聚糖生成、磷脂酰肌醇信号系统、肌动蛋白调控、溶酶体等通路（图3）。

表1 引物序列
Table 1 Primer sequences

pagenumber_ebook=194,pagenumber_book=840
表2 质控前后样本测序数据的质量
Table 2 Quality of sample sequencing data before and after QC

pagenumber_ebook=194,pagenumber_book=840
SL 表示松辽黑猪，L 表示长白猪。SL stands for Songliao black pig, L stands for Landrace pig

2.4 新lncRNA 的鉴定及差异表达分析
通过cufflinks、CPAT 和CPC 流程鉴定转录本 中的lncRNA，共得到885 个新的lncRNA。所鉴定出的新lncRNA 中，在两个组中都有表达的lncRNA有792 个，只在松辽黑猪和长白猪中表达的基因分别为49 个和44 个 (图4-A)。同样利用edgeR分析软件筛选差异表达分析，结果显示松辽黑猪和长白猪之间一共有178 个非编码基因显著差异，其中在松辽黑猪中84 个显著上调，94 个显著下调（图4-B）。

pagenumber_ebook=195,pagenumber_book=841
图1 转录本表达量丰度分布函数图
Fig. 1 Transcript distribution abundance distribution function map

横坐标为log10（FPKM）值，数值越高，表示转录本表达量越高；纵坐标为转录本的密度，即为对应横轴表达量的转录本数/检测已表达转录本的总数；分布曲线的峰值表示整个样本转录本表达量最集中的区域
The abscissa is the log10 (FPKM) value. The higher the value, the higher the transcript expression; the ordinate is the density of the transcript, which is the number of transcripts corresponding to the horizontal axis expression/detection of the total number of expressed transcripts; The peak of the distribution curve represents the region where the expression of the entire sample transcript is most concentrated

pagenumber_ebook=195,pagenumber_book=841
图2 基因的差异表达分析
Fig. 2 Results of differential expression analysis of genes.

A 图中横坐标表示样品间的差异倍数对数值，纵坐标表示两个样品的-log10（FDR）值；红色表示长白猪相对于松辽黑猪的样本表达量上调和绿色表示表达量下调，黑色的点标志没有差异。B 图中横坐标表示成对比较的样品，纵坐标表示差异表达的基因数
In the graph A, the abscissa indicates the logarithm of the difference between the samples, the ordinate indicates the -log10 (FDR) value of the two samples, and the red indicates the up-regulation of the sample expression and the green expression of the Landrace pig relative to the Songliao black pig. Down, there is no difference in the black dot mark. In the B diagram, the abscissa indicates the pairwise comparison of the samples, and the ordinate indicates the number of genes differentially expressed

表3 差异表达基因最显著富集的GO 条目
Table 3 GO entries with the most significant enrichment of differentially expressed genes

pagenumber_ebook=196,pagenumber_book=842
2.5 lncRNA 的靶基因预测
基于lncRNA 的作用模式，通过cis 作用和反义作用预测差异表达lncRNA 的靶基因。根据lncRNA 的在染色体上的位置，共在241 个新预测到的lncRNA 的上下游500kb 内找到342 个编码基因。同时，在26 个已知的lncRNA 中发现81 个编码基因。使用RNAplex软件预测lncRNA 的反义作用靶基因，共鉴定得到84个lncRNA 与编码基因存在碱基的互补配对的关系。

结合两组间差异表达分析和靶基因预测的结果发 现，几对差异表达lncRNA 和基因之间可能存在着潜在的调控关系。如新lncRNA TCONS_00041713 和TCONS_00041712 在 ZNF7 基因的上游；lncRNA ENSSSCT00000034982在ERNBP基因的下游且与HCFC1存在碱基互补配对的关系；新lncRNA TCONS_00024150分别与FGFR1 和 ENSSSCT00000034063 基因间存在碱基互补配对的关系（图5）。以上的这些潜在调控关系为研究lncRNA 调控基因表达进而影响肉质品质提供新的素材。

pagenumber_ebook=197,pagenumber_book=843
图3 差异表达基因的KEGG 通路功能富集分析结果
Fig. 3 Functional enrichment analysis of KEGG pathways of differentially expressed genes

圆的大小和颜色分别代表富集通路和富集显着性的差异表达基因的量
The size and color of the circle represent the enrichment pathway and the amount of differentially expressed genes with significant enrichment, respectively

pagenumber_ebook=197,pagenumber_book=843
图4 lncRNA 的差异表达分析
Fig. 4 Differential expression analysis results of lncRNA

A 图表示鉴定出的新lncRNA，蓝色在松辽黑猪中表达的lncRNA，红色表示在长白猪中表达的lncRNA，中间的交集表示松辽黑猪和长白猪中均表达的lncRNA；B 图表示松辽黑猪和长白猪之间差异表达的lncRNA，横坐标表示样品间的差异倍数对数值，纵坐标表示两个样品的-log10（FDR）值；红色表示长白猪相对于松辽黑猪的样本表达量上调和绿色表示表达量下调
A figure shows the identified new LNCRNA, blue shows the LNCRNA expressed in Songliao black pig, red shows the LNCRNA expressed in landrace pig, and the middle intersection shows the LNCRNA expressed in both Songliao black pig and landrace pig; B figure shows lncRNA of difference expression between Songliao black pig and landrace pig, horizontal coordinate shows logarithm of difference multiple between samples, vertical coordinate shows-log10(FDR) value of two samples; The red color showed that the expression of landrace pig was up-regulated and the green color showed down-regulated compared with Songliao black pig

2.6 测序结果的验证
通过RT-PCR对随机挑选的10个差异表达mRNA和10 个差异表达lncRNA 进行验证。结果表明，随机挑选的10 个mRNA 的RT-PCR 结果与转录组测序结果表达趋势相似（图6-A），同样随机挑选的lncRNA的验证结果也与转录组测序结果表达趋势相似（图6-B），证实了松辽黑猪和长白猪背最长肌组织转录组测序结果具有可靠性，可以用于后续的深入研究。

pagenumber_ebook=198,pagenumber_book=844
图5 lncRNA 和基因之间存在的潜在调控关系
Fig. 5 Potential regulatory relationships between lncRNA and genes

不同颜色表示基因的表达水平，黄色表示基因在松辽黑猪中表达量高，蓝色表示基因在长白猪中的表达量高。箭头表示lncRNA 与基因之间存在的潜在调控关系。单向箭头代表基因上下游关系，箭头指向方向表示位于下游的基因；双向箭头表示基因之间存在碱基互补配对关系
Different colors indicate the level of gene expression. Yellow indicates that gene expression is high in Songliao Black pig, and blue indicates that gene expression is high in Landrace pigs. Arrows indicate potential regulatory relationships between lncRNA and genes. One-way arrows represent upstream and downstream relationships of genes, arrows point to downstream genes, and two-way arrows represent complementary pairing relationships between genes

pagenumber_ebook=198,pagenumber_book=844
图6 差异表达基因的RT-PCR 验证
Fig. 6 RT-PCR verification of differentially expressed genes

通过RT-PCR 验证松辽黑猪和长白猪背最长肌组织基因的差异表达
The differential expression of genes in the longissimus dorsi transcriptome of Songliao black pig and Landrace was verified by RT-PCR

3 讨论
本研究选用在肌肉品质和育肥期期间日增重方面存在较大差异的松辽黑猪和长白猪为对象，对其背最长肌进行RNA-seq 测序，筛选与肌肉生长和脂肪沉积相关的基因。结合生物学功能分析发现了一些与松辽黑猪和长白猪背最长肌肌肉发育和脂肪沉积可能相关的mRNA 和lncRNA。

猪的脂肪性状直接影响猪肉品质，如肌间脂肪含量、胴体瘦肉率等。猪脂肪组织的沉积能力具有品种差异，且不同品种的猪脂肪沉积能力具有较大差异[24]。本研究中，通过功能富集分析发现了一些与脂类代谢差异性显著的信号通路和代谢途径，如磷脂酰肌醇信号通路、鞘脂信号通路、甘油磷脂代谢途径、醛固酮的合成与分泌，与之相关的基因包括二酰基甘油激酶α(diacylglycerol kinase α, DGKα)、磷脂酰肌醇-4-激酶α(phosphatidylinositol 4-kinase α, PI4Kα)、磷脂酰肌醇5-磷酸 4 激酶 2α(phosphatidylinositol-5-phosphate 4-kinase type 2 α, PIP4K2α)、磷脂酶D1(phospholipase D1, PLD1)、脂肪酶F(lipase F, LIPE)、二酰基甘油脂酶α(diacylglycerol lipase α, DAGLα)、蛋白激酶D2(protein kinase D2, PRKD2)、脂质 1(lipin 1, LPIN1)、3-磷酸甘油酰基转移酶4(glycerol-3-phosphate acyltransferase 4, GPAT4)。在筛选得到的差异表达基因中，发现一些基因同时参与不同信号通路和代谢途径，如DGKα，其在磷脂酰肌醇信号通路和甘油磷脂代谢中均有参与。DGKα 属于二酰基甘油激酶家族，该酶类家族成员众多，主要作用为将1,2-二酰基甘油生成磷脂酸。DGKα 是一种在含有特殊脂类和蛋白质分子的外泌体生成的过程中发挥重要调控作用的因子，其在脂类代谢活动中具有重要作用[25]。前人研究发现，DGKα 作为LIPF 的下游分子，而LIPF 与食物中的甘油三酯在胃肠道中的消化过程相关[25]。因此可将DGKα 基因作为影响猪脂类代谢活动重要候选基因之一。HE 等研究发现在猪中LPIN1 基因与肌内脂肪沉积及肌肉品质存在显著相关性，可影响肌内脂肪含量、瘦肉率和风味，故将LPIN1 基因作为改善肌肉品质的候选基因之一[26]。这些基因现在已知的功能均与能量代谢和氧化磷酸化相关，但是在猪中影响组织脂类代谢的机制还有待深入研究。

育肥期期间日增重是猪的一种重要的经济性状，其受到肌肉发育直接影响。本研究发现了一些在松辽黑猪和长白猪之间显著差异的信号通路和代谢途径与肌肉发育有关，可能是影响育肥期期间日增重差异的因素。这些信号通路和代谢途径包括肌动蛋白调控过程、磷脂酰肌醇3-激酶/丝氨酸-苏氨酸蛋白激酶（phosphoinositide-3-kinase/serine-threonine kinase, PI3K- AKt）信号通路、单磷酸腺苷活化蛋白激酶（AMP- activated protein kinase, AMPK）信号通路，参与这些生物过程的基因包括成纤维因子1（fibroblast growth factor 1, FGF1）、成纤维因子10（fibroblast growth factor 10, FGF10）、成纤维细胞生长因子受体1（fibroblast growth factor receptor 1, FGFR1）、肌醇多磷酸-1-磷酸酶（inositol polyphosphate-1-phosphatase, INPP1）、磷脂酰肌醇-3-激酶调节亚基1（phosphoinositide-3-kinase regulatory subunit 1, PIK3R1）、血管生成素2（angiopoietin 2, ANGPT2）。研究已经证实，FGF1 在细胞增殖和细胞分化等多种生物学过程中发挥作用[27]；FGF10 在多种细胞的增殖、迁移、分化和生存中起到重要的支持作用[28]。PIK3R1 基因可编码磷脂酰肌醇3 激酶（phosphoinositide-3-kinase, PIK3）的调节亚基p85α蛋白[29-30]，而PIK3 作用于多种细胞外因子的信号转导，在细胞的生长、增殖、分化和凋亡等多种生命活动中起到重要作用[31]。可发现这些基因多数与细胞增殖、生长和分化等生命活动相关。在功能方面，多数与肿瘤细胞相关，但与肌肉细胞增殖生长的相关报道很少。因此，这些基因可作为研究肌肉发育的候选 基因进行深入研究和验证。

在生物学功能分析中，本研究发现生物素代谢过程最为显著。生物素，又称为维生素H、辅酶R，是一种在机体内生物活动过程中必不可少的辅酶因子[32]。GO 功能分析发现全羧化酶合成酶（holocarboxylase synthetase, HLCS）和生物素酶（biotinidase, BTD）基因参与了生物素代谢过程。HLCS 作用是能够催化生物素与羧化酶及组蛋白结合[33-35]。DO 等在对猪饲料利用效率的研究中发现该基因与剩余采食量性状具有显著关联性，发现全羧化酶合成酶活性丧失时会对机体组织的葡萄糖、脂肪及蛋白质的代谢产生重要影响[36]。而在本研究中发现HLCS 基因在两个品种猪背最长肌组织中表达显著差异，提示HLCS 可能参与了肌肉组织的生物代谢活动，与猪生长性状具有潜在联系。BTD缺乏时会导致肠道对生物素的摄取能力下降，导致体内生物素缺乏，使依赖于生物素的各种羧化酶（丙酰辅酶A 羧化酶、丙酮酰羧化酶、乙酰辅酶A 羧化酶和甲基巴豆酰辅酶A 羧化酶）的活性下降，影响糖类、脂肪及氨基酸代谢活动[37]。

在本研究中通过对差异表达的lncRNA 进行靶基因预测和功能分析，3 个差异表达lncRNA 与基因之间可能存在着潜在的调控关系，分别是 TCONS- 00024150、TCONS-00041713、TCONS-00041712。其中新lncRNA TCONS-00041713 和TCONS-00041712位于锌指蛋白7（zinc finger protein 7, ZNF7）基因的上游。ZNF7 属于锌指蛋白家族，前人研究发现，许多锌指蛋白家族成员在脂肪生成过程中起到关键的作用[38-39]。同时，在本研究发现新lncRNA TCONS- 00024150 与FGFR1 基因间存在碱基互补配对的关系，FGFR1 在脂肪组织中最先被发现，其可能在血管新生、糖脂调节、新陈代谢等方面发挥作用[40]。本研究表明新lncRNAs 可能参与了肌肉组织的脂类代谢活动，为解释猪品种间脂肪性状差异性提供新的思路，可为脂肪性状相关的基因研究奠定理论基础。

综上所述，通过对松辽黑猪和长白猪背最长肌之间差异表达基因的鉴定与分析，筛选出多个参与肌肉发育和脂类代谢的候选基因，其包括HLCS、BTD、DGKα、LPIN1、FGF1、FGF10、FGFR1 和ZNF7 基因。这些基因可能对研究松辽黑猪和长白猪两品种间肌肉发育和肌肉脂类代谢的差异性研究具有参考价值。

4 结论
本研究通过高通量测序技术对松辽黑猪和长白猪 背最长肌肌肉组织进行RNA 测序，通过基因比对分析，筛选出4 239 个编码基因差异mRNA，其中2 023个显著上调，2 216 个显著下调；与此同时，筛选出178 个差异表达的lncRNAs，84 个显著上调，94 个显著下调。通过对差异表达基因进行GO 功能分析和KEGG 通路分析共筛选出了HLCS、BTD、DGKα、LPIN1、FGF1、FGF10、FGFR1 和ZNF7 等与脂类代谢和肌肉发育相关的基因，并且参与了相关的信号通路和生物代谢途径。结合差异表达基因分析和靶基因预测的结果发现，几对差异表达lncRNA 和mRNA 之间可能存在着潜在的调控关系。如新 lncRNA TCONS_00041713 和TCONS_00041712 在ZNF7 基因的上游；lncRNA ENSSSCT00000034982 在ERNBP基因的下游且与HCFC1 存在碱基互补配对的关系；新 lncRNA TCONS_00024150 分别与 FGFR1 和 ENSSSCT00000034063 基因间存在碱基互补配对的关系。这些结果为了解松辽黑猪和长白猪之间肌内脂肪含量和肌肉发育性状差异的分子机制提供了新的参考。

References

[1] ARYAL B, ROTLLAN N, Fernandez-Hernando C. Noncoding RNAs and atherosclerosis. Current Atherosclerosis Reports. 2014, 16(5): 407.

[2] LEE J T, BARTOLOMEI M S. X-inactivation, imprinting, and long noncoding RNAs in health and disease. Cell, 2013, 152(6): 1308-1323.

[3] YANG Y, ZHOU R, ZHU S, LI X, LI H, YU H, LL K. Systematic identification and molecular characteristics of long noncoding RNAs in pig tissues. BioMed research International, 2017, doi: 10.1155/2017/ 6152582.

[4] XING K, WANG K, AO H, CHEN S, TAN Z, WANG Y, WANG C. Comparative adipose transcriptome analysis digs out genes related to fat deposition in two pig breeds. Scientific Reports, 2019, 9(1): 12925.

[5] 刘怡冰, 敖红, 邢凯, 王楚端. 脂肪相关lncRNAs 在动物中的研究进展. 中国畜牧兽医, 2019, 46(07):2021-2029. LIU Y B, AO H, XING K, WANG C D. Advances in research of fat-related lncRNAs in animals. Chinese Animal Husbandry Veterinarian, 2019, 46(07): 2021-2029. (in Chinese)

[6] ZHANG J, ZHOU C, MA J, CHEN L, JIANG A, ZHU L, SHUAI S, WANG J, LI M, LI X. Breed, sex and anatomical location-specific gene expression profiling of the porcine skeletal muscles. BMC Genetics, 2013, 14:53. DOI: 10.1186/1471-2156-14-53.

[7] WANG L J, ZHAO Y, BAO X C, ZHU X H, KWOK Y K Y, SUN K, CHEN X N, HUANG Y H, JAUCH R, ESTEBAN M A, SUN H, WANG H T. lncRNA Dum interacts with Dnmts to regulate Dppa2 expression during myogenic differentiation and muscle regeneration. Cell Research, 2015, 25(3): 335-350.

[8] DEY B K, PFEIFER K, DUTTA A. The H19 long noncoding RNA gives rise to microRNAs miR-675-3p and miR-675-5p to promote skeletal muscle differentiation and regeneration. Genes Development, 2014, 28(5): 491-501.

[9] LU L N, SUN K, CHEN X N, ZHAO Y, WANG L J, ZHOU L, SUN H, WANG H T. Genome-wide survey by ChIP-seq reveals YY1 regulation of lincRNAs in skeletal myogenesis. EMBO Journal, 2013, 32(19): 2575-2588.

[10] 李燕, 陈明明, 张俊星, 张林林, 李新, 郭宏, 丁向彬, 刘新峰. 牛lncRNA-133a 对骨骼肌卫星细胞增殖分化的影响. 中国农业科学, 2019, 52(1): 143-153. LI Y, CHEN M M, HANG J X, ZHANG L L, LI X, GUO H, DING X G, LIU X F. Effects of bovine lncRNA-133a on the proliferation and differentiation of skeletal muscle satellite cells. College of Animal Science and Veterinary Medicine, 2019, 52(1): 143-153. (in Chinese)

[11] LAU E. Non-coding RNA: Zooming in on lncRNA functions. Nature Reviews Genetics, 2014, 15(9):574.

[12] 张树敏, 陈群, 李娜, 金鑫, 赵晓东. 松辽黑猪种质特性及其利用的研究. 当代畜禽养殖业, 2006(06): 42-45. ZHANG S M, CHEN Q, LI N, JIN X, ZHAO X D. Study on germplasm characteristics and utilization of Songliao black pig. Contemporary Livestock and Poultry Farming, 2006(06): 42-45. (in Chinese)

[13] 贺希文, 高彩霞, 姜骞. 加系SPF 大白猪和长白猪群体遗传学分析. 中国实验动物学报, 2015, 23(6): 551-556. HE X W, GAO C X, JIANG Q. Genetic Analysis of the pedigreed population genetics of SPF Yorkshire and Landrace pigs imported from Canada. Chinese Journal of Experimental Animals, 2015, 23(6): 551-556. (in Chinese)

[14] 李江华, GB/T 17236-2008《生猪屠宰操作规程》的主要内容. 肉类研究, 2014(7):彩插. LI J H, GB/T 17236-2008 Main contents of the "Pest Management Procedure for Pigs". Meat Research, 2014(7). (in Chinese)

[15] WINGETT S, ANDREWS S. FastQ Screen: A tool for multi-genome mapping and quality control. F1000 Research, 2018, 7. DOI: 10. 12688/f1000research.15931.1

[16] KIM D, PERTEA G, TRAPNELL C, PIMENTEL H, KELLEY R, SALZBERG S L. TopHat2: accurate alignment of transcriptomes in the presence of insertions, deletions and gene fusions. Genome Biology, 2013, 14(4): R36.

[17] TRAPNELL C, ROBERTS A, GOFF L, PERTEA G, KIM D, KELLEY D R, PIMENTEL H, SALZBERG S L, RINN J L, PACHTER L. Differential gene and transcript expression analysis of RNA-seq experiments with TopHat and Cufflinks. Nature Protocols, 2012, 7(3): 562.

[18] ANDERS S, PYL P T, HUBER W. HTSeq-a Python framework to work with high-throughput sequencing data. Bioinformatics, 2015; 31(2): 166-169.

[19] SUN L, LUO H, BU D, ZHAO G, YU K, ZHANG C, LIU Y, CHEN R, ZHAO Y. Utilizing sequence intrinsic composition to classify protein- coding and long non-coding transcripts. Nucleic Acids Research, 2013, 41(17): e166.

[20] KONG L, ZHANG Y, YE Z Q, LIU X Q, ZHAO S Q, WEI L, GAO G. CPC: assess the protein-coding potential of transcripts using sequence features and support vector machine. Nucleic Acids Research, 2007, 35(Suppl.2): W345-W349.

[21] TRAPNELL C, WILLIAMS B A, PERTEA G, MORTAZAVI A, KWAN G, BAREN M J, SALZBERG S L, WOLD B J, PACHTER L. Transcript assembly and quantification by RNA-Seq reveals unannotated transcripts and isoform switching during cell differentiation. Nature Biotechnology, 2010, 28(5): 511.

[22] ROBINSON D M, MCCARTHY D J, SMYTH G K. edgeR: a Bioconductor package for differential expression analysis of digital gene expression data. Bioinformatics, 2010, 26(1): 139-140.

[23] TAFER H. RNAplex: a fast tool for RNA-RNA interaction search. Bioinformatics, 2008, 24(22): 2657-2663.

[24] WOOD J D, ENSER M, FISHER A V, NUTE G R, SHEARD P R, RICHARDSON R I, HUGHES S I, WHITTINGTON F M. Fat deposition, fatty acid composition and meat quality: A review. Meat Science, 2008, 78(4): 343-358.

[25] ALONSO R, MZAAEO C, RODRIGUEZ M C, MARSH M, FRAILE R A, AVILAFLORSE A, MERIDA I, IZQUIERDO M. Diacylglycerol kinase alpha regulates the formation and polarisation of mature multivesicular bodies involved in the secretion of Fas ligand-containing exosomes in T lymphocytes. Cell death and Differentiation, 2011, 18(7): 1161-1173.

[26] HE X P, XU X W, ZHAO S H, FAN B, YU M, ZHU M J, LI C C, PENG Z Z, LIU B. Investigation of LPIN1 as a candidate gene for fat deposition in pigs. Molecular Biology Reports, 2009, 36(5): 1175-1180.

[27] GROSE R, DICKSOM C. Fibroblast growth factor signaling in tumorigenesis. Cytokine & Growth Factor Reviews, 2005, 16(2): 179-186.

[28] JOANNE W, CHIARA F. Regulation of FGF10 signaling in development and disease. Front Genetics, 2018, 9:500.

[29] KATSO R, OKKENHAUG K, AHMADI K, WHITE S , TIMMS J, WATERFIELD M D. Cellular function of phosphoinositide 3- kinases: implications for development, homeostasis, and cancer. Annual Review of Cell and Developmental Biology, 2001, 17: 615-675.

[30] OKKENHAUG K,VANHAESEBROECK B. New responsibilities for the PI3K regulatory subunit p85 alpha. Science's STKE, 2001(65):pe1. DOI: 10.1126/stke.2001.65.pe1.

[31] CARNERO A, BLANCO-APARICIO C, RENNER O, LINK W, LEAL J F M. The PTEN/PI3K/AKT signaling pathway in cancer, therapeutic implications. Current Cancer Drug Targets, 2008, 8(3): 187-198.

[32] 郑巍振. 生物素高产菌种选育的研究[D]. 杭州: 浙江工业大学, 2010． ZHENG W Z. Study on strain screening of high biotin-producing[D]. Hangzhou: Zhejiang University of Technology, 2010. (in Chinese)

[33] WOLF B. Disorders of biotin metabolism [M]//SCRIVER C R, BEAUDET A L, SLY W S. The Metabolic and Molecular Basis of Inherited Disease, New York: McGraw-Hill, 2011: 3935-3956.

[34] STANLEY J S, GRIFFIN J B, ZEMPLEMI J. Biotinylation of histones in human cells: Effects of cell proliferation. FEBS Journal, 2001, 268(20): 5424-5429.

[35] KOBZA K, CAMPOREALE G, RUECKERT B, KUEH A, GRIFFIN J B, SARATH G, ZEMPLENI J. K4, K9 and K18 in human histone H3 are targets for biotinylation by biotinidase. FEBS Journal, 2005, 272(16): 4249-4259.

[36] DO D N, OSTERSEN T, STRATHE A B, MARK T, JENSEN J, KADARMIEDEEN H N. Genome-wide association and systems genetic analyses of residual feed intake, daily feed consumption, backfat and weight gain in pigs. BMC Genetics. 2014, 15:27.

[37] WEI S, ZHANG L, ZHOU X, DU M, IIANG Z, HAUSMAN G J, BERGEN W G, ZAN L, DODSON M V. Emerging roles of zinc finger proteins in regulating adipogenesis. Cellular and Molecular Life Sciences, 2013, 70, 4569-4584.

[38] TSENG K,Y CHEN, Y H, LIN S. Zinc finger protein ZFP36L1 promotes osteoblastic differentiation but represses adipogenic differentiation of mouse multipotent cells. Oncotarget, 2017, 8, 20588-20601.

[39] HU X M, ZHOU Y F, YANG Y, PENG J, SONG T X, XU T, WEI H K, JIANG S W, PENG J. Identification of zinc finger protein Bcl6 as a novel regulator of early adipose commitment. Open Biology, 2016, 6, 1-14.

[40] SPOLCOVA A, HOLUBOVA M, MIKULASKOVA B, NAGELOVA V, STOFKOVA A, LACINOVA Z, JURCOVICOVA J, HALUZIK M, MALETINSKA L, ZELEZNA B. Changes in FGF21 serum concentrations and liver mRNA expression in an experimental model of complete lipodystrophy and insulin-resistant diabetes. Physiology Research, 2014, 63(4): 483-490.

Differential Expression of mRNA and lncRNA in Longissimus Dorsi Muscle of Songliao Black Pig and Landrace Pig Based on High-Throughput Sequencing Technique

WU YaoQun1,CHEN ShaoKang3,SHENG XiHui1,QI XiaoLong1,WANG XiangGuo1,NI HeMin1, GUO Yong1,WANG ChuDuan2,XING Kai1
(1College of Animal Science and Technology, Beijing University of Agriculture, Beijing 102206; 2College of Animal Science and Technology, China Agricultural University, Beijing 100193; 3Beijing Animal Husbandry Station, Beijing 100101)

Abstract:【Background】Muscle growth and fat deposition are important economic traits in pig production and are also complex quantitative traits. Landrace, originally from Denmark, is one of the world-famous lean-type, with fast growth and high lean meat rate, but its meat quality is poor. Songliao black pig is a black maternal local cultivar cultivated by ternary cross-breeding with the Min pig, the Landrace pig and the Duroc pig, which has the characteristics of high reproduction rate, good meat quality, strong adaptability, roughage resistant and etc. There is a big difference in meat quality between the Songliao Black Pig and Landrace Pig. 【Objective】 In this study, high-throughput sequencing technique was used to analyze the longissimus dorsi transcriptome of Songliao black pig and Landrace. The key mRNA and long non-coding RNA（lncRNA）, affecting muscle growth, meat quality and fat deposition were selected to provide new reference information for the study of pork quality. 【Method】We collected the longissimus dorsi muscle of 6 Songliao black pigs and 6 Landrace to extract RNA. Individuals in each cultivar were sampled in the same RNA mixing pool and sequenced by Illumina HiSeq 2500 high-throughput sequencing technology. Reads were compared, annotated and differentially expressed. Differentially expressed genes and lncRNAs were screened out and their biological functions were enriched. Moreover then, ten genes, including five up-regulated genes and five down-regulated genes, were selected from the differentially expressed RNA and lncRNA, and the relative expression levels of differentially expressed mRNA and lncRNA were detected by real-time PCR. 【Result】After quality control, each sample received 103 M and 150 M reads, of which more than 98% of the reads were of higher quality, meeting the requirements of further analysis. The reads were aligned to the reference genome and a total of 26 774 transcripts were obtained, of which 8 428 were new transcripts. Using edgeR for differential expression analysis, the reference standards are False Discovery Rate(FDR)＜0.05 and fold change (FC)＞2. A total of 4 239 significant differentially expressed genes were obtained, of which 2 023 were up-regulated and 2 216 were down-regulated in Songliao black pig muscle. Meanwhile, 178 significant differentially expressed lncRNAs were detected, of which 84 were up-regulated and 94 were down-regulated in Songliao black pig muscle tissue. The biological function enrichment analysis of differentially expressed genes revealed that it is mainly enriched in cell components and cell development, biological metabolism regulation, and signaling pathways and biological metabolic pathways involved in muscle development and lipid metabolism. The target gene prediction results of the screened lncRNA showed that there was a potential differential relationship between the four differentially expressed lncRNAs and the differentially expressed genes, namely TCONS-00041713, TCONS-00041712, ENSSSCT00000034982, ENSSSCT00000034269. At the same time, the results of real-time PCR analysis of 10 differentially expressed mRNAs and lncRNAs showed that the results of real-time PCR were similar to the expression of transcriptome sequencing of the longissimus dorsi muscle tissue of Songliao black pig and Landrace, indicating that high throughput sequencing results are reliable.【Conclusion】We obtained the key mRNA and lncRNA that affect the muscle growth, meat quality, and fat deposition of pigs. Meanwhile, it was found that there may be a potential regulatory relationship between several pairs of differentially expressed lncRNA and mRNA.The discovery of these differentially expressed mRNAs and lncRNAs and possible potential regulatory relationships may have a certain significance for studying the molecular mechanisms that regulate metabolic activity in muscle tissue and provide new reference information.

Key words: Songliao black pig; Landrace pig; longissimus dorsi muscles; mRNA; lncRNAs

收稿日期：2018-12-10；

接受日期：2019-12-05

基金项目：国家重点研发计划资助（2018YFD0501000）；现代农业产业技术体系北京市创新团队(BAIC02-2018)；北京市教委2019 年度科技计划一般项目（KM201910020010）；北京农学院青年教师科研基金（SXQN201807）；北京农学院“大北农青年教师科研基金”

联系方式：吴垚群，E-mail：1530321769@qq.com。

通信作者王楚端，E-mail：cdwang@cau.edu.cn。

通信作者邢凯，E-mail：xk181986@163.com

开放科学（资源服务）标识码（OSID）：

pagenumber_ebook=190,pagenumber_book=836
（责任编辑 林鉴非）