靜原雞胸肌和腿肌肌苷酸特異性沉積相關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的特異性沉積機制研究提供新思路。