基于轉錄組測序鑒定牛RFI表型差異相關circRNAs及其靶基因富集通路

王鵬飛，丁燕玲，楊朝云，馬 應，周小南，史遠剛，康曉龍

（寧夏大學農學院，銀川 750021）

飼料成本在家畜養殖成本中占比最高，提高飼料利用率、降低養殖成本為家畜規?；B殖及提高效益重要環節。剩余采食量（Residual feed intake，RFI）為家畜實際采食量與其維持生產（產蛋、產奶等）的預期采食量之差[1]，為評價動物飼料轉化率及代謝水平指標[2]，本試驗開展家畜RFI表型研究，為高飼料效率表型的家畜選育提供理論基礎。

環狀RNA（circRNAs）為一種非編碼RNA，由線性RNA反向剪接產生[3]，與其他非編碼RNA相比，circRNAs具有封閉圓形結構，更穩定，不易降解[4-5]。circRNAs由內含子和外顯子構成，可調節多種生物學功能，如疾病、肌肉發育、脂肪代謝等[6]。circRNAs最重要功能為充當miRNA的“海綿”，消除miRNAs對靶基因抑制[7]，如在沼澤水牛乳腺上皮細胞中，circ_015003和circ_014121的海綿miR-103過表達增加乳脂合成相關基因表達，導致脂滴形成及甘油三酯累積[8]；ssc_circ_0002807和ssc_circ_0009352為miRNA-433的海綿，miRNA-433能激活p38MAPK通路，而激活的p38MAPK促進骨髓間充質干細胞脂肪生成[9]。研究表明，miRNA參與調控動物機體采食及能量代謝，如miRNA-499和miRNA-15a可分別通過其靶基因TCF12和FOXO1介導AMPK-PGC-1α、TGF-β及PI3KAKT-mTOR等信號通路增強機體對營養物質吸收利用，調控動物采食[10]。circRNAs可直接參與該生物學過程，或間接通過與其他分子（如miRNAs等）互作共同調控目標性狀，表明circRNAs在調控表型變異過程中具有廣泛且獨特的生物學功能。

目前，circRNAs在畜禽主要表型性狀研究集中于脂肪沉積[11]、肌肉發育[12]、泌乳性能[13]等，但circRNAs是否參與調控牛RFI表型變異，或在牛營養物質吸收和能量代謝過程中是否具有潛在生物學作用鮮有報道。下丘腦為調節食欲和能量平衡關鍵中樞，通過調節特定神經遞質/神經肽表達對營養物質/激素信號作出反應，調控能量攝入和消耗[14]。弧形核（Arcuate nucleus，ARC）被認為是下丘腦控制食物攝入量主要中心，其兩個不同神經元群體接收外周營養信號，通過不同食欲肽表達參與攝食調節，如促食類神經肽有神經肽Y（Neuropeptide Y，NPY）、刺鼠相關蛋白（Agoutirelated protein，AGRP）等。研究發現，活化蛋白激酶（AMP-activated protein kinase，AMPK）在下丘腦神經元中扮演重要角色，作為身體能量穩態的傳感器和調節器，將代謝狀態與神經肽系統聯系起來，調節攝食[15]。鑒于上述背景，本研究采用高通量測序技術鑒定牛下丘腦中RFI表型相關差異circRNAs，并對靶基因富集分析，旨在揭示circRNAs在牛日糧攝入、能量代謝和消化吸收中潛在作用，為后續全面探究circRNAs在動物主要經濟性狀表型變異中分子功能及調控作用提供參考依據。

1 材料與方法

1.1 研究對象

選取健康狀況良好初始體重相近（280±30）kg，月齡相近（16±1）mo高RFI（High RFI，HRFI）、低RFI（Low RFI，LRFI）秦川牛各5頭。試驗牛頸部放血屠宰開顱后，立即從下丘腦前段剪取1 cm×1 cm下丘腦組織，生理鹽水沖洗裝入凍存管，立即放入液氮保存。HRFI、LRFI牛均設5個重復，分別編 號 為HRFI_1、HRFI_2、HRFI_3、HRFI_4、HRFI_5和LRFI_1、LRFI_2、LRFI_3、LRFI_4、LRFI_5。依據標準方法計算RFI[16]。試驗牛飼養方式為單欄飼喂，每日兩次（8:00和15:00），采食1.5 h，自由飲水。

1.2 RNA提取及文庫構建

總RNA提取采用Trizol法。對提取的RNA作濃度、純度（OD230，OD260和OD280）及完整性檢測，評估總RNA質量。總RNA的OD260/OD280值≥1.8，完整性值在8.6以上，表明總RNA質檢合格，用于文庫構建。合格RNA樣品使用Covaris超聲波破碎儀隨機打斷成長度約350 bp片段，經末端修復、加A尾、加測序接頭、純化、PCR擴增等步驟完成文庫制備。使用Qubit2.0作初步定量，稀釋文庫至2 ng.μL-1，隨后使用Agilent 2100檢測文庫insert size，符合預期后，使用qRT-PCR方法準確定量文庫有效濃度（文庫有效濃度高于3 ng.μL-1），保證文庫質量。文庫質檢合格后，用Illumina平臺作雙端測序，得到用于分析的原始測序片段（raw reads）。

1.3 差異表達circRNAs篩選

為確保數據分析可靠性，使用fastx_tool-kit軟件（v0.0.14）去除raw reads中帶接頭、低質量序列及接頭污染等，得到clean reads。采用BWA中BWA-MEM算法將clean reads比對到參考基因組中，采用CIRI識別circRNAs，比對結果以SAM文件輸出后掃描該文件，找到PCC（Paired chiastic clipping）和PEM（Paired-endmapping）位點及剪接信號GT-AG，最后為確保識別circRNAs可靠性，將具有junction位點的序列重新比對。用TPM（Transcripts per million）對circRNAs表達量作歸一化處理，DEGseq分析樣本間circRNAs表達差異，以|Fold Change|≥1和P＜0.05作為表達標準確定差異表達的circRNAs。

1.4 差異表達circRNAs功能富集分析

差異表達circRNAs的GO及KEGG富集分析采用DAVID數據庫；GO富集分析包括3個層次，分為生物學過程（Biological process，BP）、細胞組分（Cellular component，CC）、分子功能（Molecular function，MF）。KEGG數據庫是將差異表達circRNAs所對應基因作通路富集，確定來源基因主要參與的代謝途徑和信號通路。

1.5 差異表達circRNAs驗證

為驗證測序結果準確性，以GAPDH為內參基因，隨機選取6個差異表達circRNAs（novel_circ_0031230、novel_circ_0038739、novel_circ_0022168、novel_circ_0005798、novel_circ_0018798、novel_circ_0034790）作qRT-PCR驗證，引物由上海生工合成（見表1）。取RNA樣本反轉錄為cDNA，以此為模板作qRT-PCR檢測，具體操作方法依照TaKaRa試劑盒說明書。qRT-PCR反應程序為95℃預變性30 sec，95℃變性10 sec，退火30 sec，4℃保存，共40個循環。采用公式2-ΔΔCt計算差異circRNAs相對表達量。

表1 基因及對應引物Table 1 Gene and corresponding primer sequences

2 結果與分析

2.1 下丘腦組織circRNAs測序數據質量分析

去除測序數據接頭序列、低質量及錯配reads等，得到HRFI、LRFI組中clean reads（見表2）。由表2可知在所測10個個體中clean reads占總raw reads比例均在97%以上，表示測序數據質量良好，可用于后續分析。

表2 數據過濾統計情況Table 2 Data filtering statistics

2.2 circRNA來源及長度統計

circRNAs可分為exon、intergenic、intron 3種類型，通過與參考基因組比對，發現circRNAs在exon中占比最高（約87%），intron中占比最少（約5%），結果見圖1A。證實circRNAs為閉合環狀RNA，大部分由外顯子構成，故無多聚A尾巴；長度統計表明大部分circRNAs分布在100～500 nt（見圖1B）。

圖1 circRNAs來源（A）及長度（B）Fig.1 Source(A)and length(B)of circRNAs

2.3 差異circRNAs篩選

采用TPM對各樣本中已知和新circRNAs表達量作歸一化處理（見圖2A），通過差異倍數（|log2（Fold change）|＞1）和顯著水平（P＜0.05）兩個水平篩選兩組中差異表達circRNAs，差異基因分布情況以火山圖表示（見圖2B），從HRFI組中篩選得到736個與LRFI差異表達的circRNAs，其中表達上調有315個，下調有421個。差異顯著上調和下調前10個circRNAs相關信息見表3，其中novel_circ_0031230與bta-miR-103具有靶向關系。據報道，bta-miR-103及其靶基因（Dicer、PANK）可調節脂肪生成[17]等。

表3 顯著差異前10的circRNAsTable3 Top 10 significant differential expressed circRNAs

2.4 差異表達circRNAs靶基因GO富集分析

對736個差異表達circRNAs靶基因作GO分類（見圖3）。

由圖3可知，在BP中，差異表達circRNAs對應基因主要涉及細胞器組織、細胞定位、蛋白質定位及大分子定位等；在MF中主要涉及GTP酶結合、Ras GTP酶結合、酸性胺酸連接酶活性等過程；在CC中主要富集在細胞器、細胞內膜、細胞質及細胞核組成等；其中細胞內細胞器、膜結合細胞器和細胞內膜細胞器注釋高，說明與細胞器生成密切相關，而線粒體調控能量代謝，內質網促使蛋白質加工及脂質合成，核糖體將氨基酸合成蛋白質，最終通過高爾基體加工合成多糖[18]，這對機體能量代謝發揮重要作用。以上結果表明，差異表達circRNAs可在不同RFI表型牛脂肪合成和能量代謝中起重要作用。

圖3 高、低RFI牛下丘腦差異circRNAs的GO功能富集分析Fig.3 GO enrichment analysis of hypothalamusdifferential circRNAs in bovine with high and low RFI

2.5 差異表達circRNAs靶基因的KEGG富集分析

使用KEGG分析探討差異表達circRNAs生物學作用。KEGG富集分析表明（見圖4），差異表達circRNAs靶基因主要富集在AMPK信號通路、mTOR信號通路及肌動蛋白細胞骨架形成通路中。AMPK信號通路在骨骼肌和肝臟中發揮重要作用，如丁酸通過增加肌肉和肝臟中AMP/ATP比率，直接激活AMPK信號通路[19]，該通路激活后可抑制糖原和蛋白質合成，葡萄糖6-磷酸酶和磷酸烯醇式丙酮酸羧激酶基因表達下降，最終葡萄糖轉運和脂肪酸氧化增加[20]。這些信號通路的富集為后續開展circRNAs在牛食物攝入及能量代謝過程中分子功能提供重要信息。

圖4 高、低RFI牛下丘腦差異circRNAs的KEGG功能富集分析Fig.4 KEGG enrichment analysis of hypothalamus differential circRNAs in bovine with high and low RFI

2.6 差異表達circRNAs驗證

為驗證測序數據可靠性，對4個差異表達circRNAs（novel_circ_0031230、novel_circ_0038739、novel_circ_0022168、novel_circ_0005798、novel_circ_0018798、novel_circ_0034790）作qRT-PCR分析。結果表明（見圖5），在高RFI牛下丘腦組織中novel_circ_0031230、novel_circ_0038739、novel_circ_0018798表達上調，novel_circ_0022168、novel_circ_0005798、novel_circ_0034790表達下調，定量結果與測序結果趨勢一致，說明測序數據相對可靠，相關差異分子可用于后續功能驗證分析。

圖5 高、低RFI牛下丘腦6個差異表達circRNAs的qRT-PCR驗證Fig.5 qRT-PCR verification of six differentially expressed circRNAs in hypothalamus of brovine with high and low RFI

2.7 miRNA靶點預測及circRNA-miRNA-mRNA網絡分析

circRNAs富含miRNA結合位點，通過miRNA海綿效應吸附miRNA，解除其對靶基因抑制作用。HRFI組315個上調和421個下調的差異表達circRNAs中共篩選到1 018個miRNAs靶點，其中在HRFI牛下丘腦差異表達上調最大的novel_circ_0037052靶向結合34個miRNAs，差異表達下調最大的novel_circ_0000923結合31個miRNAs靶點，提示上述miRNAs可能在牛RFI調控方面具有重要作用。通過篩選的miRNAs靶點制作circRNA-miRNA-mRNA靶向關系圖，圖6展示顯著差異表達的4個circRNAs及其靶向miRNAs，mRNA。結果顯示，差異circRNAs靶向結合的bta-miR-103、btamiR-33 a/b、bta-miR-499、bta-miR-15a等miRNAs調節動物脂肪沉積及能量代謝等相關mRNA。

圖6 差異表達circRNA-miRNA-mRNA調控網絡Fig.6 Differentially expressed of circRNA-miRNA-mRNA regulatory network

3 討論

本研究分析高、低RFI牛下丘腦組織中circRNAs，共得到736個差異表達的circRNAs，其中315個在HRFI牛中表達上調，421個表達下調。通過qRT-PCR驗證在HRFI牛下丘腦組織中novel_circ_0031230和novel_circ_0038739顯著上調；novel_circ_0022168和novel_circ_0005798顯著下調，與測序結果趨勢一致，說明本次測序數據真實可靠。對差異circRNAs的來源基因作GO和KEGG富集分析發現，差異circRNAs來源基因主要參與細胞質及細胞器組成，能量代謝相關的AMPK信號通路、mTORx信號通路等生物學過程。

本研究預測circRNAs-miRNAs靶標關系，得出bta-miR-103、bta-miR-33 a/b、bta-miR-499及bta-miR-15a分別為novel_circ_0031230、novel_circ_0020172、novel_circ_0002558、novel_circ_0019875及novel_circ_0020172的海綿。bta-miR-103屬于高度保守miRNA家族，已被證明bta-miR-103及其靶基因（如Dicer、PANK）可調節脂肪生成[21]。研究發現bta-miR-103可增加細胞內乙酰輔酶A含量，并將乙酰輔酶A導入三羧酸循環，bta-miR-103也可通過與PANK蛋白共生作用調節細胞乙酰輔酶A和脂肪酸水平，影響機體能量代謝和脂肪沉積[22]。在山羊乳腺上皮細胞中，bta-miR-103被認為是circ_015003和circ_014121的海綿，過表達bta-miR-103可增加乳脂合成相關基因表達，并抑制與脂肪分解和β-氧化相關基因（如瘦素、AMP活化蛋白激酶亞單位α），導致脂肪滴形成、甘油三酯積累增加[23]；Zhang等研究發現，circARF3可作為bta-miR-103的海綿而發揮作用，bta-miR-103抑制脂肪組織中有絲分裂，而bta-miR-103的海綿circARF3及bta-miR-103的靶基因TRAF3可促進脂肪組織中有絲分裂，即circARF3通過阻斷btamiR-103效應，導致TRAF3基因表達增加，調控小鼠脂肪合成[24]。Najafi-Shoushtari等發現miR-33 a/b為膽固醇控制的miRNA，其對膽固醇平衡的調節為通過與宿主基因協作，進而抑制三磷酸腺苷結合盒轉運體A1（ATP binding cassette transporter A1，ABCAl）轉運蛋白表達[25-26]。當膽固醇增加時，miR-33通過調節高密度脂蛋白消除血液中多余膽固醇[27]，改善動脈粥樣硬化；當膽固醇減少時，miR-33下調ABCAl表達，減少膽固醇外流，升高細胞內膽固醇水平。miR-499和miR-15a可分別通過其靶基因TCF12和FOXO1介導AMPKPGC-1α、TGF-β及PI3K-AKT-mTOR等信號通路調控動物采食，增強機體吸收利用，參與調控動物機體采食及能量代謝[28-29]。盡管上述miRNA與circRNA間共同作用參與眾多脂肪沉積調控通路，但其是否對牛食物攝入、能量代謝等生物學過程具有調控作用還需進一步探討和驗證。

為深入探索circRNAs在能量代謝方面功能，對差異表達circRNAs的靶基因作GO及KEGG富集分析，結果表明靶基因主要富集在AMPK信號通路、mTORx信號通路及肌動蛋白細胞骨架形成通路中。AMPK信號通路在能量代謝中具有核心地位，該通路被激活后可促進葡萄糖攝取和脂質氧化，當能量缺乏時產生能量，減緩能量消耗[30]；mTOR控制細胞代謝，包括氨基酸生物合成、葡萄糖穩態和其他方面[31]，AMPK可通過mTOR通路調節細胞生長和蛋白質翻譯過程，在能量充足條件下，AMPK無活性，但mTOR活躍，而在能量缺乏時，AMPK活性增加導致mTOR活性降低，機體通過加快攝取葡萄糖提供能量。在下丘腦中，AMPK和mTOR信號通路為食物攝入的決定因素，下丘腦AMPK的調節依賴于攝食狀態，攝入食物可降低下丘腦AMPK活性，而禁食導致下丘腦多個區域AMPK活性增加。下丘腦AMPK和mTOR活性改變可導致代謝紊亂，影響食物攝入量和體重及機體能量平衡[32]。

瘦素為一種脂肪細胞衍生的激素，可激活下丘腦mTOR信號通路活性，因而在減少食物攝入和長期控制體重方面發揮關鍵作用[33]；注射瘦素可抑制下丘腦AMPK活性，給予促食欲肽可增加AMPK活性[34]，如在自由采食的大鼠中，AMPK過度表達降低神經肽Y的mRNA表達，從而控制攝食行為[35]。生物信息學分析顯示，circRNAs可能參與脂肪沉積及能量代謝等生物學過程，但關于circRNAs在脂肪沉積及能量代謝功能及circRNAmiRNA-mRNA互作關系還需進一步在細胞水平進行探討與驗證。

4 結論

本研究采用高通量測序技術獲得牛RFI極端表型個體下丘腦組織中circRNAs表達譜和差異表達信息。在樣本中檢測到豐富circRNAs，并統計其基因類型、長度分布及基因組特征。差異表達circRNAs靶基因主要參與細胞質及細胞器組成，多富集于AMPK、mTOR等食物攝入和能量代謝相關信號通路。為后續探索circRNAs參與牛RFI的生物學過程和分子功能，研究circRNAs在哺乳動物食物攝入方面的調控機制提供數據支持。