鄧小書，衛(wèi)秋陽，譚發(fā)銀，郭連安，王嘉，邢康康，陳仕江*

1.重慶市中藥研究院，重慶 400065；

2.重慶中醫(yī)藥學院，重慶 402760

冬蟲夏草Cordyceps sinensis（BerK.）Sacc.是我國傳統(tǒng)名貴中藥材，具有抗腫瘤、免疫調(diào)節(jié)、抗菌等作用[1-2]。隨著消費者對冬蟲夏草認可度的提升和消費需求的日益增長，過度采挖、資源破壞、氣候變化等原因?qū)е乱吧x夏草產(chǎn)量急劇下降，市場供需失衡[3]。近年來，人工培殖冬蟲夏草技術(shù)已取得關鍵性突破，如人工室內(nèi)感染蝙蝠蛾幼蟲形成冬蟲夏草子實體[4]、低海拔規(guī)?；斯わ曫B(yǎng)寄主昆蟲蝙蝠蛾[5]等。隨著高通量技術(shù)的飛速發(fā)展，基于高通量測序的轉(zhuǎn)錄組分析已成為研究功能基因組的重要工具，是挖掘生物學性狀的關鍵功能基因和揭示不同生物學性狀分子機制的重要手段，在菌物類、動物類中藥材研究中均有應用[6-7]。本課題組在人工培殖冬蟲夏草過程中發(fā)現(xiàn)，重慶人工培殖的冬蟲夏草顏色較康定基地培殖的冬蟲夏草更黑。本研究對康定和重慶人工培殖冬蟲夏草僵蟲進行轉(zhuǎn)錄組測序，分析差異表達基因，為進一步探討人工培殖冬蟲夏草蟲體顏色變化提供參考。

1 材料

1.1 樣品

康定人工培殖冬蟲夏草僵蟲和重慶常規(guī)人工培殖冬蟲夏草僵蟲分別由重慶市中藥研究院冬蟲夏草研究所康定實驗基地和重慶實驗基地冬蟲夏草人工培殖，溫度為（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 公司）；反轉(zhuǎn)錄試劑盒HiScript?Ⅱ 1stStrand cDNA Synthesis Kit、SYBR Green Master Mix熒光定量試劑盒（南京諾唯贊生物有限公司）。

2 方法

2.1 RNA提取與cDNA文庫構(gòu)建

參照mirVana? miRNA Isolation Kit 操作流程提取冬蟲夏草僵蟲的總RNA，NanoDrop 2000 型超微量分光光度計和生物分析儀檢測RNA 樣品的純度、濃度和完整性。按照試劑盒TruSeq Stranded mRNA LT Sample Prep Kit操作說明書進行cDNA 文庫構(gòu)建。使用2100 Bioanalyzer 型生物分析儀進行質(zhì)檢，質(zhì)檢合格后于Illumina HiSeq 2500 型高通量測序平臺對文庫進行測序。

2.2 序列篩選與組裝

測序后產(chǎn)生的原始數(shù)據(jù)（raw reads）使用Trimmomatic[8]軟件進行質(zhì)量控制，獲得高質(zhì)量的clean reads。使用Hisat2[9]將clean reads比對到物種參考基因組。參考基因組來源于美國國家生物技術(shù)信息中心（NCBI，https://www.ncbi.nlm.nih.gov/assembly/GCA_000448365.1/）數(shù)據(jù)庫，軟件參數(shù)為默認參數(shù)，通過基因組比對率評估樣本的情況。

2.3 基因功能注釋及差異基因篩選

對得到的unigene 利用基因本體（GO）數(shù)據(jù)庫（http://geneontology.org/）和京都基因與基因組百科全書（KEGG）數(shù)據(jù)庫（http://www.kegg.jp/）進行GO分類和注釋、KEGG 代謝通路分析，進而獲得所有unigene的功能注釋信息。利用DESeq軟件對差異表達基因進行分析，結(jié)合Benjamini-Hochberg 方法對原有假設檢驗得到的顯著性P值進行校正，并最終采用校正后的P值，即錯誤發(fā)現(xiàn)率（FDR）<0.05為差異表達基因篩選的關鍵指標，以降低對大量基因的表達值進行獨立的統(tǒng)計假設檢驗帶來的假陽性[10]。FC 為差異倍數(shù)。以FDR<0.05且|log2FC|>1作為篩選標準，并對篩選到的差異表達基因進行GO 富集及KEGG通路富集分析。

2.4 實時熒光定量PCR（RT-qPCR）分析

為了驗證轉(zhuǎn)錄組數(shù)據(jù)的可靠性，采用真菌18S rRNA[7]作為內(nèi)參基因，對篩選到的9個差異表達基因進行RT-qPCR驗證，設計RT-qPCR特異引物（表1），利用反轉(zhuǎn)錄試劑盒HiScript?Ⅱ 1stStrand cDNA Synthesis Kit 獲得cDNA，采用SYBR Green Master Mix 熒光定量試劑盒，在PCR 儀上進行定量檢測，每個樣品3 個重復，采用2-ΔΔCt法計算差異基因的相對表達量[11]。

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

2.5 數(shù)據(jù)分析

采用Microsoft Excel 2020 進行數(shù)據(jù)整理及作圖，運用SPSS 19.0 軟件對數(shù)據(jù)進行顯著性分析。RTqPCR 驗證的統(tǒng)計分析采用Student′st-test，數(shù)據(jù)值為（）。

3 結(jié)果

3.1 轉(zhuǎn)錄組測序結(jié)果

RNA 樣品質(zhì)量濃度為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），經(jīng)質(zhì)檢表明RNA質(zhì)量符合測序要求，可進行cDNA 文庫構(gòu)建及轉(zhuǎn)錄組測序。測序總計得到raw reads 301.42 M，過濾產(chǎn)生clean reads 293.43 M，clean bases共40.58 G，Q30堿基為94.59%～94.95%，平均鳥嘌呤和胞嘧啶（GC）占比為59.9%（表3）。從測序產(chǎn)出數(shù)據(jù)質(zhì)量評估結(jié)果可見，本次測序數(shù)據(jù)飽和度較高，為后續(xù)轉(zhuǎn)錄本的拼接提供了良好的數(shù)據(jù)支撐。

表2 康定和重慶人工培殖冬蟲夏RNA樣品檢測結(jié)果

表3 康定和重慶人工培殖冬蟲夏轉(zhuǎn)錄組測序數(shù)據(jù)質(zhì)量統(tǒng)計情況

3.2 差異基因表達水平分析

從6個樣本文庫中，獲得21 942條unigenes。6個樣品轉(zhuǎn)錄本的表達豐度主要集中在FPKM（fragments per kilobase million）≥10（高水平表達）和FPKM 1～10（低水平表達）2個區(qū)域，F(xiàn)PKM≤1（極低表達水平）基因數(shù)量則較少[12]，結(jié)果見圖2。

圖2 康定和重慶人工培殖冬蟲夏草的基因表達量分布

3.3 人工培殖冬蟲夏草僵蟲差異表達基因分析

利用DESeq 對重慶和康定人工培殖冬蟲夏草僵蟲轉(zhuǎn)錄組進行比較（圖3），篩選閾值為FDR<0.05且|log2FC|>1，共獲得1575 個差異表達基因。由火山圖、柱狀圖展示差異基因數(shù)目，重慶與康定人工培殖冬蟲夏草僵蟲存在1575 條差異表達基因，這些差異表達基因中上調(diào)表達有626 條、下調(diào)表達有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 個差異表達基因）；在分子功能中，差異表達基因主要富集在結(jié)合（421 個差異表達基因）、催化活性（464個差異表達基因）。差異表達基因的主要功能分類包含細胞過程、代謝過程、催化活性等，說明不同海拔改變了人工培殖僵蟲的代謝和酶的催化功能。

圖4 康定和重慶人工培殖冬蟲夏差異表達基因GO富集分析

3.5 差異表達基因KEGG 通路富集分析

對重慶和康定人工培殖冬蟲夏草僵蟲的差異表達基因進行KEGG 富集分析，注釋到273 條unigenes，差異基因顯著富集（P<0.05）的代謝通路有氨基酸代謝（amino acid metabolism）、折疊-分選-降解（folding,sorting and degradation）、糖代謝（carbohydrate metabolism）、轉(zhuǎn)運和分解代謝（transport and catabolism）、輔酶因子和微生物代謝（metabolism of cofactors and vitamins）。其中，氨基酸代謝、糖代謝、折疊-分選-降解中差異基因數(shù)目改變最為明顯（圖5）。

圖5 康定和重慶人工培殖冬蟲夏差異表達基因KEGG富集分析

3.6 參與冬蟲夏草僵蟲顏色變化的差異表達基因分析

冬蟲夏草為昆蟲、真菌復合體，而昆蟲與食藥用菌顏色變化主要受酚氧化酶、過氧化物酶、活性氧、黑化反應、色素沉積、褐變等影響[13-19]。為挖掘人工培殖冬蟲夏草僵蟲顏色變化成因的關鍵功能基因，在差異表達基因中篩選酚氧化酶類、過氧化物酶類、氧化還原酶類，以及活性氧代謝、黑化反應、色素沉積、褐變相關基因，共篩選到66 個差異表達基因（35 個上調(diào)基因和31 個下調(diào)基因），含46 個氧化還原酶活性相關基因、12 個過氧化物酶體相關基因、3 個單加氧酶活性調(diào)節(jié)相關基因、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，F(xiàn)BP1），進一步分析這9 個關鍵酶基因的表達水平。由圖6可知，與康定人工培殖冬蟲夏草相比，重慶人工培殖冬蟲夏草中patE、TYR、FBP1基因表達上調(diào)，CYP51、fdh、carD、fap1、icdA、tyr1基因表達下調(diào)。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個差異表達基因的表達水平與轉(zhuǎn)錄組測序分析結(jié)果相一致，即與康定人工培殖冬蟲夏草相比，重慶人工培殖冬蟲夏草中patE、TYR、FBP1基因表達上調(diào)，CYP51、fdh、carD、fap1、icdA、tyr1基因表達下調(diào)。由此表明，轉(zhuǎn)錄組測序數(shù)據(jù)可靠。

圖7 康定和重慶人工培殖冬蟲夏草RT-qPCR驗證結(jié)果（,n=3）

4 討論

近年來，隨著高通量測序技術(shù)的不斷發(fā)展，利用轉(zhuǎn)錄組測序成為研究動物類、菌類中藥材功能基因調(diào)控機制的新思路。本研究開展了不同海拔人工培殖冬蟲夏草僵蟲的轉(zhuǎn)錄組測序，以求從分子層面探討人工室內(nèi)培殖冬蟲夏草外觀品質(zhì)的調(diào)控生物學研究基礎。

本研究基于高通量轉(zhuǎn)錄組學技術(shù)對重慶和康定人工室內(nèi)培殖冬蟲夏草僵蟲進行測序分析，在差異表達基因中挖掘冬蟲夏草僵蟲顏色變化相關的關鍵功能基因并分析其表達情況，共篩選到46 個氧化還原酶活性相關基因、12 個過氧化物酶體相關基因、3 個單加氧酶活性調(diào)節(jié)相關基因、3 個黑化反應相關基因、2個ROS相關基因；carD、tyr1、TYR、patE、FBP1、CYP51、fdh、fap1、icdA能顯著富集到KEGG通路中。在這些差異表達基因中，tyr1和TYR屬于酪氨酸酶基因家族基因，而酪氨酸酶廣泛存在于原核、真核微生物及動植物中，屬于多酚氧化酶，是植物及真菌褐變、動物黑色素形成的關鍵酶，參與多個調(diào)控通路，如甜菜素生物合成途徑、酪氨酸代謝、黑化反應等[20-23]。其作用是催化酪氨酸轉(zhuǎn)化為多巴，進而啟動一系列生化反應，最終合成各種色素前體物質(zhì)（包含黑色、白色、黃色）[24-25]。此外，酪氨酸酶是黑色、白色、黃色3 種色素前體物質(zhì)合成過程中的關鍵限速酶，具有酪氨酸羥化酶和多巴氧化酶2 種酶活性，酪氨酸酶活性或其基因表達受影響時，會造成底物大量積累并沿著某一色素合成方向進行反應，不同的比例導致不同的色素合成方向[26-27]。例如，香菇在貯存過程中由于酪氨酸酶Letyr的前體蛋白受蛋白酶水解變成有活性的蛋白，發(fā)生褐變變黑反應，由此表明，酪氨酸酶的活性與褐變相關[28]。Auffret 等[29]利用多階段全轉(zhuǎn)錄組測序追蹤從幼年到養(yǎng)殖珍珠收獲階段與珠母貝白化表型相關的關鍵基因，其研究表明酪氨酸酶Tyr-1基因與黑色素生物合成相關。Bahraman 等[30]研究表明，TYR是催化黑素形成的一種重要限速調(diào)節(jié)酶。黑色素可分為真黑素、類黑色素和異黑色素，通過酚氧化酶激活系統(tǒng)合成，而酪氨酸、γ-谷氨酰-4-羥基苯、兒茶酚等是合成黑色素的前體物質(zhì)。Zhang等[31]對平菇菌蓋顏色調(diào)控機制研究表明，菌蓋中的色素由真黑色素和棕黑色素組成，并且菌蓋顏色不同是由于兩者比例不同造成的，同時，還鑒定得到調(diào)控白黃側(cè)耳菌蓋顏色的關鍵基因——酪氨酸酶基因（PcTYR）。在鱗翅目昆蟲家蠶中，已克隆多個黑色素合成信號通路基因，如酪氨酸羥化酶（TH）[32]、多巴脫羧酶（DDC）[33]、漆酶（yellow/laccase2）[34]等。研究表明，小鼠中出現(xiàn)的白化現(xiàn)象就是由酪氨酸酶基因發(fā)生突變所造成的[35]。同樣，米曲霉酪氨酸酶基因也是黑化過程中的關鍵酶基因[36]。最新研究表明，平菇菌蓋中的色素由真黑色素和棕黑色素組成，并且菌蓋顏色不同是由于兩者比例不同造成的[37]。本研究發(fā)現(xiàn)，酪氨酸酶基因tyr1和TYR在重慶與康定人工培殖冬蟲夏草僵蟲中差異表達，且顯著富集到甜菜素生物合成、異喹啉生物堿生物合成、酪氨酸代謝通路，而這幾個通路最終均能產(chǎn)生黑色素，由此表明tyr1和TYR可能是冬蟲夏草蟲體顏色變化的關鍵酶基因，并為改良人工培殖蟲草品質(zhì)提供理論基礎。