曹庆雷
疲劳是一种由外周和中枢系统双重因素影响的复杂现象,如糖原的耗竭、心血管劳损等(Neher et al.,2006)。运动性疲劳的机制存在多种假说,包括中枢神经递质失衡、内稳态失调、氨基酸、离子代谢紊乱、自由基增多、能源衰竭和神经内分泌、免疫系统的平衡打破等(张蕴琨等,2006)。Matsui等(2011)报道,长期的低血糖运动会引起脑糖原降低,导致单胺代谢的激活引起中枢疲劳。长时间运动产生的疲劳以中枢神经系统(central nervous sys‐tem,CNS)出现的保护性抑制为主,其中CNS的神经递质参与了疲劳的产生与发展。过去的运动性中枢疲劳研究多集中在脑和脊髓的形态结构与功能方面,而运动性中枢疲劳是如何激活、是否存在于人类运动性中枢疲劳的相关基因群中,尚不清楚。运动性中枢疲劳是机体对外界刺激做出的保护性反应,其调控机制非常复杂,且基因表达存在着时间和组织特异性。因此,采用传统的核酸印记杂交以及RT-PCR等方法难以获得足够的信息,特别是无法了解基因间的相互关系(Fernstrom et al.,2006)。
热预适应是指将实验体短时间内多次放置在高温条件下进行运动,从而使实验体的表皮和体内温度升高,为后续运动或实验做准备,这也成为研究运动疲劳的一种实验手段(Guy et al.,2015;Tyler et al.,2016;Weller et al.,2007)。热预适应作为一种安全有效且操作简便的实验手段,可减轻实验体在实验过程中产生的机体损伤,比如缺血预适应后可保护脑组织、肌肉组织(刘亮等,2001;张银环等,1995)。陈开祥等(2007)研究表明,热预适应可减轻心肌细胞损伤,从而起到一定的保护作用。热预适应作用的主要原理为:环境温度的改变使得皮肤等表皮器官作为接收器,将热应激信号传递至下丘脑等相关区域,下丘脑整合为热敏元信号释放相关激素,影响中枢神经甚至整个机体,引发中枢性疲劳或运动性疲劳。有研究显示,当机体达到中枢性疲劳或运动性疲劳时,会提高多巴胺(dopamine,DA)、去甲肾上腺素(norepinephrine,NE)的活性来抑制中枢性疲劳,提高运动耐受力,从而延长运动时间(张念云 等 ,2012;Hasegawa et al.,2008;Takatsu et al.,2010)。此外,机体还会产生一系列的反应,例如释放相应神经递质、改变相应受体表达从而调整受体数量,由此进行全身性调节等以适应环境的改变(胡琰茹 等,2015;Fuller et al.,1998;Rodrigues et al.,2003;Walters et al.,2000)。机体在长时间产生此种反应调节后,必定会改变自身原有调节机制,适应环境变化,推测其原有调节机制可能与某些基因的变化相关,变化的基因可能存在于上调基因或改变的转录本中。
在机体维持自身稳态以适应外界环境的过程中,必然会对作为运动承受方及发出方的肌肉组织产生影响。主要的过程为:经下丘脑区域激发产生的神经递质作用于相应效应器以刺激其释放激素与肌肉细胞上的受体结合,影响肌肉细胞内部基因的表达,若肌肉细胞内部失衡,下丘脑分泌相关激素作用于大脑皮层,抑制神经递质释放,使运动能力下降、运动时间减少(卢文彪,2010)。大鼠经热预适应运动后产生的神经和中枢疲劳的信号被大脑皮层接收从而释放神经递质,进行神经-体液调节,使机体逐渐适应外界环境,达到自身稳态平衡(Neher et al.,2006)。大鼠神经递质和激素在全身性调节下趋于平衡状态,有可能适应了实验的热应激条件,推迟了运动疲劳产生的时间(曹庆雷等,2019)。
实验通过转录组分析获得的差异基因与相关参考的注释信息进行比对,对因热预适应产生改变的转录组进行分析总结(窦金焕,2016),了解热预适应所影响的相关转录组,该生物学方法在医学、运动学领域的研究,今后可能广泛应用于大鼠甚至人体。肌肉细胞作为运动的承载者,与运动及一切生命活动的关系最为紧密(谢永涛等,2011),因此本研究将大鼠肌肉细胞作为实验对象,通过高通量测序技术,发现并分析转录组差异表达的基因。
成年雄性SD大鼠24只,体质量230~260 g(购自北京维通利华实验动物公司),饲养于首都师范大学实验动物中心大鼠饲养房,自由饮水、饮食,饲养房温度18℃~24℃,空气相对湿度40%~60%,昼夜节律正常。
1.2.1 热预适应处理
大鼠随机分为对照组和热预适应组,每组12只。将热预适应组大鼠置于42℃水中,使其不间断水浴运动,使大鼠内脏器官温度升高至肛温42℃维持15 min,期间保持大鼠呼吸通畅。运动结束后,将大鼠移出水浴锅,全身擦干,置于常温下1 h,再次水浴运动至肛温42℃维持15 min。每天进行3次,连续7天。对照组大鼠正常饲养在室温25℃的环境中,不做任何处理。所有的实验程序均经首都师范大学动物伦理审查委员会批准,实验方案符合首都师范大学实验动物管理和使用规程。
1.2.2 RNA提取
取热预适应组大鼠和对照组大鼠的肌肉组织,利用Trizol试剂(Life Techonology公司,美国)提取总RNA,置于-80℃条件下保存备用。提取的各RNA样品均以无DNase的DNaseI(TaKaRa公司,大连)进行处理,以消除可能存在的基因组DNA污染,然后分别等量混合3个对照组和热预适应组的RNA样品(30 μg)。使用Oligotex mRNA小量提取试剂盒(Qiagen公司,德国)分离和纯化mRNA。通过 Agilent 2100 Bioanalyzer(Agilent RNA 6000 Nano Kitb Agilent公司,美国)检测total RNA的浓度、RIN值、28S/18S和片段大小。纯度使用紫外分光光度计NanoDrop TM进行检测。
1.2.3 cDNA文库构建及高通量测序
采用mRNA富集法对total RNA进行处理,用带有Ol‐igodT的磁珠富集有polyA尾巴的mRNA以获得片段化的RNA,随机的N6引物进行反转录,合成cDNA。采用琼脂糖凝胶电泳筛选目的片段,通过特异的引物进行PCR,扩增富集所得到的cDNA文库。高通量测序由广州华大基因生物科技有限公司提供测序服务,使用BGISEQ-500进行文库的测序。
1.2.4 数据处理与分析
使用SOAPnuke软件对测序原始数据进行过滤,将其包含的低质量、接头污染以及未知碱基N含量过高的测序片段进行过滤筛选,经数据过滤后获得高质量测序片段。使用Bowtie2将过滤后测序片段与参考序列进行比对,以统计基因比对率,之后再使用RSEM计算基因和转录本的表达水平。使用R语言软件里的cor函数计算每两个样品之间的皮尔逊相关系数。根据基因本体(gene ontology,GO)注释结果以及官方分类,将差异基因进行功能分类,根据KEGG注释结果以及官方分类,将差异基因进行生物通路分类。使用PossionDis算法进行差异基因检测,差异表达基因默认定义为FDR≤0.001且差异倍数在2倍以上。
用于转录组测序RNA先进行质量检测,结果如表1所示。对照组和热预适应组样品的OD260/280值均分别为2.00和2.02;RIN值分别为8.80和8.40,28S/18S值分别为1.20和1.10。这表明本研究所用的样品组织RNA质量较好,满足转录组测序要求。
由于测序仪、试剂、样品质量等均能影响碱基质量测序,因此原始数据包含低质量、接头污染以及未知碱基N含量过高的reads,为保证测序reads的准确性需要确定碱基质量分布。正常情况下,reads中的前几个碱基质量值不高,是因为反转录时随机引物不能很好地结合RNA模板;测序长度的增加在一定范围内有助于提高测序质量,但长度达到一定阈值后,由于测序试剂的消耗,高质量碱基的比例会降低。从整体上看,如果低质量(Quality<20)的碱基比例较低,说明测序质量较好(表2)。
表2 测序数据评估统计Table 2 Satistical Results of Sequencing Data
从上述结果可知,通过测序,共获得了219.72 Mb条的reads,总碱基数为10.99 Gb,每个样品均获得了10 Gb以上的碱基,Q20在97%以上,Q30在90%以上。这表明转录组测序得到的结果质量可靠,可用于后续结果分析。得到过滤后reads之后,使用HISAT序列比对软件将过滤后reads比对到参考基因组序列,结果如表3所示。两个样品分别获得110 454 472和109 271 956个reads,以大鼠基因组为参考进行比对,结果表明,平均每个样品的比对率达到92.53%。比对到参考基因组某一位置过滤后reads平均值为67.50%。两个样品的reads与参考基因组的比对效率较高,说明本次测序结果具有较强可靠性。
表3 参考基因组对比结果统计表Table 3 Statistical Table of Comparison Results of Reference Genome
2.2.1 基因比对与表达
使用Bowtie2将过滤后reads比对到参考序列,之后再使用RSEM转录组分析软件计算基因和转录本的表达水平。以大鼠的表达基因序列为参考基因,对比及表达结果表明,reads比对到参考基因的比率平均为84.26%,reads比对到参考基因唯一位置的比例平均值为65.22%(表4)。说明两组样本的reads对表达基因序列的比对率较高,可以用于后期表达差异分析。从表5可以看出,在基因数、转录本数量上,两组样本间都存在明显差异,热预适应组样本的基因数量较对照组多,特别是新转录本(新转录本是指在参考注释信息中未能找到已知注释信息的转录本,即新转录本可能属于已知基因的新的剪接亚型,或者属于未知基因的新转录本)。这些新的转录本预示新的蛋白编码,这些蛋白很可能在延缓肌肉疲劳的过程中发挥重要作用。
表4 转录组测序数据质量检测Table 4 Quality Testing of Transcriptome Sequencing Data
表5 基因和转录本统计Table 5 Statistics of Gene and Transcripts
2.2.2 差异表达基因的筛选
为了展示基因表达量的分布,对FPKM(FPKM<1、FPKM 1~10、FPKM>10)的3种情况进行了基因数目统计。通过公式计算得到热预适应组及对照组大鼠肌肉组织FPKM值(图1)。对照组FPKM<1的基因占36.63%,FPKM>10的基因占23.25%。热预适应组FPKM<1的基因占34.49%,FPKM>10的基因占28.18%。同时看出,已知基因和新基因在热预适应组中表达量高于对照组,说明大鼠经热预适应后刺激了某些已知基因提高其表达量,或是诱导了某些新基因的表达,参与到延缓疲劳的过程中。
图1 基因表达图Figure 1.Gene Expression Map
2.2.3 差异表达基因检测
为检测大鼠经热预适应后的基因表达差异,进一步使用IDEG6软件,从差异倍数和显著水平进行评估,将差异倍数在2倍以上且显著水平为P<0.05定义为基因表达差异显著。根据火山图(图2)所示,对照组和热预适应组存在较多的差异表达基因。与对照组相比,在热预适应后的肌肉组织中共获得2 267个差异表达基因,其中表达下调的基因555个,表达上调的基因1 712个(图3)。表达上调的基因数目大约是表达下调的基因数目的3倍,表明经过热预适应处理后,跟对照组相比,更多的基因表达上调去适应这种温度的变化。换言之,热预适应处理过程中很多基因需要上调表达。表达量上调最多的基因是α-微管蛋白特异性1c链(tubulin specific 1C chain,TUBA1C),与对照组相比上调432倍,TUBA1C是组成微管的重要蛋白,微管蛋白在细胞运动中具有重要作用,热预适应可以促进微管蛋白的表达升高,表明热预适应也可以改变肌肉组织相关蛋白的表达进而促进对热应激的适应。表达量下调最多的基因是乙酰辅酶A酰基转移酶(acetyl-Co A C-acetyl‐transferase,AACT),是参与细胞营养物质代谢的重要功能酶,与对照组相比,其在热预适应后表达量显著降低533倍,表明热预适应也会影响肌肉组织中相关蛋白酶的表达。这些结果显示,热预适应处理组和对照组差异基因表达模式相差较大。在高温处理后,大量基因表达发生明显的变化,上调的差异表达基因数量较多,这可能参与了大鼠高温应答的互作通路。
图2 差异表达基因火山图Figure 2. Volcanic plot of Differentially Expressed Genes
图3 差异表达基因的统计Figure 3.Statistics of Differentially Expressed Genes
2.2.4 差异表达基因GO功能分析
根据差异基因检测结果,对其GO功能进行分类及富集分析。对GO分析的结果表明,差异表达基因按功能可以分为细胞组分、分子功能和生物过程三大类,生物学过程所占比例最多。在生物学过程中,与生物体代谢调节相关的过程(如代谢过程、生物过程的调节、刺激应答等过程)最多,表明在热预适应后大鼠的肌肉组织中有大量的代谢调控相关基因发生应答反应,从而适应此过程。在细胞组分分类中,细胞部分和细胞器、细胞膜所占比例最多,总体而言,组分类基因变化少于生物学过程基因,通过热预适应处理1周,该过程中肌肉细胞发生响应,但是组分变化的作用低于代谢调控类基因变化的作用。在分子功能中,上调基因发挥的主要功能是催化细胞中相关酶的活性,说明热预适应过程中相关的催化活性蛋白或者结合蛋白发挥了重要作用,这些蛋白可能通过分子功能的改变适应外界环境的变化。此外,在GO功能分类中,与运动相关的生物学过程包括细胞过程、生物调节过程、代谢过程、刺激应答过程、分子或组织的生物合成过程、信号转导过程等。由此可见,经热预适应后的大鼠产生的差异转录组会作用到各种基因类型,共同参与细胞或者机体适应环境。
根据差异基因检测结果,对GO功能进行分类及富集分析,GO功能对应的差异表达基因上、下调统计结果见图4。通过“基因本体功能富集”模块对差异基因进行富集分析,对每类基因数量进行统计,筛选出PDF≤0.01的通路。大鼠在热预适应后整体基因都处于上调状态,下调基因数目较少,可以看出,在生物过程中,上调基因的GO terms与细胞转化、生物调节及代谢过程有关,主要参与细胞的生长发育,下调基因主要参与基因调控。在细胞组分中,上、下调基因主要参与生成细胞器、细胞膜等细胞结构的过程。在分子功能中,上调基因主要参与催化细胞中相关酶的活性与结合活性,下调基因参与细胞器的构成。应激反应中的差异表达基因富集较为显著,表明大鼠肌肉在热预适应过程中产生了大量参与调控代谢的活动,为后续延缓运动疲劳提供了可能性。
图4 差异表达基因GO功能分类统计Figure 4.GO Function Classification Statistics of Differentially Expressed Genes
2.2.5 差异表达基因KEGG富集分析
根据差异基因检测结果,我们对其进行KEGG生物通路分类以及富集分析。将基因根据参与的KEGG代谢通路分为6个部分:细胞过程、环境信息处理、遗传信息处理、人类疾病(仅限动物)、代谢、有机系统(图5)(袁姣,2015)。差异基因在细胞过程中的主要功能是运输和分解代谢,在环境信息处理中参与信号传导和信号分子交互过程,在遗传信息处理中主要参与遗传信息的折叠、排序和退化以及翻译过程,在代谢方面主要参与脂质代谢,在有机系统中主要参与免疫系统和内分泌系统。由此可见,经热预适应后的大鼠转录组发生的改变,主要体现在提高细胞的运输和代谢、增强信号分子传导效率、参与遗传信息表达过程和提高脂质代谢的过程中。热预适应对大鼠肌肉组织转录组的影响主要在于提高细胞代谢、增强传导效率,从而延缓运动疲劳的产生(马国强等,2014)。
图5 差异基因的通路分类Figure 5.Pathway Classification of Differential Genes
使用R语言软件中的phyper函数进行富集分析,选取10条差异基因显著的信号通路进行解释(表6)。可以看出,大鼠经热预适应后差异基因主要参与蛋白酶体、HIF-1信号通路、PI3K-Akt信号通路、细胞外基质(extracellular matrix,ECM)受体相互作用、膀胱癌、肝纤维化、血管重建、促进卵巢癌细胞增殖与侵袭等过程。其中,蛋白酶体富集最为显著,其功能主要涉及细胞周期控制、细胞凋亡和应激反应等。
表6 差异表达基因通路功能富集部分结果示例Table 6 Examples of the Results of Pathway Functional Enrichment of Differentially Expressed Genes
将富集明显的30条信号通路中基因的上、下调情况进行了统计。实验样本富集通路差异基因上调的基因主要存在于蛋白酶体和谷胱甘肽(glutathione,GSH)代谢过程。蛋白质是生命功能的体现者,而蛋白酶体直接影响某些蛋白质的更新,其中包括错误折叠蛋白和许多在生命活动中起重要作用的蛋白质,如p53、cyclin等。蛋白酶体对蛋白质的降解通过泛素介导,泛素识别要被降解的蛋白质,然后将这种蛋白质送入蛋白酶体进行降解。这就说明,参与蛋白合成和降解的基因表达都是上调的,换言之,经热预适应后,机体在运动过程中很可能通过提高代谢率缓解疲劳。
GSH作为一种细胞内重要的调节代谢物质,既是甘油醛磷酸脱氢酶的辅基,又是乙二醛酶及丙糖脱氢酶的辅酶,参与体内三羧酸循环及糖代谢,并能激活多种酶,如巯基酶-辅酶等,从而促进糖类、脂肪和蛋白质代谢。同时,GSH通过巯基与体内的自由基结合,可直接使自由基还原成酸性物质,从而加速自由基的排泄,并对抗自由基,减轻其对重要脏器的损害。Haddad等(2002)研究发现,GSH参与了脂多糖诱导的细胞因子转录的调节及I-κB/NF-κB信号通路的调节。NF-κB是具有多向性调节作用的蛋白质分子,参与调控多种因子的基因表达,在免疫调控、炎症、应激反应及细胞凋亡中起重要作用。NF-κB还可进行反馈调节,通过迅速合成I-κB终止转录。NF-κB表达广泛,功能及调节较为复杂,而且能够诱导NF-κB和DNA结合的诱导物很多,对应激反应、细胞黏附和增生均有重要的调节作用。
本研究对热预适应的大鼠肌肉组织转录组文库进行分析,经测序得到17 951个基因,已知的基因17 616个,预测的新基因435个;同时,测序得到的9 825个新转录本中,8 325个属于已知蛋白编码基因的新的可变剪接亚型,435个属于新的蛋白编码基因,剩下的1 065个属于长链非编码RNA。在热预适应后的大鼠肌肉组织中共获得2 267个差异表达基因,其中下调基因555个,上调基因1 712个。进一步将这些差异基因进行GO功能注释和KEGG代谢途径注释。GO功能注释分类表明,肌肉组织不同生理状态差异表达基因主要涉及细胞相关类别(细胞组分、细胞器等)、膜相关类别(包括膜、膜部分)、分子绑定、催化活性、分子传感器活性、细胞过程、代谢过程等。这些结果表明,参与大鼠肌肉生理活动的差异基因除了集中于细胞内,还广泛分布于生物调节和代谢过程、刺激应答过程、多生物体过程或发育过程、生物过程正反馈调节、细胞组分或生物合成信号通路等。在生物过程中,上调基因的GO terms与细胞转化、生物调节及代谢过程有关,主要参与细胞的生长发育。下调基因主要参与基因调控。在细胞组分中,上、下调基因主要参与生成细胞器、细胞膜等细胞结构的过程。在分子功能中,上调基因主要参与催化细胞中相关酶的活性、结合活性,下调基因参与细胞器的构成。应激反应中的差异表达基因富集较为显著,表明大鼠肌肉在热预适应过程中存在大量的代谢调控。此外,在肌肉细胞的表达基因中,分子绑定功能类别占有较高比例,进一步证实分子绑定功能类别在肌肉中具有显著表达,推测分子绑定在大鼠肌肉的生理活动中可能具有重要的作用。对差异基因进行通路富集分析,可以确定差异表达基因参与的主要代谢途径和信号传导途径,以及与其他基因的相互作用。本研究KEGG通路分析表明,差异表达基因共涉及326个通路,其中蛋白酶体、HIF-1信号通路、PI3K-Akt信号通路和ECM受体相互作用等4个通路被显著富集。
蛋白酶体主要参与细胞周期控制、细胞凋亡、应激反应、DNA修复、基因转录、抗原提呈、信号转导、癌症、炎症、神经退行性疾病的发生等(Javitt et al.,2019)。低氧诱导因子(hypoxia inducible factor,HIF)由HIF-α和HIF-β两个亚基组成,在低氧下能激活多种靶基因的转录,在细胞、组织生长发育和生理应激,以及某些病理过程中具有重要的作用。近年来,有关HIF-1在细胞低氧应激时的信号转导通路,特别是HIF-α介导的基因转录调控的研究取得了巨大进展。HIF-1可以上调血管系统相关蛋白的基因表达,如促进血管生成的VEGF及其受体的表达,这些基因的表达能够诱导某些血管收缩因子的产生,如一氧化氮合酶(iNOS)ET-1和血红素氧化酶1(HO-1)等,从而增加血流(王萍萍等,2011),为肌肉运动提供更多氧气,延缓运动疲劳。
PI3K(phosphatidylinositide 3-kinases)是一种胞内磷脂酰肌醇激酶,由调节亚基p85和催化亚基p110构成,且与癌基因的产物相关。PI3K本身具有丝氨酸/苏氨酸(Ser/Thr)激酶的活性,也具有磷脂酰肌醇激酶的活性。AKT又称PKB(protein kinase B),是一种丝氨酸/苏氨酸特异性蛋白激酶,在多种细胞生长过程中发挥关键作用,是PI3K重要的下游分子,包括至少3种形式,分别为AKT1、AKT2和AKT3。它们对于调控细胞的生长、增殖、迁移、存活以及糖代谢起到重要的作用,而且此信号通路只存在于经刺激后的个体(黄秀兰等,2008),这就解释了大鼠在热预适应后,经此信号通路既可以调控肌肉细胞生长,又可增加糖代谢为机体提供能量,从而达到延缓疲劳的作用。
细胞外基质是组成间质和上皮血管中基质的不溶性结构成分,主要有胶原蛋白、弹性蛋白、蛋白多糖和糖蛋白等。研究表明,ECM可影响细胞分化、增殖、黏附、形态发生和表型表达等生物学过程。
已有研究发现,中枢性疲劳与脑内神经递质分泌有关(曹庆雷等,2019;杨瑾,2013),了解热预适应产生的新转录本与脑内神经递质分泌有何关系,或许可以为延缓疲劳产生、研究预防措施和治疗方法有一定的推动作用。有研究显示,热预适应可改善神经递质、延缓运动性疲劳出现(王峰,2008)。虽然改变神经递质的具体机制还不清楚,但实验已证明,热预适应对脑内神经递质活动网络造成了影响,使实验体预先进入调节机制(邓聪颖,2007)。而本研究的实验结果表明,热预适应也可通过外界环境造成的体液调节影响激素分泌,甚至造成新转录本的出现影响实验体本身。但研究结果都表明,外界刺激使得实验体产生了适应环境的改变。在转录组测序中,可变剪接需要用到已知的参考基因组序列和已经注释的可变剪接事件(王正志等,2006),由于注释基因组对可变剪接注释不够全面,因此所得结论和研究的深度并不能完全相符,随着科学技术的进步及实验研究的深入,研究结果也将不断推进。
本研究利用高通量测序技术对大鼠经热预适应后的肌肉组织进行转录组测序和分析,探讨了大鼠不同生理状态下肌肉组织差异表达基因的数量,获得了差异表达基因的功能、分类和代谢通路,为丰富大鼠经运动处理后肌肉组织转录组信息、探究热预适应对基因转录组的影响以及揭示神经递质在运动疲劳过程中的作用机制奠定了基础。