静原鸡胸肌和腿肌肌苷酸特异性沉积相关circRNA的联合分析

王卫振 邓占钊 辛国省 虎红红 禹宝军 蔡正云 顾亚玲 张 娟*

(1.宁夏大学 农学院，银川 750021；2.彭阳县畜牧技术推广服务中心，宁夏 固原 756599；3.宁夏大学 生命科学学院/宁夏饲料工程技术研究中心，银川 750021)

随着人们生活水平的提高，膳食结构改善，如何提高肉品质已成为当前养禽业发展面临的重要课题[1]。研究表明，鸡肌肉中多种致鲜肽类和核苷酸是鲜味的主要来源，而肌苷酸(Inosine monphosphate, IMP)是其中最强的鲜味物质，并且与谷氨酸钠具有很强的协同作用，已成为评价肉质鲜味和新鲜程度的重要指标[2]。目前，已经有很多研究发现IMP的沉积受多种不同因素影响，如品种、饲养管理和营养因素等[3-5]。对于养殖业而言，利用遗传育种方法提高禽肉中IMP含量无疑是更经济高效的手段。

目前，对IMP沉积机制的研究已取得显著研究成果。对不同品种鸡的腺苷琥珀酸裂解酶(Adenylosuccinate lyase，ADSL)、甘氨酰胺核苷酸-5-氨基咪唑核苷酸合成酶(Glycinamide ridonucleotide synthetase-aminoimidazole ribonucleotide synthetase-glycinamide ribonucleotide formyltransferase，GARS-AIRS-GART)进行SNPs检测发现，ADSL基因TT型个体鸡的胸肌IMP含量极显著高于CC型个体，GARS-AIRS-GART基因TT型个体鸡的胸肌IMP含量极显著高于CC和CT型个体[6]。灵山鸡和大青麻鸡腺苷磷酸脱氨酶1(Adenosine monophosphate 1,AMPD1)基因第4、6、8外显子中分别检测出3个SNP位点，且SNP 6805A/G多态性与IMP含量呈显著正相关，可作为IMP沉积的候选基因[7]。Zhang等[8]在转录组测序中发现60日龄散养鸡AMPD1、胞质5′-核苷酸酶1A(Cytosolic 5-nucleotidase 1A,NT5C1A)和外核苷三磷酸二磷酸水解酶8(Ectonucleoside triphosphate diphosphohydrolase 8,ENTPD8)基因与120日龄散养鸡GART、GARS和ADSL基因表达量均显著高于同日龄的笼养鸡。

静原鸡作为宁夏地区优质地方品种肉质鲜美，研究其与肉质鲜味和新鲜程度紧密相关的IMP沉积机制是非常有必要的。课题组前期对静原鸡胸肌IMP研究发现，静原鸡公鸡和母鸡胸肌中IMP含量都显著高于腿肌；此外，在对其他品种公母鸡胸肌和腿肌IMP的测定中均发现胸肌中IMP的含量显著高于腿肌[10-12]；转录组测序发现腺苷酸激酶1(Adenylate kinase 1,AK1)基因在胸肌和腿肌中显著差异表达，并且AK1表达量与母鸡腿肌和胸肌中IMP含量呈正相关，与公鸡胸肌IMP含量呈负相关,与腿肌IMP含量呈正相关[9]。这些研究在基因水平揭示了IMP沉积的潜在调控机制，但基因的表达受多种因素的影响，尤其是非编码RNA的转录后调控和表观遗传修饰是影响特定基因表达量和翻译产物形成的重要因素。目前，对鸡肉中IMP沉积多集中于关键基因的筛选与功能验证，而对非编码RNA的联合分析尚未见报道。

本研究以RNA-seq技术对静原鸡(母鸡)胸肌和腿肌进行全转录组测序，从转录组水平上对IMP在胸肌和腿肌中特异性沉积机制进行系统的、全面的量化分析。通过对筛选出的静原鸡胸肌和腿肌中差异表达的circRNA进行生物信息学分析，进一步筛选出可能调控IMP特异性沉积的关键circRNA，并构建circRNA-miRNA-mRNA共表达互作网络，对IMP在胸肌和腿肌中特异性沉积的关键调控机制进行探究，为静原鸡优良地方鸡品种肉质性状调控机制研究提供参考。

1 材料与方法

1.1 试验材料

本研究所用的实验动物由宁夏彭阳县朝那鸡繁育中心提供，选择具有相同遗传背景同批孵化的静原鸡雏鸡，所用鸡在统一条件下饲养管理至180日龄，屠宰前12 h禁食,放血法屠宰，采集胸肌和腿肌肌肉组织切割成小块置于液氮中速冻，-80 ℃冰箱保存备用。

1.2 总RNA的提取与质控

用Trizol法提取胸肌和腿肌样品总RNA，总RNA用1%的凝胶电泳检测完整性，用Nanodrop ND-2000(Therom, 美国)和Agilent 2100 bioanalyzer(Agilent,美国)检测RNA的浓度及纯度，检测合格后-80 ℃保存备用。

1.3 cDNA文库的构建与转录组测序

RNA质检合格后，使用反转录试剂盒(TaKaBa, 大连)分别构建胸肌和腿肌共6个cDNA文库。从总RNA中去除rRNA、线性RNA，然后以片段化的RNA为模板合成cDNA，加A尾并连接测序接头，筛选350～400 bp的cDNA。对文库有效浓度准确定量，文库有效浓度≥2 nmol/L。库检合格后用Illumina HieqTM2 500平台进行双端测序。

1.4 测序数据处理和分析

对于原始测序数据通过FastQC软件进行质控，质控后的有效测序数据通过Hisat2和Bowtie2比对到鸡的参考基因组。获得的转录本使用String Tie进行组装拼接，通过FPKM算法确定基因表达定量，用EdgeR包筛选差异表达基因。利用Finder_circ和CIRI2识别、注释circRNA，用TPM对样本中circRNA表达量归一化处理，用DEseq2对样本中circRNA的表达量进行差异分析，差异表达circRNA的筛选从2个水平进行评估|log2Fold change|≥1，P<0.05。通过MiRBase、MiREvo和Mirdeep2数据库识别、鉴定、注释已知和未知miRNA。

1.5 差异表达基因的功能和通路富集分析

利用基因本体(Gene ontology，GO)数据库对差异表达基因(DEGs)进行细胞成分、生物学过程和分子功能进行分析。通过KEGG信号通路确定基因参与的主要生化代谢途径和信号转导途径。

1.6 circRNA和miRNA靶向位点预测

使用MiRanda和RNAhybrid数据库预测miRNA与mRNA间的靶向关系，使用MiRanda软件对鉴定的circRNA进行miRNA结合位点分析，MiRanda数据库是基于circRNA-miRNA结合的自由能进行靶向关系预测，为保证预测结果的可靠性，以得分值>145和自由能值<-12 Kcal /mol为筛选条件。

1.7 静原鸡转录组数据qRT-PCR验证

为确保转录组数据的准确性和可靠性，对6个差异表达circRNA和6个差异表达的miRNA进行qRT-PCR验证。荧光染料为TB Green Premix Ex TaqTMⅡ(TaKaBa, 大连)，使用Primer 5设计所挑选基因的引物，使用鸡actin-β(登录号：NM_205518.1)作为circRNA的内参，5S(登录号：NR_046276)作为miRNA的内参，所用引物均由陕西致研生物科技有限公司合成(表1)。

1.8 数据统计与分析

以 2-ΔΔCt平均值表达出各个基因的相对表达量，利用SPSS 23.0软件对circRNA和miRNA表达量进行单因素方差分析(P<0.05)。

2 结果与分析

2.1 静原鸡全转录组数据统计分析

本研究通过Illumina HieqTM2 500平台对静原鸡胸肌和腿肌进行深度测序，每个样本文库平均获得Raw reads约97 800 000，质控后平均获得clean reads 约975 00 000，每个样本平均获得约14.6G的有效数据量。测序数据Q20在98%以上，Q30在94%以上，测序质量良好，可用于进一步分析。测序数据结果显示每个样本GC碱基含量基本相等，碱基组成稳定均衡。比对结果显示约87% clean_reads比对到参考基因组，有79%的clean_reads唯一比对到参考基因组(表2)。为保证测序结果的准确性，排除异常样本带来的误差，计算评估各样本间pearson积矩相关系数R2均≥0.92。

表1 qRT-PCR引物序列Table 1 qRT-PCR primer sequence

表2 质控后基因比对分析Table 2 Analysis of gene alignment after quality control

2.2 circRNA鉴定与表达量分析

通过Finder_circ和CIRI2软件识别、鉴定转录本中的circRNA，共鉴定出circRNA 3 283个，在对转录本识别中，不同肌肉组织circRNA TPM值总体水平有差异(图1(a))。对比胸肌和腿肌circRNA的表达水平发现，仅在腿肌中表达的circRNA为80个，仅在胸肌中表达的circRNA为162个，共表达circRNA为3 041个，说明circRNA具有组织特异性。以P<0.05，|log2Fold change|≥1为筛选条件，共筛选出差异表达circRNA 446个，其中294个上调，152个下调(图1(b))。

蓝色表示上调，紫色表示下调，黑色表示没有差异。Blue indicates up-regulation, purple means down-regulation and black indicates no difference(a)不同样本之间circRNA TPM值分布情况;(b)circRNA火山图(a) Distribution of circRNA TPM values among different samples; (b) circRNA volcano plot图1 circRNA表达量分析Fig.1 circRNA expression analysis

2.3 差异表达基因和miRNA分析

对静原鸡胸肌和腿肌统计发现共鉴定出转录本30 252，其中已知转录本22 314，未知转录本7 938。利用FPKM算法标准化样本中转录本表达量，对转录本表达量的统计结果显示，FPKM>1 000的转录本有68，10

图2 差异表达mRNA (a)和miRNA (b)火山图Fig.2 Differential expression mRNA and miRNA volcano plot

表3 IMP合成相关基因表达量Table 3 Expression of genes related IMP synthesis

2.4 差异表达基因和差异表达circRNA的GO功能和KEGG信号通路富集分析

富集的GO条目结果显示，共有5 496个条目得到富集，有740个条目显著富集(P<0.05)，富集到生物过程497个条目、细胞组成82个条目、分子功能161个条目。2种肌肉组织差异表达circRNA显著富集于277个生物过程条目、67个细胞组成条目和81个分子功能条目。DEGs和差异表达circRNA显著富集条目结果显示，DEGs和差异表达circRNA都显著富集到碳水化合物代谢等生物过程，并且差异表达circRNA显著富集到基因的转录后调控等生物过程(表4)。KEGG通路富集结果显示，DEGs共富集到111个通路，有12条通路显著富集(图3(a))；差异表达circRNA共富集到87个通路，有6条通路显著富集(图3(b))。DEGs和差异表达circRNA都同时富集到糖酵解/糖异生通路(gga00010)、心肌细胞的肾上腺素信号传导(gga04261)，并且DEGs富集通路中氨基酸的生物合成(gga01230)、磷酸戊糖途径(gga00030)、丙氨酸、天冬氨酸和谷氨酸代谢(gga00250)等通路都与IMP的合成密切相关；差异表达circRNA富集通路中嘌呤代谢(gga00230)显著富集，此外还有多个信号通路被显著富集，如MAPK信号通路(gga04010)、FoxO信号通路(gga04068)。结合GO和KEGG分析可发现，静原鸡IMP的合成可能受到多种机制的调节，包括多种能量代谢通路、信号转导通路和氨基酸代谢通路。此外，在对差异表达circRNA的筛选中，筛选得到37个具有碳水化合代谢、IMP合成和嘌呤核苷酸合成等功能描述的circRNA，并根据其在胸肌和腿肌中的转录水平绘制了聚类分析热图(图4)，可以看出这些circRNA聚类结果良好，并且在样本中显著上调或下调。

表4 差异表达基因和circRNA最显著富集的GO条目Table 4 The most significantly enriched GO items of differentially expressed genes and circRNA

(a)差异表达基因KEGG通路富集分析，x轴表示富集度，y轴表示通路名称，气泡大小表示富集到通路中的基因数量，颜色表示富集到该通路富集显著性(b)差异表达circRNA的KEGG通路富集分析。(a) Enrichment analysis of KEGG pathway of differentially expressed genes, x axis represents enrichment degree, y axis represents path name, bubble size indicates the number of genes enriched in the pathway, color indicates enrichment significance (b) Enrichment analysis of KEGG pathway of differentially expressed circRNA.图3 差异表达基因和circRNA的 KEGG通路富集分析Fig.3 Enrichment analysis of differentially expressed genes and circRNA by KEGG pathway

2.5 circRNA和miRNA靶向关系分析

为分析circRNA、miRNA和mRNA之间的调控网络，挖掘出具有潜在价值的关键候选基因，通过miRanda和RNAhybrid数据库预测miRNA-mRNA靶向关系，共找到71 619对互作关系。KEGG通路富集分析发现，共有595对miRNA-mRNA互作关系与嘌呤代谢通路有关。对筛选出的miRNA-mRNA互作关系进行GO功能富集分析发现GMPS、PKM2、APRT、ITPA和IMPDH显著富集于前20个GO条目中(P<0.05)。使用miRanda数据库对circRNA-miRNA结合位点分析显示，共有90 603对circRNA-miRNA结合，筛选到446个差异表达circRNA均与多个miRNA互作，其中有234个circRNA-miRNA互作对与IMP合成代谢相关，构建circRNA-miRNA网络互作关系，可直观的反应非编码RNA之间的互作关系(图5)。circRNA-miRNA结合位点得分值越高、自由能越低表明靶向结合可能性越大，在234对circRNA-miRNA互作关系中依据得分值和自由能大小筛选到novel_circ_0001398、novel_circ_0013580、novel_circ_0012931、novel_circ_0011886和novel_circ_0013588 5个差异表达circRNA，表明这些circRNA可能在IMP的合成中发挥潜在的调控作用。

绿色表示下调，红色表示上调。Green indicates down-regulation, red indicates up-regulation.图4 差异表达circRNA的聚类分析热图Fig.4 Cluster analysis heat map of differentially expressed circRNA

黄色圆圈表示miRNA，三角形表示circRNA，绿色表示下调，红色表示上调；连接线表示两者之间互作关系。Yellow circle represents miRNA, triangle represents circRNA, green represents down-regulation, red represents up-regulation; the connecting line represents the interaction between two.图5 差异表达circRNA与miRNA的靶向关系Fig.5 Targeting relationship between differentially expressed circRNA and miRNA

2.6 构建circRNA-miRNA-mRNA共表达互作网络

对筛选到的circRNA、miRNA和mRNA构建了具有靶向关系的共表达互作网络。互作网络中6个差异表达circRNA与差异表达miRNA和差异表达mRNA共形成213对circRNA-miRNA-mRNA互作网络(图6)。结合KEGG和GO分析发现novel_409、gga-miR-456-3p、gga-miR-6660-3p、gga-miR-6631-5p、gga-miR-301b-3p、gga-miR-130b-3p、gga-miR-130c-3p、HIBCH、PKM2、PFKM、ESD和GART在共表达网络中发挥重要的调控作用。来源于亲本基因GPD2的novel_circ_0013580和novel_circ_0013588能够同时结合nove_409作用于靶基因GART，并且novel_circ_0013580还可以通过结合gga-miR-456-3p作用于靶基因HIBCH，来源于亲本基因ENSGALG00000032550的novel_circ_0012931能够分别结合gga-miR-6660-3p和gga-miR-6631-5p并作用于靶基因PKM2和HIBCH。对亲本基因功能分析发现，主要富集于内体运输、3-磷酸甘油代谢过程和糖异生等生物过程，基于以上结果推测circRNA通过调节物质运输和能量代谢影响IMP合成代谢。

紫色表示circRNA，蓝色表示miRNA，绿色表示mRNA，连接线表示存在互作关系。Purple indicates circRNA, blue indicates miRNA, green indicates mRNA, the connecting line indicates interaction.图6 circRNA-miRNA-mRNA网络互作分析图Fig.6 Analysis of circRNA-miRNA-mRNA network interaction

2.7 差异表达circRNA和miRNA的qRT-PCR检测

对6个circRNA和6个miRNA的qRT-PCR验证结果分析，所有circRNA和miRNA的溶解曲线均为单峰，引物设计合理。circRNA和miRNA胸肌和腿肌中相对表达量趋势与全转录组测序结果一致(图7)，表明测序结果可靠，可用于后续功能验证。对circRNA和miRNA的差异分析结果显示，novel_circ_0013580在胸肌和腿肌中显著差异表达、novel_circ_0013588、novel_circ_0012931、novel_circ_0011886和novel_circ_0008274在胸肌和腿肌中极显著差异表达；novel_409、gga-miR-6660-3p、gga-miR-1628、gga-miR-107-3p和gga-miR-6516-5p 在胸肌和腿肌中显著差异，gga-miR-196-1-3p在胸肌和腿肌极显著差异。

(a)circRNA qRT-PCR定量；(b)circRNA RNA-seq表达量;(c)miRNAqRT-PCR定量;(d)miRNA RNA-seq表达量；*P<0.05 表示差异显著，**P<0.01表示差异极显著;TPM表示每百万转录本，用于量化circRNA和miRNA的表达量(a) circRNA qRT-PCR quantification； (b) circRNA RNA-seq expression quantity； (c) miRNA qRT-PCR quantification； (d) miRNA RNA-seq expression quantity; *P<0.05 indicates a significant difference， **P<0.01 indicates that the difference is extremely significant; TPM represents transcripts per million， and used to quantify the expression of circRNA and miRNA图7 差异表达circRNA和miRNA的qRT-PCR和RNA-seq定量分析比较Fig.7 Quantitative analysis of comparison differentially expressed circRNA and miRNA by qRT-PCR and RNA-seq

3 讨 论

为探究胸肌中IMP的含量高于腿肌的分子机制，对静原鸡胸肌和腿肌进行全转录组测序，筛选调控IMP特异性沉积的关键circRNA、miRNA和mRNA，结合生物信息学对筛选到的基因进行功能分析和调控网络构建。本研究发现，筛选到的差异表达circRNA显著富集到IMP合成的嘌呤代谢通路与糖酵解/糖异生通路等，GO生物过程分析表明差异表达circRNA参与基因转录后调控和单一生物碳水化合物分解代谢过程等；对DEGs的分析中发现，DEGs显著富集到氨基酸的生物合成、磷酸戊糖途径、丙氨酸、天冬氨酸和谷氨酸代谢等通路。IMP有2条合成途径，从头合成途径和补救合成途径，在从头合成途径中以5-磷酸核糖焦磷酸、一碳单位、谷氨酸、天冬氨酸、CO2和ATP等原料在10种酶作用参与下合成IMP[13]。可见IMP的合成与上述通路中富集到的差异表达circRNA与DEGs具有显著的相关性。对胸肌和腿肌差异表达circRNA筛选发现共有17个差异表达circRNA与IMP的合成直接相关，并且这些差异表达circRNA呈现出组织特异性，有6个circRNA仅在胸肌中表达，有5个circRNA仅在腿肌中表达，这些特异性表达的circRNA可能与IMP的特异性沉积相关。结合GO和KEGG功能分析对转录本筛选，得到差异基因GART、NME2、POLR2H、PKM2、AMPD1、AMPD3和AK1等19个关键基因，这些基因涉及到嘌呤代谢、氨基酸生物合成、磷酸戊糖途径、糖酵解/糖异生等通路，这与Zhang等[8]对笼养和散养鸡的转录组学测序分析结果一致。

竞争性内源RNAs(ceRNAs)假说认为共享miRNA结合位点的转录本可以通过螯合影响miRNA的活性，解除miRNA对靶基因的抑制，从而调节miRNA靶标表达水平[14-15]。这些ceRNAs通过由miRNA识别元件(MRE)介导的一种新的“语言”进行通信，维持细胞内稳态。circRNA作为ceRNAs中的成员，能以海绵吸附的方式调节miRNA的活性和表达量[16]。在本研究中通过对circRNA-miRNA-mRNA的共表达分析，发现一些circRNA能够靶向miRNA调控基因的表达，如novel_circ_0006579、novel_circ_0013794和novel_circ_0006792下调而靶向的gga-miR-205c-5p上调则抑制POLR2H基因的表达，而novel_circ_0009612、novel_circ_0007952和novel_circ_0013239上调而靶向的gga-miR-490-3p下调则顺势调控PKM2基因的表达；这种miRNA与circRNA和mRNA表达水平相反的调控机制与前人的研究成果一致[17-18]。但值得注意的是，这种调控机制并非适用于所有的靶向关系，如novel_409上调靶基因GART同样上调，gga-miR-6660-3p上调PKM2基因也上调，并且gga-miR-490-3p和gga-miR-6660-3p能同时与novel_circ_0013239结合，但gga-miR-490-3p与novel_circ_0013239结合自由能低靶向关系更稳定，表明miRNA之间可能存在竞争同一circRNA或circRNA同一结合位点的关系，这可能是导致miRNA与circRNA或mRNA上下调关系相一致的原因。circRNA-miRNA-mRNA的共表达互作关系分析结果表明circRNA在调控IMP沉积中发挥着重要的调控机制。

PKM存在2种同工酶PKM1型和PKM2型是糖代谢的关键调控因子，糖酵解的限速酶[19-20]。PKM2转化催化磷酸烯醇式丙酮酸完成糖酵解途径生成丙酮酸和ATP[21]。PKM2是唯一受变构调控的基因，通过变构调节，PKM2重新编程细胞代谢，将糖酵解或转录过程的PKM2从四聚体切换成二聚体状态，四聚体PKM2定位于细胞质中起丙酮酸激酶的作用，产生ATP和丙酮酸，而二聚体PKM2定位于参与转录级联的细胞核中[22]。此外，发现PKM2敲除后，糖酵解和氧化磷酸化在PKM2敲除的细胞中减弱，糖代谢和脂代谢产物减少[23]。AMPD基因编码的AMP脱氨酶能够将AMP转化为IMP，在PKM2和AMPD3双缺陷型小鼠中提升了细胞内ATP、AMD和ADP的浓度[24]，而AMP的合成则是以IMP和天冬氨酸为底物。这从侧面验证了本研究成果，即胸肌中IMP的沉积量高于腿肌可能与胸肌中PKM2、AMPD1和AMPD3基因的上调相关。

GART是编码多功能酶基因家族的一部分位于鸡的1号染色体，与GARS和AIRS共同编码一个110 KDa三功能蛋白催化IMP从头合成第二、三、五步，在IMP的从头合成中至关重要[25-26]。GART具有四氢叶酸脱氢酶/环水解酶，NAD结合域，需要10-甲酰四氢叶酸作为辅助因子在IMP合成中具有显著特征，GARS和AIRS的结构域则不需要叶酸或其他辅助因子[27]。ADSL和ATIC基因的缺失导致严重嘌呤代谢紊乱，畸形和神经缺陷，GARS-AIRS-GART同是IMP合成的关键酶，在人类的唐氏综合征GARS-AIRS-GART的缺陷表达引起幼儿的先天愚型，嘌呤代谢水平显著提升[28]。仓鼠卵巢细胞中GART的突变能导致GARS和AIRS活性丧失而蛋白的稳定性和表达水平不变[29]。在对优质肉鸡IMP的研究中发现GARS-AIRS-GART的纯合基因型个体IMP含量显著高于其他非完全聚合基因型个体[30]。此外，研究还发现GARS-AIRS-GART基因第4外显子的突变对鲜味氨基酸和IMP的含量有显著的影响[31]。这些研究结果表明，在IMP的合成中GART扮演着重要的调控角色，本研究中GART在胸肌中的上调可能是IMP特异性沉积的重要原因。综上，可以推测PKM2和GART可能是调控IMP沉积的潜在标志性基因。

4 结 论

本研究基于全转录组学测序揭示了静原鸡胸肌和腿肌circRNA和mRNA表达谱，构建IMP特异性沉积的调控网络，通过对circRNA进行靶基因的预测和功能分析，得到6个circRNA与靶基因之间调控关系，筛选到的调控IMP特异性沉积核心调控因子novel_circ_0013580和novel_circ_0013588，明确了PKM2、GART2个靶基因与IMP的合成直接相关，本研究表明这些新发现的circRNA可能参与IMP合成中的调控网络，为地方鸡品种IMP的特异性沉积机制研究提供新思路。