長(zhǎng)白山豬苓與陜西豬苓的菌絲體轉(zhuǎn)錄組差異比較分析

任 潔， 李太元， 李艷茹， 梁運(yùn)江， 趙洪顏， 許廣波

(延邊大學(xué)農(nóng)學(xué)院，吉林 延吉 133002)

豬苓[Polyporusumbellatus(Pers.) Fries.]為非褶菌目，多孔菌科，多孔菌屬[1-2]，是一種常用的傳統(tǒng)藥用真菌[3-4]。中醫(yī)中有大量關(guān)于豬苓的記載，“豬苓利水盡，則口益干，而欲其口舌之生津”[5]，豬苓的藥用部位為地下生長(zhǎng)的多年生菌核[6-7]，在臨床上可用于治療小便不利、全身水腫、急性腎炎和其他疾病[8]。菌核中最早分離出來(lái)的活性成分是豬苓多糖，現(xiàn)代醫(yī)學(xué)研究發(fā)現(xiàn)，豬苓多糖還具有抑制腫瘤生長(zhǎng)[9-10]，增強(qiáng)機(jī)體免疫力，降血糖[11-12]等藥理作用。

豬苓在我國(guó)分布廣泛[13]，主要分布于河北、山西、陜西、云南等省區(qū)和東北的長(zhǎng)白山區(qū)，其中以云南的產(chǎn)量最大，陜西的質(zhì)量最好[6]。長(zhǎng)白山豬苓在形態(tài)和生物學(xué)特性上與陜西、山西等地的豬苓有很大差異[14]，長(zhǎng)白山地區(qū)的豬苓菌體表面多褶皺、分枝較多、菌體較小，因其形如雞爪而俗稱雞爪苓[15]。許廣波等[16]提出雞爪苓為長(zhǎng)白山區(qū)特有的野生豬苓種屬。邢曉科等[17]利用nrDNA-ITS區(qū)序列分析表明，長(zhǎng)白山區(qū)豬苓與陜西等地區(qū)的豬苓差異較大。在實(shí)驗(yàn)室培養(yǎng)條件下，陜西豬苓的菌絲體長(zhǎng)勢(shì)緩慢，生長(zhǎng)周期長(zhǎng)，并且很容易黃化；相反長(zhǎng)白山豬苓菌絲體生長(zhǎng)速度比較快，菌絲呈白色，在較短時(shí)間內(nèi)就可以鋪滿整個(gè)培養(yǎng)基表面，還能夠生長(zhǎng)出純培養(yǎng)菌核。

該研究通過(guò)對(duì)長(zhǎng)白山豬苓與陜西豬苓菌絲體進(jìn)行轉(zhuǎn)錄組測(cè)序以及數(shù)據(jù)分析，篩選出2種豬苓菌絲體的差異表達(dá)基因[18]，并對(duì)其進(jìn)行生物信息學(xué)分析[19-20]，旨在探討豬苓的基因表達(dá)特性和種源間差異性，為深入研究豬苓種質(zhì)資源的系統(tǒng)進(jìn)化關(guān)系提供理論依據(jù)。

1 材料與方法

1.1 材料與試劑

1.1.1 試驗(yàn)材料

試驗(yàn)材料為延邊大學(xué)微生物實(shí)驗(yàn)室保存的長(zhǎng)白山豬苓和陜西豬苓的菌絲體。

1.1.2 主要試劑與儀器

主要儀器：Qubit2.0熒光計(jì)，微型漩渦混合儀，臺(tái)式高速低溫離心機(jī)，電泳儀，生物電泳圖像分析系統(tǒng)，PCR儀。

1.2 方法

1.2.1 總RNA的提取與檢測(cè)

利用Trizol法提取2種豬苓菌絲體樣品的總RNA，Qubit2.0檢測(cè)RNA濃度，1%瓊脂糖凝膠檢測(cè)RNA完整性以及基因組污染情況[21]。

1.2.2 文庫(kù)構(gòu)建及測(cè)序

1) 利用Qubit2.0 RNA檢測(cè)試劑盒對(duì)Total RNA 精確定量，以確定文庫(kù)構(gòu)建所加入總RNA的量。

2) 利用信使RNA 3′端polyA結(jié)構(gòu)及相關(guān)分子生物學(xué)技術(shù)，對(duì)6個(gè)樣本完整總 RNA進(jìn)行mRNA分離、片段化、雙鏈cDNA合成、cDNA片段化學(xué)修飾、磁珠純化及片段化分選、文庫(kù)擴(kuò)增等處理。經(jīng)過(guò)檢測(cè)和質(zhì)控后，最終獲得能夠適用于Illumina平臺(tái)的測(cè)序文庫(kù)。

1.2.3 測(cè)序數(shù)據(jù)處理

將構(gòu)建好的文庫(kù)送上海生工進(jìn)行測(cè)序，然后根據(jù)測(cè)序結(jié)果，整理、過(guò)濾原始數(shù)據(jù)，獲取優(yōu)質(zhì)序列，并進(jìn)行對(duì)比和統(tǒng)計(jì)分析[22]。

2 結(jié)果與分析

2.1 RNA質(zhì)量檢測(cè)

對(duì)長(zhǎng)白山豬苓以及陜西豬苓菌絲體6個(gè)樣品的總 RNA 進(jìn)行濃度檢測(cè)(表1)，OD260/OD280均在1.8～2.2之間，瓊脂糖凝膠電泳檢測(cè)結(jié)果如圖1所示，從6個(gè)樣品中提取的總RNA都具有清晰的28 S和18 S條帶，且比值約為2∶1，證明所提取的RNA質(zhì)量良好，這與轉(zhuǎn)錄組數(shù)據(jù)庫(kù)的測(cè)序要求一致，可以進(jìn)入下一步試驗(yàn)。

表1 總RNA濃度檢測(cè)

注：M為Marker

2.2 轉(zhuǎn)錄組測(cè)序數(shù)據(jù)組裝

將6個(gè)樣品的原始數(shù)據(jù)進(jìn)行混合，使用Trinity將clean數(shù)據(jù)de novo組裝成轉(zhuǎn)錄本，參數(shù)min_kmer_cov 2，其余默認(rèn)。對(duì)Trinity拼裝得到的轉(zhuǎn)錄本Transcript去冗余，取每個(gè)轉(zhuǎn)錄本聚類中最長(zhǎng)的轉(zhuǎn)錄本作為Unigene，以此作為后續(xù)分析的參考序列。

拼接結(jié)果如表2所示：組裝得到133 662個(gè)Transcript，平均長(zhǎng)度為1 870.76 bp， N50長(zhǎng)度為3 136 bp。其中，長(zhǎng)度≥500 bp的Transcript有97 327個(gè)，占總體的72.82%；≥1 000 bp的Transcript有77 431個(gè)，占總體的57.93%。

在Transcript的數(shù)據(jù)基礎(chǔ)上，進(jìn)一步進(jìn)行組裝，共獲得20 999個(gè)Unigene，平均長(zhǎng)度為1 235.7 bp，N50為3 135 bp。其中，長(zhǎng)度≥500 bp的Unigene有9 368個(gè)，占總體的44.61%；≥1 000 bp的Unigene有6 527個(gè)，占總體的31.08%。

表2 拼接結(jié)果統(tǒng)計(jì)

注：N50/N90的定義為按照長(zhǎng)度將拼接轉(zhuǎn)錄本從大到小排序，累加轉(zhuǎn)錄本的長(zhǎng)度，到不小于總長(zhǎng)50%/90%的拼接轉(zhuǎn)錄本的長(zhǎng)度。

2.3 與參考基因組比對(duì)結(jié)果

將拼接得到的轉(zhuǎn)錄本作為參考序列，使用Bowtie2將質(zhì)控后的測(cè)序序列與參考序列進(jìn)行比對(duì)，通過(guò)RSeQC統(tǒng)計(jì)比對(duì)結(jié)果。

測(cè)序序列經(jīng)過(guò)測(cè)序數(shù)據(jù)過(guò)濾后的高質(zhì)量測(cè)序標(biāo)簽(Total reads)分別為50307178、53538364、55874460、55892258、52633804和50705500。將長(zhǎng)白山豬苓和陜西豬苓菌絲體的6個(gè)文庫(kù)中的所有高質(zhì)量測(cè)序標(biāo)簽與參考基因序列進(jìn)行比對(duì)，顯示在參考序列上有唯一比對(duì)位置的測(cè)序序列數(shù)分別為2020634、2139209、2109861、3053043、3200648和2745108，分別占高質(zhì)量標(biāo)簽總數(shù)的4.02%、4.00%、3.78%、5.46%、6.08%和5.41%。整段比對(duì)到外顯子的測(cè)序序列數(shù)為2020634、2139209、2109861、3053043、3200648和2745108(表3)。

表3 與參考基因組比對(duì)結(jié)果統(tǒng)計(jì)

注：Total Reads: 測(cè)序序列經(jīng)過(guò)測(cè)序數(shù)據(jù)過(guò)濾后的數(shù)量統(tǒng)計(jì)(Clean data)；Total Mapped：能定位到參考序列上的測(cè)序序列的數(shù)量的統(tǒng)計(jì)；Mutiple mapped：在參考序列上有多個(gè)比對(duì)位置的測(cè)序序列的數(shù)量統(tǒng)計(jì)；Unique Mapped：在參考序列上有唯一比對(duì)位置的測(cè)序序列的數(shù)量統(tǒng)計(jì)；Read1Mapped，Read2Mapped：測(cè)序序列Read-1比對(duì)到參考序列上的數(shù)量統(tǒng)計(jì)(僅計(jì)算Unique Mapped序列)；Mapped to '+'， Mapped to '-'：測(cè)序序列比對(duì)到參考序列上正鏈和負(fù)鏈的數(shù)量統(tǒng)計(jì) (僅計(jì)算Unique Mapped序列)；Non-splice reads：整段比對(duì)到外顯子的測(cè)序序列的統(tǒng)計(jì)；Splice reads：分段比對(duì)到兩個(gè)外顯子上的測(cè)序序列(也稱為Junction reads)的統(tǒng)計(jì)；Reads mapped in proper pairs：雙端reads同時(shí)比對(duì)上的測(cè)序序列統(tǒng)計(jì)。

2.4 Unigene功能注釋

各數(shù)據(jù)庫(kù)及功能注釋所用到的軟件及方法：

CDD、KOG、COG、NR、NT、PFAM、Swissprot、TrEMBL注釋：NCBI Blast+；

GO功能注釋：基于Swissprot和TrEMBL 2部分的蛋白注釋結(jié)果根據(jù)Uniprot的注釋信息得到GO注釋；

KEGG注釋：KAAS，KEGG Automatic Annotation Server。

統(tǒng)計(jì)結(jié)果如表4所示。共有20 999個(gè)基因被注釋到，注釋到不同數(shù)據(jù)庫(kù)之中的基因有重疊。其中注釋到NR數(shù)據(jù)庫(kù)的基因最多，為10 824，匹配比率為51.55%；其次是TrEMBL，注釋到10 603個(gè)基因，匹配比率為50.49%。

表4 基因注釋比率統(tǒng)計(jì)

注：Annotated in Database：在數(shù)據(jù)庫(kù)中成功注釋的基因數(shù)量及其占基因總數(shù)的比例；Annotated in at least one Database：在至少一個(gè)數(shù)據(jù)庫(kù)中成功注釋的基因數(shù)量及其占基因總數(shù)的比例；Annotated in all Databases：在所有數(shù)據(jù)庫(kù)中都成功注釋的基因數(shù)目及其占總基因數(shù)的比例。

2.5 差異表達(dá)分析

對(duì)2種材料進(jìn)行組間基因表達(dá)差異分析，記作SX-VS-CBS，采用DESeq進(jìn)行分析。為了得到顯著差異的基因，將篩選條件設(shè)為：qValue<0.05 且差異倍數(shù)|FoldChange|>2。根據(jù)組間比較得到的差異基因，統(tǒng)計(jì)上下調(diào)基因個(gè)數(shù)詳見表5。

與陜西豬苓菌絲體相比，長(zhǎng)白山豬苓菌絲體共有5 881個(gè)表達(dá)差異基因出現(xiàn)，其中上調(diào)基因有2 616個(gè)，下調(diào)基因有3 265個(gè)。

表5 表達(dá)差異分析統(tǒng)計(jì)結(jié)果表

注：Comp.為差異表達(dá)基因集名稱；Up為上調(diào)基因數(shù)；Down為下調(diào)基因數(shù)；Total為差異表達(dá)基因總數(shù)。

2.6 差異表達(dá)基因GO功能富集

基于GO數(shù)據(jù)庫(kù)的注釋結(jié)果，可以看出長(zhǎng)白山豬苓和陜西豬苓菌絲體差異表達(dá)基因和所有基因的注釋情況(圖2)。分析可知，共有9 389個(gè)基因注釋到GO條目中，參與的細(xì)胞組分?jǐn)?shù)目為22個(gè)；分子功能為20個(gè)；生物學(xué)過(guò)程為26個(gè)。差異表達(dá)基因在GO條目中共注釋到2 470個(gè)，上調(diào)基因和下調(diào)基因分別有970和1 500個(gè)，參與生物學(xué)過(guò)程、細(xì)胞組分以及分子功能的基因分別有1 538、1 751、2 015個(gè)。對(duì)差異基因進(jìn)行GO富集分析發(fā)現(xiàn)，在細(xì)胞組分中，膜條目富集到861個(gè)，膜部分條目富集到777個(gè)，細(xì)胞和細(xì)胞部分條目都富集到 1 228個(gè)，細(xì)胞器條目富集到965個(gè)；分子功能中，催化活性條目富集到1 438個(gè)，綁定條目富集到1 306個(gè)；而在生物學(xué)過(guò)程中，代謝過(guò)程條目富集到1 125個(gè)，細(xì)胞過(guò)程條目富集到1 197個(gè)。

注：橫軸是功能分類，縱軸是分類中基因數(shù)(右)及其占被注釋上基因總數(shù)的百分比(左)。不同顏色代表不同的分類。柱狀圖和坐標(biāo)軸上淺色代表差異基因，深色代表所有基因。

圖2 差異基因GO注釋分類柱狀圖

Fig.2 Differential gene GO annotation classification histogram

2.7 差異表達(dá)基因KOG功能分類

根據(jù)KOG數(shù)據(jù)庫(kù)的注釋結(jié)果，注釋到KOG類別中的基因共有4 976個(gè)，差異表達(dá)基因在KOG類別中共注釋到1 006個(gè)，上調(diào)基因和下調(diào)基因分別有322和684個(gè)(圖3)?？梢钥闯鲩L(zhǎng)白山雞爪苓和陜西豬苓菌絲體差異表達(dá)基因所涉及到的KOG功能類別比較全面，涉及了大多數(shù)的生命過(guò)程。其中，一般功能基因預(yù)測(cè)包含基因數(shù)目最多有172個(gè)；其次是蛋白質(zhì)翻譯后修飾與轉(zhuǎn)運(yùn)、伴侶蛋白有108個(gè)，次生代謝產(chǎn)物的生物合成、轉(zhuǎn)運(yùn)和分解代謝有94個(gè)，信號(hào)轉(zhuǎn)導(dǎo)機(jī)制有82個(gè)，脂質(zhì)運(yùn)輸和新陳代謝有74個(gè)，氨基酸轉(zhuǎn)運(yùn)和代謝有67個(gè)，翻譯，核糖體結(jié)構(gòu)和生物發(fā)生有67個(gè)，能源生產(chǎn)和轉(zhuǎn)換有66個(gè)，碳水化合物的運(yùn)輸和新陳代謝有62個(gè)。而參與細(xì)胞活性的基因數(shù)目最少，只有1個(gè)。

注：橫軸為KOG各功能分類，縱軸為該分類內(nèi)基因個(gè)數(shù)(右)及其占被注釋上基因總數(shù)的百分比(左)。柱狀圖和坐標(biāo)軸上淺色代表差異基因，深色代表所有基因。A，RNA加工和修飾;B，染色質(zhì)結(jié)構(gòu)與動(dòng)力學(xué);C，能量生產(chǎn)與轉(zhuǎn)化;D，細(xì)胞周期控制，細(xì)胞分裂，染色體分裂;E，氨基酸轉(zhuǎn)運(yùn)與代謝;F，核苷酸轉(zhuǎn)運(yùn)和代謝;G，碳水化合物轉(zhuǎn)運(yùn)和代謝;H，輔酶轉(zhuǎn)運(yùn)和代謝;I，脂質(zhì)轉(zhuǎn)運(yùn)和代謝;J，翻譯、核糖體結(jié)構(gòu)、生物發(fā)生;K，轉(zhuǎn)錄;L，復(fù)制、重組和修復(fù);M，細(xì)胞壁/膜/被膜生物發(fā)生;N，細(xì)胞活性;O，蛋白質(zhì)翻譯后修飾與轉(zhuǎn)運(yùn)、伴侶蛋白;P，無(wú)機(jī)離子轉(zhuǎn)運(yùn)與代謝;Q，次生代謝物的生物合成、轉(zhuǎn)運(yùn)和分解代謝;R，一般功能基因預(yù)測(cè);S，未知功能;T，信號(hào)轉(zhuǎn)導(dǎo)機(jī)制;U，細(xì)胞內(nèi)運(yùn)輸、分泌和囊泡轉(zhuǎn)運(yùn);V，防御機(jī)制;W，細(xì)胞外結(jié)構(gòu);Y，核結(jié)構(gòu);Z，細(xì)胞骨架。

圖3 差異表達(dá)基因KOG注釋分類統(tǒng)計(jì)

Fig.3 Differential expression gene KOG annotation classification statistics

2.8 差異表達(dá)基因KEGG功能注釋

根據(jù)KEGG中的途徑類型對(duì)差異表達(dá)基因KEGG的注釋結(jié)果進(jìn)行分類，分類結(jié)果如圖4所示。

KEGG代謝通路可被分為細(xì)胞過(guò)程、環(huán)境信息處理、遺傳信息處理、新陳代謝和有機(jī)系統(tǒng)5大類。注釋到KEGG中的基因共有1 959個(gè)，差異表達(dá)基因在KEGG中共注釋到305個(gè)，一共涉及到5大類33個(gè)小類216個(gè)代謝通路，上調(diào)基因和下調(diào)基因分別有76和229個(gè)。其中細(xì)胞過(guò)程中的細(xì)胞生長(zhǎng)與死亡注釋到23個(gè)差異基因，運(yùn)輸和分解代謝注釋到28個(gè)差異基因；環(huán)境信息處理中的信號(hào)轉(zhuǎn)導(dǎo)注釋到27個(gè)差異基因；遺傳信息處理中的折疊、分揀和降解注釋到27個(gè)差異基因，翻譯注釋到38個(gè)；新陳代謝中的氨基酸代謝注釋到37個(gè)差異基因，碳水化合物代謝注釋到48個(gè)差異基因；有機(jī)系統(tǒng)中的內(nèi)分泌系統(tǒng)注釋到13個(gè)差異基因，神經(jīng)系統(tǒng)注釋到10個(gè)差異基因。

注：橫軸為KEGG各代謝通路，縱軸為注釋到該通路下的基因個(gè)數(shù)(右)及其占被注釋上基因總數(shù)的百分比(左)。柱狀圖和坐標(biāo)軸上淺色代表差異基因，深色代表所有基因。

圖4 差異表達(dá)基因KEGG分類

Fig.4 Differential expression gene KEGG classification

3 討論與結(jié)論

近年來(lái)有關(guān)豬苓的遺傳多樣性、功能基因克隆分析等方面的分子生物學(xué)研究正在逐步展開[23]，但是關(guān)于豬苓轉(zhuǎn)錄組方面的研究鮮有報(bào)道。而作為相同分類地位的多孔菌科樹花屬的大型真菌灰樹花關(guān)于轉(zhuǎn)錄組學(xué)方面的研究則相對(duì)較多。王偉科等[24]人利用Illimina高通量測(cè)序技術(shù)對(duì)灰樹花菌絲體和原基進(jìn)行了全轉(zhuǎn)錄組測(cè)序，共發(fā)現(xiàn)兩者具有顯著差異表達(dá)的基因4 094個(gè)，其中，1 886個(gè)基因在原基中上調(diào)表達(dá)，有2 208個(gè)基因在原基中下調(diào)表達(dá)。本研究采用新一代Illimina高通量測(cè)序技術(shù)對(duì)長(zhǎng)白山豬苓和陜西豬苓菌絲體進(jìn)行了轉(zhuǎn)錄組測(cè)序，通過(guò)對(duì)轉(zhuǎn)錄組原始數(shù)據(jù)進(jìn)行數(shù)據(jù)質(zhì)控共得到133 662個(gè)Transcript和20 999個(gè)Unigene。差異表達(dá)分析結(jié)果表明，長(zhǎng)白山豬苓與陜西豬苓菌絲體相比共有5 881個(gè)表達(dá)差異基因出現(xiàn)，其中在長(zhǎng)白山豬苓菌絲體中上調(diào)表達(dá)的基因有2 616個(gè)，下調(diào)表達(dá)的基因有3 265個(gè)。

楊楊[25]對(duì)灰樹花菌絲體進(jìn)行高通量轉(zhuǎn)錄組測(cè)序，拼接共得到14 347條Unigene。其中，共有7 366個(gè)Unigene注釋到GO 數(shù)據(jù)庫(kù)的54個(gè)類別中；有4 347個(gè)Unigene在 KOG中獲得具體蛋白功能注釋；2 799個(gè)Unigene與KO庫(kù)有同源信息，注釋到255個(gè)代謝通路。本研究中，對(duì)差異表達(dá)基因進(jìn)行GO注釋分類，結(jié)果表明長(zhǎng)白山豬苓與陜西豬苓菌絲體的差異表達(dá)基因一共涉及到3大類56個(gè)分支。在細(xì)胞組分中，差異表達(dá)基因主要分布在膜、膜部分、細(xì)胞、細(xì)胞部分和細(xì)胞器；分子功能中，差異表達(dá)基因主要與催化活性、結(jié)合有關(guān)；而在生物學(xué)過(guò)程中，主要與代謝過(guò)程和細(xì)胞過(guò)程有關(guān)。差異表達(dá)基因KOG注釋分類表明，長(zhǎng)白山豬苓與陜西豬苓菌絲體的差異表達(dá)基因主要與一般功能基因預(yù)測(cè)、蛋白質(zhì)翻譯后修飾與轉(zhuǎn)運(yùn)、伴侶蛋白，次生代謝產(chǎn)物的生物合成、轉(zhuǎn)運(yùn)和分解代謝，信號(hào)轉(zhuǎn)導(dǎo)機(jī)制，脂質(zhì)運(yùn)輸和新陳代謝，氨基酸轉(zhuǎn)運(yùn)和代謝翻譯，核糖體結(jié)構(gòu)和生物發(fā)生，能源生產(chǎn)和轉(zhuǎn)換，碳水化合物的運(yùn)輸和新陳代謝，細(xì)胞內(nèi)運(yùn)輸，分泌和囊泡運(yùn)輸有關(guān)。差異表達(dá)基因KEGG分類表明，長(zhǎng)白山豬苓與陜西豬苓菌絲體的差異表達(dá)基因主要涉及到碳水化合物代謝，翻譯，氨基酸代謝，脂質(zhì)代謝，運(yùn)輸和分解代謝，信號(hào)轉(zhuǎn)導(dǎo)，折疊、分揀和降解等。相對(duì)于陜西豬苓菌絲體，長(zhǎng)白山豬苓菌絲體上調(diào)表達(dá)的轉(zhuǎn)錄本在KEGG中顯著富集的是氨基酸代謝，翻譯，碳水化合物代謝，脂質(zhì)代謝，能量代謝；而在陜西豬苓菌絲體中上調(diào)表達(dá)的轉(zhuǎn)錄本在KEGG中顯著富集的是碳水化合物代謝，運(yùn)輸和分解代謝，翻譯，信號(hào)轉(zhuǎn)導(dǎo)，細(xì)胞生長(zhǎng)和死亡。