冷驯化下中缅树鼩BAT中主要内分泌信号通路转录组分析

王慧娟， 贾 婷， 朱万龙， 王政昆

(1. 云南师范大学生命科学学院云南省高校西南山地生态系统动植物生态适应进化及保护重点实验室生物能源持续开发利用教育部工程研究中心云南省生物质能与环境生物技术重点实验室， 昆明650500；2. 云南经济管理学院医学院， 昆明650106)

随着高通量测序快速发展，转录组学成为筛选和定量分析差异表达基因的主要研究方法之一[1]，也涉及对小鼠等哺乳动物的脂肪、肌肉和肝脏研究[2]。

胰岛素(insulin, INS)由胰腺β细胞合成和分泌，在调节葡萄糖稳态、代谢，器官生长发育等方面发挥作用[3]。胰岛素受体(insulin receptor, IR)由α、β两个亚基组成，INS与IR结合导致IR酪氨酸激酶的激活和β亚基酪氨酸磷酸化，β亚基是酪氨酸磷酸化胰岛素受体底物1/2(insulin receptor substrate, IRS1/2)的对接位点。IRS1/2对IR的募集激活两个主要下游信号通路：磷脂酰肌醇三激酶(phosphatidylinositol-3-kinase, PI3K)和丝裂原活化蛋白激酶(mitogen-activated protein kinase, MAPK)通路。PI3K信号通路的激活能促进磷酸肌醇依赖激酶1(phosphoinositide dependent kinase-1, PDK1)磷酸化和蛋白激酶B(protein kinase B, PKB/Akt)的激活，进而磷酸化叉头盒蛋白O(forkhead box protein O, FOXO)和糖原合成激酶-3β(glycogen synthase kinase, GSK-3β)等，在产热、脂肪形成、糖原合成和蛋白质合成等方面发挥作用[4]。MAPK级联的激活是对胰岛素调节有丝分裂细胞增殖、分化起关键作用[5]。

甲状腺激素(thyroid hormone, TH)对维持机体能量平衡、生长发育和增殖分化具有重要作用[6]。甲状腺激素受体(thyroid hormone receptor, TR)主要分为TRα(如TRα1、TRα2)和TRβ(如TRβ1、TRβ2)两个亚型，它们在组织发育、细胞生长和代谢中起作用[7]。TH对下游基因的调控机制可分为两类：经典的基因组途径和非经典的非基因组途径[8]。前者是TH通过与TR结合，再与辅激活因子、辅抑制因子相互作用，刺激或抑制基因的转录，参与调节能量平衡[9]。后者包括TH结合TR调节通路活性和TH与膜整合素受体相互作用。T3和TRβ结合调节PI3K信号通路活性，从而诱导下游基因的表达；整合素αVβ3能同时激活PI3K通路和MAPK途径[10]，这种细胞表面受体通过两个位点(S1和S2)结合TH，能激活以上两种通路[11]，在脂质代谢、细胞增殖分化等方面发挥作用。

脂肪组织根据其功能不同主要分两大类：褐色脂肪(brown adipose tissue，BAT)和白色脂肪(white adipose tissue, WAT)[12]。WAT被经典地认为是一个主要的能量储存库，能量在体内以脂质的形式积累；而BAT在适应性产热过程中具有将多余能量转化为热量的显著能力在冷暴露条件下，BAT大量产热，维持机体体温[13]。本研究组先前的研究表明冷驯化下中缅树鼩(Tupaiabelangeri)产热增加，BAT质量和解偶联蛋白-1(uncoupling protein 1, UCP1)含量增加，以及脂肪转化因子表达量等方面会发生变化[14]。

中缅树鼩属于攀鼩目(Scadentia)，树鼩科(Tupaiidae)，树鼩属(Tupaia)，为攀鼩目在国内仅有物种，热带、亚热带为主要分布区，在中国主要分布于云南、海南和贵州等地[15]。由于其与现生灵长类有较近的亲缘关系，在生物医学研究中作为人类代谢性疾病动物模型。先前冷驯化下中缅树鼩的研究主要从个体水平、细胞水平、分子水平分析其冷适应机制，而本文基于RNA-seq研究冷驯化条件下中缅树鼩BAT中基因的差异表达，从而提供基因功能富集等相关信息，实现对基因功能的解析，对深入研究小型哺乳动物冷适应对策具有重要作用。

1 材料与方法

1.1 试验动物

本研究以成年健康的雄性中缅树鼩为试验对象，采集于云南省昆明市禄劝县(25°25′～26°22′N，102°13′～102°57′ E，海拔1 679 m)。在云南师范大学生命科学学院动物房饲养，室温下适应一个月，平均体重为(115.36±3.65) g，根据其体重分为两组：对照组(n=6，温度25 ℃±1 ℃)；冷驯化组(n=6，温度5 ℃±1 ℃)，两组均给予中等光照，水食自取，驯化28 d。驯化后用乙醚麻醉后将其处死，迅速取动物肩胛部的BAT，-80 ℃保存备用。

1.2 转录组测序

根据miRNeasy Mini Kit试剂盒说明书分别提取总RNA后检测样品。分别对12个个体进行cDNA文库构建及文库质控，完成后选择北京百迈客生物公司Illumina novaseq 6000测序平台对文库进行测序。

1.3 测序数据处理

为保证得到高质量的Clean Data，需将测序所获得的原始数据去除含有接头和低质量的读长Reads，同时计算Q20，Q30值和碱基含量来控制测序数据的质量。使用HISAT2将过滤后得到Clean reads与指定的基因组作为参考进行序列比对，参考基因组来源于Ensembl数据库，网址为(http://www.ensembl.org/Tupaia_belangeri/Info/Index?redirect=no)；再进行文库质量评估；对样本中表达的基因进行归一化处理后计算FPKM值，得到基因的表达量。

1.4 差异表达分析

根据对照组和冷驯化组间各样品中基因的表达量筛选出差异表达基因(differentially expressed genes, DEGs)。筛选条件为差异倍数(fold change, FC)≥2，错误发现率(false discovery rate, FDR)<0.01，其中FC表示两样品组间表达量的比值，而FDR是由差异显著性P值校正后得到的。绘制差异表达基因火山图和聚类图。

1.5 差异表达基因功能注释和富集分析

对差异表达基因做基因信息注释，为了解差异表达基因的生物学功能，将差异表达基因进行基因本体(gene ontology, GO)数据库分析，以KS<0.01作为显著性富集标准，筛选出3个分支显著富集到GO term；为进一步解读基因功能，利用京都基因与基因组百科全书(kyoto encyclopedia of genes and genomes, KEGG)数据库将差异表达基因进行通路富集分析，计算富集显著性Q值(Q-value为多重假设检验校正之后的P-value)，找出显著富集的pathways，并筛选出与胰岛素和甲状腺激素信号通路相关的差异基因。再结合FPKM值以及相关文献的报道，筛选出关键基因。

1.6 转录组测序可靠性分析

为验证转录组测序结果的可靠性，选取5个DEGs进行RT-PCR定量，内参基因引物序列为F：GAGAGGGAAATCGTGCGTGAC；R：CATCTGCTGGAAGGTGGACA，差异基因引物设计和实时定量操作步骤详见文献[14]，引物序列见表1。

表1 实时PCR定量引物

2 结果与分析

2.1 转录组数据评估

各个样本RNA总量>3 μg，完整度(RIN值)均≥7，样本检测合格。转录组测序数据统计见表2，各样品中Q30碱基百分比均不小于93.62%，GC含量均在50%左右，样品间重复相关性良好，测序所得数据质量可靠。

表2 各样品数据产出统计

2.2 差异表达分析

根据差异表达基因的FPKM值以及筛选标准，冷驯化组较对照组，共2 879个基因表达存在显著差异，差异表达火山图见图1。与对照组相比，冷驯化组中有1 181个基因上调，1 698个基因下调。图2显示了2 879个差异表达基因在各样品中的表达量，左侧分支越接近说明基因表达量越接近。

2.3 差异表达基因的功能注释和富集分析

将差异表达基因与COG、GO、KEGG、KOG、NR、Pfam、Swiss-Prot、eggNOG数据库进行比对，统计注释到的基因数量(表3)。GO注释结果统计：细胞组分、分子功能和生物学过程分别注释到1 787、1 736和1 759个差异基因(图3)；以富集条件KS<0.01为标准，KS值越小富集越显著，细胞组分、分子功能、生物学过程分别富集到24、28和55个term，主要富集在生物学过程。表4列出了3个分支富集最显著的前10个条目：“生物学过程”有5个与代谢过程相关(GO:0006695、GO:0055114、GO:0006658、GO:0006509、GO:0008299)，“细胞组分”包括细胞突触和神经元，“分子功能”包含催化活性和酶活性。

横轴为log2(FC)，偏离0越远说明差异倍数越大；纵轴为-log10(FDR)，数值越大说明结果越可靠；红点、绿点分别代表差异基因上调、下调，黑点则为无差异基因。图1 差异表达火山图Figure 1 Volcano plot of differentially expressed gene

每列代表一个样本，每行代表一个基因，颜色代表了基因在样品中的表达量水平log10(FPKM+0.000001)。图2 差异表达基因聚类图Figure 2 Analysis of differentially expressed genes

图3 差异基因GO分类注释统计Figure 3 GO annotation classification statistics of differential expression genes

表3 注释的差异表达基因数量统计表

表4 差异表达基因GO富集分析最显著的前10条目

续表4(Continued Table 4)

通过KEGG注释系统分析，本研究发现差异表达基因参与细胞过程(cellular processes)、环境信息处理(environmental information processing)、代谢(metabolism)、生物系统(organismal systems)等6大类以及47小类(图4)；KEGG显著富集到20条pathways，其中胰岛素和甲状腺激素信号通路差异表达基因的表达情况见图5，在胰岛素信号通路中筛选出的关键基因有PIK3CB、PIK3CD、PDPK1、AKT2、ACACA和SLC2A4等；甲状腺激素信号通路中筛选出的关键基因有THRB、PIK3CB、PIK3CD和PDPK1等。

图4 差异表达基因KEGG分类图Figure 4 KEGG classification map of differentially expressed genes

(a)胰岛素信号通路；(b)甲状腺激素信号通路。图5 通路的DEGs和log2FCFigure 5 DEGs and log2FC were significantly enriched by thetwo pathways

2.4 测序结果可靠性

对5个DEGs进行验证，Real-time PCR与RNA-Seq比较结果见图6，两者之间的上下调表达趋势一致，说明转录组测序质量可靠。

图6 RNA-Seq与 Real-time PCR结果比较Figure 6 Comparison of the RNA-Seq and Real-time PCR results

3 讨论

胰岛素信号通路主要参与调节机体的糖和脂肪代谢，INS激活PI3K通路会导致葡萄糖转运体-4(glucose transporter 4, GLUT4)从细胞内转移到细胞膜表面，促进脂肪对葡萄糖摄取，从而抑制糖原合成[16-17]。PDK1和AKT都是丝/苏氨酸激酶，属于AGC蛋白激酶家族成员，脂肪细胞缺失PDK1可显著抑制胰岛素诱导的PKB/AKT的激活，从而抑制葡萄糖转运体GLUT4的细胞膜转位和葡萄糖的摄取[18]。Hosooka等[19]研究表明脂肪细胞特异性缺乏PDK1的小鼠，由胰岛素诱导的葡萄糖摄取和脂肪生成被抑制。研究者采用shRNA方式抑制AKT2的表达阻碍了胰岛素诱导的脂肪生成过程[20]。小鼠PKBβ基因座的靶向破坏导致脂肪细胞中胰岛素刺激的葡萄糖摄取受损[21]。本研究中，冷驯化组PDPK1和AKT2基因表达都显著下调，说明中缅树鼩可能通过下调PDPK1和AKT2的表达来抑制脂肪的生成和葡萄糖摄取，其中葡萄糖的摄取被抑制，从而促进糖原的合成，为机体提供能量。溶质载体家族2(soluate carrier family 2,SLC2A4)基因又称GLUT4，据报道SLC2A4基因表达下降降低了脂肪组织对葡萄糖的吸收能力[22]。本研究冷驯化组SLC2A4基因表达显著下调，说明中缅树鼩可能通过下调SLC2A4基因的表达来抑制葡萄糖的吸收，从而促进糖原的合成，为机体提供能量。

INS在调节脂质代谢时由PI3K通路中两个主要转录因子介导：固醇反应元件结合蛋白1c(steriod response element binding protein-1c, SREBP-1c)和FoxO1，能决定乙酰辅酶A羧化酶α(acetyl-CoA carboxylase alpha,ACACA)、脂肪酸合成酶(fatty acid synthase,FASN)等脂肪细胞特异性基因的转录活性[23-24]。ACACA活性和表达水平能够调控脂肪酸的合成速度从而影响脂肪的沉积[25]。本研究中与对照组相比，冷驯化组编码ACC和FAS的基因ACACA和Tupaia_belangeri_newGene_30753表达显著下调，说明中缅树鼩ACACA和Tupaia_belangeri_newGene_30753基因可能在调控脂肪合成方面发挥作用。FoxO1转录因子属于叉头盒O家族的重要成员，干扰FOXO1基因在脂肪细胞(前体细胞)中表达，抑制了脂肪细胞分化，从而阻碍脂肪生成[26-27]。本研究中编码FoxO1的ENSTBEG00000011894基因下调，说明冷驯化下中缅树鼩可能通过下调ENSTBEG00000011894基因来影响脂肪的分化和形成。

TH通过与TR结合的基因组作用能调节脂质与能量代谢，与TRβ结合后能降低胆固醇[9]。TRβ1属于核受体超家族成员，是由THRB基因编码的[28]。Hashimoto等[29]证明TRβ突变的大鼠在甲状腺素刺激下增强降胆固醇效应，有研究表明小鼠脂肪前体细胞中胆固醇含量不足，最终为其他组织或器官提供能量来源[30]。本研究中冷驯化组THRB基因表达下调，说明中缅树鼩BAT中胆固醇含量降低了，从而为体内其他组织或器官提供能量来源。整合素αVβ3属于整合素α链家族成员，有证据证实TH结合整合素αvβ3在ERK1/2信号通路中发挥了促血管生成作用[31]。整合素αV的表达与生存率具有相关性[32]。本研究中冷驯化组编码整合素αVβ3的ENSTBEG00000015637基因下调，可能影响中缅树鼩血管生成和生存，中缅树鼩本就属于热带区生存的动物，可能需要调控与存活有关的基因来应对寒冷的环境。蛋白激酶A(protein kinase, PKA)属于丝/苏氨酸蛋白激酶，编码PKA的下调能抑制细胞的增殖[33]，本研究编码PKA的蛋白激酶cAMP激活的催化亚基β(protein kinase cAMP-activated catalytic subunit beta,PRKACB)的下调可能导致冷驯化下中缅树鼩细胞的增殖被抑制。研究表明PI3K抑制剂可减弱T3对脂肪酸合成酶mRNA的诱导效应，从而减少脂肪的生成[34]。冷驯化组编码PI3K的基因PPIK3CB、PIK3CD、ENSTBEG00000008364下调，说明中缅树鼩可能通过下调PPIK3CB、PIK3CD、ENSTBEG00000008364基因减弱脂肪酸合成酶的表达，最终导致脂肪生成减少。Hashimoto等[29]研究表明，T3能通过激活PI3K通路调节SREBP-1c的表达，最终调节脂质的合成。本研究中INS和TH都能通过激活PI3K通路调控SREBP-1c的表达，从而调节细胞内脂肪的生成，说明冷驯化下中缅树鼩BAT中INS和TH在脂质代谢方面可能存在协同作用。胰岛素与甲状腺激素信号在脂质合成代谢方面存在协同效应[35]。

综上表明，胰岛素信号通路和甲状腺激素信号通路在中缅树鼩应对寒冷环境过程中发挥重要作用，中缅树鼩通过调控基因的表达，从而影响能量代谢、生存、血管生成细胞增殖分化来应对寒冷环境。