杨永刚, 张梅, 胡耀中, 赵雨, 赵大庆, 幺宝金*
(1. 长春中医药大学附属医院,长春130117; 2. 长春中医药大学创新实践中心,长春130117;3.长春中医药大学吉林省人参科学研究院,长春130117)
鹿茸为鹿科Cervidae动物梅花鹿Cervusnippon或马鹿Cervuselaphus的雄鹿未骨化密生茸毛的幼角,为传统名贵药材,具有壮肾阳、益精血和强筋骨等功效(国家药典委员会,2015)。鹿茸具有循环再生和快速生长两大特点,其生长机制一直是研究热点。快速生长期的鹿茸每天可增加2 cm,伴随着骨化程度的不断加深和最终茸皮的脱落,进而形成彻底钙化的鹿角(Li,2012;Chuetal.,2017)。鹿茸的生长主要由鹿茸顶端的生长中心调控,在结构上主要划分为茸皮、间充质和软骨组织等(Lietal.,2002)。茸皮与鹿茸其他部位的皮肤不同,主要表现为表皮层较厚、皮脂腺较大、立毛肌和汗腺缺失、毛囊组织存在不同发育阶段及无疤痕性损伤修复等(Li & Suttie,2000;Li,2010)。但有关茸皮在鹿茸生长发育过程中的基因表达变化的研究仍然几近空白。因此,研究不同生长时期鹿茸生长中心茸皮的基因表达变化,对于揭示鹿茸再生及快速生长机制具有重要意义。同时,对于皮肤、软骨及骨组织损伤修复及相关疾病的治疗具有重要的指导意义。
近年来,高通量测序技术的不断发展促使生物学领域的研究手段发生了重大的变革,尤其是以基因组学、转录组学和蛋白质组学为代表的技术在科学研究中发挥着巨大作用(Bourgardetal.,2018)。其中,转录组测序(RNA-Seq)更是在基因表达分析研究方面起着重要作用(Hrdlickovaetal.,2017)。RNA-Seq技术具有通量高、灵敏度高和成本低等特点,与基因芯片和传统测序方法相比,可以对基因表达变化进行系统研究,进而对组织或细胞在某种生理和病理条件下的生物学机制进行较为全面的分析和诠释。同时,RNA-Seq除了可以对模式生物进行转录组测序,还可以对无参考基因组的非模式生物进行从头测序(Costa-Silvaetal.,2017;Sudhagaetal.,2018)。本研究采用RNA-Seq技术对快速生长期和骨化期梅花鹿鹿茸生长中心的茸皮进行高通量转录组测序,通过生物信息学分析方法对测序结果进行组装和拼接,筛选不同生长时期茸皮中的差异表达基因,这些研究结果对于揭示鹿茸生长机制及组织损伤修复研究具有重要指导意义。
雄性梅花鹿东北亚种Cervusnipponhortulorum(吉林省双阳鹿乡);DEPC-treated Water(AM9916,Thermo);TRIzol Reagent(15596018,Thermo);氯仿(B0102004,北京化工);异丙醇(B0301007,北京化工);无水乙醇(B0301002,北京化工);TruSeq Stranded mRNA Library Prep Kit(RS-122-2101,Illumina)。
1.2.1 不同生长时期鹿茸生长中心茸皮组织采集选取健康的4周岁雄性梅花鹿6只,分别于快速生长期(生茸后60 d)和骨化期(生茸后90 d)采集两侧鹿茸生长中心组织(距离茸尖5 cm)。按照Li等(2002)的实验方法分离鹿茸茸皮,置于纯化水中反复漂洗去除血液,切成小块后置于液氮罐中冻存。
1.2.2 总RNA提取制备及测序文库构建采用TRIzol(Invitrogen,USA)法从鹿茸生长中心茸皮中提取总RNA,通过Bioanalyzer 2100系统(Agilent Technologies,USA)考察RNA完整性。采用TruSeq Stranded mRNA试剂盒(Illumina,USA)进行测序文库构建:首先,采用带有Oligo(dT)的偶联磁珠(Life Technologies,USA)从总RNA中纯化富集mRNA,并进行mRNA片段化;以片段化的mRNA为模板,经逆转录酶和随机引物合成双链cDNA并经核糖核酸酶处理降解mRNA;最后,经末端修复和接头连接后,通过PCR扩增获得最终的测序文库。
1.2.3 文库测序和差异基因表达分析文库测序在Illumina Hiseq 2000平台(Illumina,USA)上进行,测序结果采用Illumina HCS 1.1进行数据转换和质量评价。去除载体序列和低质量序列(碱基质量值Q≤10的碱基数占整个reads的20%以上),获得clean reads,并采用Trinity(Grabherretal.,2011)进行序列从头组装。组装数据通过BLASTX比对程序与nr(Non-redundant)和Swiss-Prot蛋白数据库比对,进行基因注释,以E<10-5为评价标准。基因表达量分析采用FPKM法(Trapnelletal.,2010)进行计算,差异基因表达分析通过DEGseq(Wangetal.,2010)进行分析,差异表达基因判定标准为:快速生长期和骨化期茸皮组织差异基因表达倍数在2倍以上[即log2(Fold change)≥1或≤-1,Fold change=骨化期FPKM/快速生长期FPKM],且假发现率(FDR)≤0.001。
1.2.4 功能分类及代谢通路富集分析针对筛选出的差异表达基因,通过Blast2GO比对到GO数据库(http://www.geneontology.org/)进行功能分类富集分析,通过KOBAS比对到KEGG数据库(http://www.genome.jp/kegg/pathway.html)进行代谢通路富集分析,采用Fisher检验和FDR校正获得Q<0.05的功能条目或代谢通路判定为显著性富集。
通过高通量RNA-Seq转录组测序,分别从快速生长期和骨化期鹿茸生长中心茸皮中获得4.42 Gb和4.45 Gb的测序数据。其中,测序raw reads分别为45 422 254条(快速生长期)和45 421 562条(骨化期),去除载体序列和低质量序列后,分别获得44 199 056条(快速生长期)和44 530 430条(骨化期)clean reads,测序质量值(Q20值)分别为98.50(快速生长期)和98.39(骨化期)。通过Trinity组装,分别获得44 352条(快速生长期)和43 194条(骨化期)组装序列(Unigenes),平均长度分别为977 nt(快速生长期)和1 135 nt(骨化期)(图1)。
通过BLASTX程序比对,将Unigenes与nr、Swiss-Prot蛋白数据库进行比对,以E<10-5为评价标准,与2个数据库同时比对上的Unigenes共有28 388条。针对这些Unigenes,采用FPKM法进行基因表达量分析,采用DEGseq进行差异基因表达分析,共筛选出5 246条差异表达基因[log2(Fold change)≥1或≤-1;FDR≤0.001]。与快速生长期比较,骨化期显著上调的差异表达基因有 2 627条,显著下调的差异表达基因有2 619条。
将差异表达基因与GO数据库进行比对,差异表达基因显著富集的GO条目特征为:细胞组分分类主要包括细胞部分(cell part)、细胞(cell)和细胞器(organelle)等;分子功能分类主要包括结合活性(binding)和催化活性(catalytic activity)等;生物进程分类主要包括细胞进程(cellular process)、单细胞进程(single-organism process)和代谢进程(metabolic process)等(图2)。
将差异表达基因与KEGG数据库进行比对,差异表达基因显著富集的代谢通路主要包括细胞通讯(cell communication)、信号分子及相互作用(signaling molecules and interaction)、信号转导(signaling transduction)和免疫系统(immune system)等(图3)。
根据功能富集和代谢通路分析结果,推测这些差异表达基因主要为生长因子、转录因子和胶原类成分。因此,根据差异表达基因的注释结果,结合目前已发现和报道的这3类差异表达基因进行进一步的筛选和分析。
图1 快速生长期和骨化期梅花鹿茸皮组装序列长度分布Fig. 1 Length distribution of Unigenes of Cervus nippon antler velvet skins at rapid growth and ossification stages
图2 差异表达基因GO功能富集分析Fig. 2 GO enrichment of differentially expressed genes (DEGs)
图3 差异表达基因KEGG代谢通路富集分析Fig. 3 KEGG pathway enrichment of differentially expressed genes (DEGs)
2.4.1 生长因子类差异表达基因在快速生长期和骨化期茸皮差异表达基因中,共筛选出7种生长因子(表1)。与快速生长期比较,骨化期茸皮中表达量显著上调的生长因子主要包括成纤维细胞生长因子7(Fgf7)和12(Fgf12)以及血管内皮生长因子A(Vegfa);表达量显著下调的生长因子主要包括血小板衍生生长因子A(Pdgfa)、转化生长因子β3(Tgfβ3)、表皮生长因子样蛋白7(Egfl7)和血管内皮生长因子B(Vegfb)。
2.4.2 转录因子类差异表达基因在快速生长期和骨化期茸皮差异表达基因中,共筛选出25种转录因子(表2)。与快速生长期比较,骨化期茸皮中表达量显著上调的转录因子有12种(Pou2f3、Elf2、Bclaf1、Sox6、Runx1、Nfyb、Carf、Nr2f2、Elf1、Atf1、Atf6α和Bbx),表达量显著下调的有13种(Sox12、Ebf3、E4f1、Scx、Tcf4、Nr2f1、Mafk、Ubtf、Crebzf、Elf4、Tfe3、Sox4和Sox18)。
表1 不同生长时期梅花鹿茸皮生长因子类差异表达基因Table 1 Differentially expressed growth factor-related genes of Cervus nippon antler velvet skins at different growth stages
2.4.3 胶原类差异表达基因在快速生长期和骨化期茸皮差异表达基因中,共筛选出6种胶原分子(表3)。与快速生长期比较,骨化期茸皮中表达量显著上调的胶原分子包括Col24a1和Col12a1;表达量显著下调的包括Col13a1、Col27a1、Col3a1和Col1a1。
组织损伤修复是涉及多种细胞和生物大分子参与的复杂的功能组织重建过程,由于人体皮肤、软骨和骨组织的再生能力差,损伤后很难修复和再生,因此,组织损伤修复一直是临床治疗上的难点(Huetal., 2014)。鹿茸生长过程非常独特,具有循环再生和快速生长的特点,鹿茸的茸皮、软骨和骨等组织在结构和功能上与正常组织存在巨大的差别(Lietal.,2014)。因此,探寻鹿茸生长过程中不同组织基因表达变化,将有利于揭示鹿茸独特生长过程的内在分子机制,为组织损伤修复的临床研究提供理论指导。本研究针对快速生长期和骨化期鹿茸生长中心茸皮组织进行深度转录组测序,采用生物信息学分析方法对测序结果进行组装拼接和差异基因表达分析,筛选出7种生长因子、25种转录因子和6种胶原分子,这些调控因子在组织损伤修复过程中发挥重要作用。
表2 不同生长时期梅花鹿茸皮转录因子类差异表达基因Table 2 Differentially expressed transcription factor-related genes ofCervus nippon antler velvet skins at different growth stages
表3 不同生长时期梅花鹿茸皮胶原类差异表达基因Table 3 Differentially expressed collagen-related genes of Cervus nippon antler velvet skins at different growth stages
研究结果表明,鹿茸在不同生长过程中,茸皮中基因表达水平差异显著。在快速生长过程中,茸皮生长主要由生长因子Pdgfa、Tgfβ3、Egfl7和Vegfb调控,并通过上调转录因子Sox12、Ebf3、E4f1、Scx、Tcf4、Nr2f1、Mafk、Ubtf、Crebzf、Elf4、Tfe3、Sox4和Sox18,促进胶原分子Col13a1、Col27a1、Col3a1和Col1a1的表达;而在骨化过程中,茸皮生长主要由生长因子Fgf7、Fgf12和Vegfa调控,并通过上调转录因子Pou2f3、Elf2、Bclaf1、Sox6、Runx1、Nfyb、Carf、Nr2f2、Elf1、Atf1、Atf6α和Bbx,促进胶原分子Col24a1和Col12a1的表达。
在组织损伤修复过程中,生长因子Pdgf、Egf、Vegf、Fgf和Tgf等通过调控细胞的趋化、有丝分裂、血管生成以及细胞外基质的合成,进而调控组织损伤修复(Baeetal., 2014)。研究表明,转录因子Sox4和Sox12同属于SoxC家族,在胚胎发育和组织器官形成过程中发挥重要作用,通过调控细胞的存活、增殖和分化促进组织损伤修复(Lefebvre & Bhattaram,2016;Chang & Hertz,2017;Muetal.,2017;Kavyanifaretal.,2018)。在组织损伤修复过程中,Tgfβ通过激活并上调Scx表达,促进细胞基质分泌进而加速创伤愈合(Sakabeetal.,2018)。此外,研究也表明转录因子E4f1、Tcf4和Elf4在人和小鼠组织损伤修复过程中发挥着重要的作用(Sivinaetal.,2011;Luetal.,2012;Goguet-Rubioetal.,2016)。鹿茸每年可以采集2次,生茸后60 d可以采集1次,俗称头茬茸,鹿茸断面伤口愈合后可以再次长出鹿茸,俗称二茬茸(Suttie & Fennessy,1985;张嵩,李峰,2013)。本研究结果表明,在快速生长期,由于茸皮中促进组织损伤修复的基因显著上调,有利于鹿茸损伤后的创口愈合进而再生出新鹿茸。因此,本研究在梅花鹿茸皮中鉴定出的相关调控因子对于组织损伤修复的临床研究具有重要的指导意义。