基于转录组测序鉴定牛RFI表型差异相关circRNAs及其靶基因富集通路

王鹏飞，丁燕玲，杨朝云，马 应，周小南，史远刚，康晓龙

（宁夏大学农学院，银川 750021）

饲料成本在家畜养殖成本中占比最高，提高饲料利用率、降低养殖成本为家畜规模化养殖及提高效益重要环节。剩余采食量（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在哺乳动物食物摄入方面的调控机制提供数据支持。