基于高通量测序的人工培殖冬虫夏草僵虫颜色变化转录组分析△

邓小书，卫秋阳，谭发银，郭连安，王嘉，邢康康，陈仕江*

1.重庆市中药研究院，重庆 400065；

2.重庆中医药学院，重庆 402760

冬虫夏草Cordyceps sinensis（BerK.）Sacc.是我国传统名贵中药材，具有抗肿瘤、免疫调节、抗菌等作用[1-2]。随着消费者对冬虫夏草认可度的提升和消费需求的日益增长，过度采挖、资源破坏、气候变化等原因导致野生冬虫夏草产量急剧下降，市场供需失衡[3]。近年来，人工培殖冬虫夏草技术已取得关键性突破，如人工室内感染蝙蝠蛾幼虫形成冬虫夏草子实体[4]、低海拔规模化人工饲养寄主昆虫蝙蝠蛾[5]等。随着高通量技术的飞速发展，基于高通量测序的转录组分析已成为研究功能基因组的重要工具，是挖掘生物学性状的关键功能基因和揭示不同生物学性状分子机制的重要手段，在菌物类、动物类中药材研究中均有应用[6-7]。本课题组在人工培殖冬虫夏草过程中发现，重庆人工培殖的冬虫夏草颜色较康定基地培殖的冬虫夏草更黑。本研究对康定和重庆人工培殖冬虫夏草僵虫进行转录组测序，分析差异表达基因，为进一步探讨人工培殖冬虫夏草虫体颜色变化提供参考。

1 材料

1.1 样品

康定人工培殖冬虫夏草僵虫和重庆常规人工培殖冬虫夏草僵虫分别由重庆市中药研究院冬虫夏草研究所康定实验基地和重庆实验基地冬虫夏草人工培殖，温度为（8±1）℃，相对湿度为50%±5%，无光照。上述冬虫夏草的寄主昆虫均为小金蝙蝠蛾Thitarodes xiaojinensis（Tu.），侵染菌种均为冬虫夏草菌Ophiocordyceps sinensi，冬虫夏草僵虫采集时期均为冬虫夏草菌丝将虫体完全包裹，2 种人工培殖冬虫夏草僵虫均采用3个生物学重复（分别命名为KD1、KD2、KD3、CQ1、CQ2、CQ3），康定和重庆人工培殖冬虫夏草僵虫表型见图1。

图1 康定和重庆人工培殖冬虫夏草僵虫表型

1.2 仪器

NanoDrop 2000 型超微量分光光度计（美国Thermo Scientific 公司）；2100 Bioanalyzer 型生物分析仪（美国Agilent Technologies 公司）；HiSeq 2500型高通量测序仪（美国Illumina公司）；CFX Connect型荧光定量聚合酶链式反应（PCR）仪（美国Bio-Rad公司）。

1.3 试药

mirVana™ miRNA Isolation Kit（美国Ambion 公司）；TruSeq Stranded mRNA LT Sample Prep Kit（美国Illumina 公司）；反转录试剂盒HiScript®Ⅱ 1stStrand cDNA Synthesis Kit、SYBR Green Master Mix荧光定量试剂盒（南京诺唯赞生物有限公司）。

2 方法

2.1 RNA提取与cDNA文库构建

参照mirVana™ miRNA Isolation Kit 操作流程提取冬虫夏草僵虫的总RNA，NanoDrop 2000 型超微量分光光度计和生物分析仪检测RNA 样品的纯度、浓度和完整性。按照试剂盒TruSeq Stranded mRNA LT Sample Prep Kit操作说明书进行cDNA 文库构建。使用2100 Bioanalyzer 型生物分析仪进行质检，质检合格后于Illumina HiSeq 2500 型高通量测序平台对文库进行测序。

2.2 序列筛选与组装

测序后产生的原始数据（raw reads）使用Trimmomatic[8]软件进行质量控制，获得高质量的clean reads。使用Hisat2[9]将clean reads比对到物种参考基因组。参考基因组来源于美国国家生物技术信息中心（NCBI，https://www.ncbi.nlm.nih.gov/assembly/GCA_000448365.1/）数据库，软件参数为默认参数，通过基因组比对率评估样本的情况。

2.3 基因功能注释及差异基因筛选

对得到的unigene 利用基因本体（GO）数据库（http://geneontology.org/）和京都基因与基因组百科全书（KEGG）数据库（http://www.kegg.jp/）进行GO分类和注释、KEGG 代谢通路分析，进而获得所有unigene的功能注释信息。利用DESeq软件对差异表达基因进行分析，结合Benjamini-Hochberg 方法对原有假设检验得到的显著性P值进行校正，并最终采用校正后的P值，即错误发现率（FDR）<0.05为差异表达基因筛选的关键指标，以降低对大量基因的表达值进行独立的统计假设检验带来的假阳性[10]。FC 为差异倍数。以FDR<0.05且|log2FC|>1作为筛选标准，并对筛选到的差异表达基因进行GO 富集及KEGG通路富集分析。

2.4 实时荧光定量PCR（RT-qPCR）分析

为了验证转录组数据的可靠性，采用真菌18S rRNA[7]作为内参基因，对筛选到的9个差异表达基因进行RT-qPCR验证，设计RT-qPCR特异引物（表1），利用反转录试剂盒HiScript®Ⅱ 1stStrand cDNA Synthesis Kit 获得cDNA，采用SYBR Green Master Mix 荧光定量试剂盒，在PCR 仪上进行定量检测，每个样品3 个重复，采用2-ΔΔCt法计算差异基因的相对表达量[11]。

表1 康定和重庆人工培殖冬虫夏草RT-qPCR 引物

2.5 数据分析

采用Microsoft Excel 2020 进行数据整理及作图，运用SPSS 19.0 软件对数据进行显著性分析。RTqPCR 验证的统计分析采用Student′st-test，数据值为（）。

3 结果

3.1 转录组测序结果

RNA 样品质量浓度为0.440 5～1.695 2 μg·μL-1，260 nm 处吸光度（A260nm）/A280nm为2.14～2.17，A260nm/A230nm为1.58～2.37，RNA 完整值（RIN）均高于7，28S/18S 为1.2～2.0（表2），经质检表明RNA质量符合测序要求，可进行cDNA 文库构建及转录组测序。测序总计得到raw reads 301.42 M，过滤产生clean reads 293.43 M，clean bases共40.58 G，Q30碱基为94.59%～94.95%，平均鸟嘌呤和胞嘧啶（GC）占比为59.9%（表3）。从测序产出数据质量评估结果可见，本次测序数据饱和度较高，为后续转录本的拼接提供了良好的数据支撑。

表2 康定和重庆人工培殖冬虫夏RNA样品检测结果

表3 康定和重庆人工培殖冬虫夏转录组测序数据质量统计情况

3.2 差异基因表达水平分析

从6个样本文库中，获得21 942条unigenes。6个样品转录本的表达丰度主要集中在FPKM（fragments per kilobase million）≥10（高水平表达）和FPKM 1～10（低水平表达）2个区域，FPKM≤1（极低表达水平）基因数量则较少[12]，结果见图2。

图2 康定和重庆人工培殖冬虫夏草的基因表达量分布

3.3 人工培殖冬虫夏草僵虫差异表达基因分析

利用DESeq 对重庆和康定人工培殖冬虫夏草僵虫转录组进行比较（图3），筛选阈值为FDR<0.05且|log2FC|>1，共获得1575 个差异表达基因。由火山图、柱状图展示差异基因数目，重庆与康定人工培殖冬虫夏草僵虫存在1575 条差异表达基因，这些差异表达基因中上调表达有626 条、下调表达有949条。

图3 康定和重庆人工培殖冬虫夏差异表达基因的火山图和柱状图

3.4 差异基因GO富集分析

对重庆和康定人工培殖冬虫夏草僵虫的差异表达基因进行GO 富集分析，GO 富集被分为3 类，生物过程（biological process）有20个小类，细胞组分（cellular component）有11 个小类，分子功能（molecular function）有12 个小类（P<0.05）。由图4 可知，在生物过程类别中，差异表达基因主要富集在细胞过程（532 个差异表达基因）、代谢过程（436 个差异表达基因）、单一生物过程（388个差异表达基因）；在细胞组分中，差异表达基因主要富集在细胞（564 个差异表达基因）、细胞部分（561 个差异表达基因）、细胞器（408 个差异表达基因）、细胞器部分（221 个差异表达基因）、膜（206 个差异表达基因）；在分子功能中，差异表达基因主要富集在结合（421 个差异表达基因）、催化活性（464个差异表达基因）。差异表达基因的主要功能分类包含细胞过程、代谢过程、催化活性等，说明不同海拔改变了人工培殖僵虫的代谢和酶的催化功能。

图4 康定和重庆人工培殖冬虫夏差异表达基因GO富集分析

3.5 差异表达基因KEGG 通路富集分析

对重庆和康定人工培殖冬虫夏草僵虫的差异表达基因进行KEGG 富集分析，注释到273 条unigenes，差异基因显著富集（P<0.05）的代谢通路有氨基酸代谢（amino acid metabolism）、折叠-分选-降解（folding,sorting and degradation）、糖代谢（carbohydrate metabolism）、转运和分解代谢（transport and catabolism）、辅酶因子和微生物代谢（metabolism of cofactors and vitamins）。其中，氨基酸代谢、糖代谢、折叠-分选-降解中差异基因数目改变最为明显（图5）。

图5 康定和重庆人工培殖冬虫夏差异表达基因KEGG富集分析

3.6 参与冬虫夏草僵虫颜色变化的差异表达基因分析

冬虫夏草为昆虫、真菌复合体，而昆虫与食药用菌颜色变化主要受酚氧化酶、过氧化物酶、活性氧、黑化反应、色素沉积、褐变等影响[13-19]。为挖掘人工培殖冬虫夏草僵虫颜色变化成因的关键功能基因，在差异表达基因中筛选酚氧化酶类、过氧化物酶类、氧化还原酶类，以及活性氧代谢、黑化反应、色素沉积、褐变相关基因，共筛选到66 个差异表达基因（35 个上调基因和31 个下调基因），含46 个氧化还原酶活性相关基因、12 个过氧化物酶体相关基因、3 个单加氧酶活性调节相关基因、3 个黑化反应相关基因、2个活性氧（ROS）相关基因（表4）。其中，有9 个差异表达基因能显著富集到KEGG 通路中，由表5可知，细胞色素P450 51B MSTRG.2450.3（cytochrome P450 51B，CYP51）、甲酸脱氢酶MSTRG.4165.2（formate dehydrogenase，fdh）、乙醛脱氢酶3B1MSTRG.6110.3（aldehyde dehydrogenase 3B1，carD）、GMC 氧化还原酶MSTRG.6971.2（glucose-methanol-choline oxidoreductase，patE）、假定蛋白CDD83_1433 MSTRG.10284.2（hypothetical protein CDD83_1433，fap1）、异柠檬酸脱氢酶NADP MSTRG.8498.1（isocitrate dehydrogenase NADP，icdA）、酪氨酸酶MSTRG.10500.1（tyrosinase，tyr1）、酪氨酸酶MSTRG.304.1（tyrosinase，TYR）、1,6-果糖二磷酸酶MSTRG.7986.1（fructose-1,6-bisphosphatase，FBP1），进一步分析这9 个关键酶基因的表达水平。由图6可知，与康定人工培殖冬虫夏草相比，重庆人工培殖冬虫夏草中patE、TYR、FBP1基因表达上调，CYP51、fdh、carD、fap1、icdA、tyr1基因表达下调。carD在2种人工培殖冬虫夏草中高表达（图6），参与类胡萝卜生物合成（ko00906，carotenoid biosynthesis）；酪氨酸酶基因tyr1、TYR能富集的KEGG 通路有酪氨酸代谢（ko00350，tyrosine metabolism）、异喹啉生物碱生物合成（ko00950，isoquinoline alkaloid biosynthesis）、甜菜素生物合成（ko00965，betalain biosynthesis）。

表4 康定和重庆人工培殖冬虫夏草66个差异表达基因的名称、编号及功能注释

表5 康定和重庆人工培殖冬虫夏草差异关键基因的名称、编号及KEGG富集通路

图6 康定和重庆人工培殖冬虫夏草僵虫颜色变化相关基因表达分析

3.7 RT-qPCR验证

对筛选到的9 个差异表达基因进行RT-qPCR 验证。由图7可知，9个差异表达基因的表达水平与转录组测序分析结果相一致，即与康定人工培殖冬虫夏草相比，重庆人工培殖冬虫夏草中patE、TYR、FBP1基因表达上调，CYP51、fdh、carD、fap1、icdA、tyr1基因表达下调。由此表明，转录组测序数据可靠。

图7 康定和重庆人工培殖冬虫夏草RT-qPCR验证结果（,n=3）

4 讨论

近年来，随着高通量测序技术的不断发展，利用转录组测序成为研究动物类、菌类中药材功能基因调控机制的新思路。本研究开展了不同海拔人工培殖冬虫夏草僵虫的转录组测序，以求从分子层面探讨人工室内培殖冬虫夏草外观品质的调控生物学研究基础。

本研究基于高通量转录组学技术对重庆和康定人工室内培殖冬虫夏草僵虫进行测序分析，在差异表达基因中挖掘冬虫夏草僵虫颜色变化相关的关键功能基因并分析其表达情况，共筛选到46 个氧化还原酶活性相关基因、12 个过氧化物酶体相关基因、3 个单加氧酶活性调节相关基因、3 个黑化反应相关基因、2个ROS相关基因；carD、tyr1、TYR、patE、FBP1、CYP51、fdh、fap1、icdA能显著富集到KEGG通路中。在这些差异表达基因中，tyr1和TYR属于酪氨酸酶基因家族基因，而酪氨酸酶广泛存在于原核、真核微生物及动植物中，属于多酚氧化酶，是植物及真菌褐变、动物黑色素形成的关键酶，参与多个调控通路，如甜菜素生物合成途径、酪氨酸代谢、黑化反应等[20-23]。其作用是催化酪氨酸转化为多巴，进而启动一系列生化反应，最终合成各种色素前体物质（包含黑色、白色、黄色）[24-25]。此外，酪氨酸酶是黑色、白色、黄色3 种色素前体物质合成过程中的关键限速酶，具有酪氨酸羟化酶和多巴氧化酶2 种酶活性，酪氨酸酶活性或其基因表达受影响时，会造成底物大量积累并沿着某一色素合成方向进行反应，不同的比例导致不同的色素合成方向[26-27]。例如，香菇在贮存过程中由于酪氨酸酶Letyr的前体蛋白受蛋白酶水解变成有活性的蛋白，发生褐变变黑反应，由此表明，酪氨酸酶的活性与褐变相关[28]。Auffret 等[29]利用多阶段全转录组测序追踪从幼年到养殖珍珠收获阶段与珠母贝白化表型相关的关键基因，其研究表明酪氨酸酶Tyr-1基因与黑色素生物合成相关。Bahraman 等[30]研究表明，TYR是催化黑素形成的一种重要限速调节酶。黑色素可分为真黑素、类黑色素和异黑色素，通过酚氧化酶激活系统合成，而酪氨酸、γ-谷氨酰-4-羟基苯、儿茶酚等是合成黑色素的前体物质。Zhang等[31]对平菇菌盖颜色调控机制研究表明，菌盖中的色素由真黑色素和棕黑色素组成，并且菌盖颜色不同是由于两者比例不同造成的，同时，还鉴定得到调控白黄侧耳菌盖颜色的关键基因——酪氨酸酶基因（PcTYR）。在鳞翅目昆虫家蚕中，已克隆多个黑色素合成信号通路基因，如酪氨酸羟化酶（TH）[32]、多巴脱羧酶（DDC）[33]、漆酶（yellow/laccase2）[34]等。研究表明，小鼠中出现的白化现象就是由酪氨酸酶基因发生突变所造成的[35]。同样，米曲霉酪氨酸酶基因也是黑化过程中的关键酶基因[36]。最新研究表明，平菇菌盖中的色素由真黑色素和棕黑色素组成，并且菌盖颜色不同是由于两者比例不同造成的[37]。本研究发现，酪氨酸酶基因tyr1和TYR在重庆与康定人工培殖冬虫夏草僵虫中差异表达，且显著富集到甜菜素生物合成、异喹啉生物碱生物合成、酪氨酸代谢通路，而这几个通路最终均能产生黑色素，由此表明tyr1和TYR可能是冬虫夏草虫体颜色变化的关键酶基因，并为改良人工培殖虫草品质提供理论基础。