楊宇婷， 張 強

(成都中醫(yī)藥大學附屬醫(yī)院，四川成都 610000)

同源四倍體具有產(chǎn)量高、抗逆性強、品質(zhì)優(yōu)良等優(yōu)點，在藥用植物育種方面?zhèn)涫荜P注[1]。多倍體植物細胞內(nèi)的染色體加倍導致細胞核和細胞體積增大，在外形上引起植物形態(tài)發(fā)生明顯變化，葉片變大變厚，莖變粗，植株的生物產(chǎn)量和藥用活性成分均有不同程度的增加[2-3]。與二倍體相比，甜葉菊四倍體植株莖稈粗壯、葉片加厚、葉面積增大、植株變矮，糖苷含量等均發(fā)生顯著變化，在農(nóng)藝性狀上，四倍體甜葉菊株系也表現(xiàn)出更高的品質(zhì)和更強的抗逆性[4-5]。

甜葉菊[Steviarebaudiana(Bertoni) Hemsl.]是菊科甜葉菊屬多年生宿根草本植物，它所含有的甜菊糖苷的甜度是蔗糖的250～300 倍，但所含熱量僅是蔗糖的1/300[6]。甜葉菊四倍體株系的高品質(zhì)和高產(chǎn)量在甜葉菊優(yōu)勢品種選育中應用廣泛。目前，關于四倍體甜葉菊的研究大多集中在四倍體誘導鑒定[7]、生物學特性[8]、遺傳物質(zhì)分子標記[9]等方面。隨著生物學的發(fā)展，下一代高通量測序技術使得基因信息處理速度得到極大提升。RNA測序技術已成功應用于白樺[10]、人工合成蕓薹屬[11]、十字花科[12]等倍性植物株系轉(zhuǎn)錄組測序。為探究甜葉菊同源四倍體的變異機制，本研究以甜葉菊二倍體植株為對照，利用BGI PE150測序平臺和生物信息學分析方法，分別取四倍體和二倍體葉片進行轉(zhuǎn)錄組測序，獲得2組差異表達基因(DEGs)，通過對這些差異表達基因數(shù)目進行功能富集分析，為甜葉菊同源四倍體表型變異形成和有效成分基因調(diào)控機制相關研究提供基礎資料。

1 材料與方法

1.1 試驗材料

以成都中醫(yī)藥大學中藥品種質(zhì)量資源研究所種植的甜葉菊二倍體及其四倍體1年生植株為試驗材料，樣本植株經(jīng)成都中醫(yī)藥大學附屬醫(yī)院楊宇婷藥師鑒定為菊科植物甜葉菊。分別隨機選取二倍體和四倍體甜葉菊栽培苗各5株，采樣時間為2020年5月13日，采集甜葉菊葉片進行轉(zhuǎn)錄組測序，做3次生物學重復，試驗時間為2020年6月11日至7月21日。

1.2 方法

1.2.1 總RNA提取 樣品送至成都百納特生物科技有限公司，采用 Trizol 法提取總 RNA，利用 Nanodrop 2000 檢測 RNA 的純度 (吸光度D260 nm/D280 nm>2.0)、Agilent 2100 檢測 RNA 的完整性(RIN>8.5)，檢測合格后構建 cDNA 文庫。

1.3 轉(zhuǎn)錄組數(shù)據(jù)分析

1.3.1 數(shù)據(jù)處理 通過去除污染的reads、未知堿基N含量>5%的reads和低質(zhì)量的reads，得到高質(zhì)量可用片段(clean reads)。使用Trinity軟件對clean reads進行組裝，之后利用Tgicl[13]對轉(zhuǎn)錄本聚類去冗余，最終得到非重復序列基因(Unigene)作為后續(xù)分析的參考序列。

1.3.2 基因表達水平分析 利用RSEM，調(diào)用bowtie2的比對結果進行統(tǒng)計，得到每個樣本比對到每個轉(zhuǎn)錄本上的reads數(shù)目，并對其進行FPKM(Fragments Per Kilobase per Million bases)轉(zhuǎn)換，來自同一個片段(fragment)的成對末端讀取序列(paired-end reads)計數(shù)為1個fragment，進而得到基因和轉(zhuǎn)錄本的表達水平。

1.3.3 差異基因表達分析 采用log2(Fold Change) ≥2.00，Adj. 假定值(P值)≤0.01或0.001的方式篩選差異基因，基于負二項分布原理并根據(jù)Michael 等所述方法[14]進行差異表達基因檢測?；?GOseq所述方法[15]對篩選到的DEGs進行GO富集分析并進行GO功能統(tǒng)計?；诰┒蓟蚝突蚪M百科全書(KEGG)注釋結果，使用KOBAS，設置參數(shù)——fdr為BH(即使用BH校正)進行Pathway富集分析。

2 結果與分析

2.1 RNA-Seq測序原始數(shù)據(jù)分析

通過過濾去除低質(zhì)量、接頭污染以及未知堿基含量過高的reads，6個樣品(T2-1、T2-2、T2-3、T4-1、T4-2和T4-3)測序得到的clean reads的數(shù)量分別為10 542 231、19 757 350、12 758 156、20 225 185、21 463 543、24 421 653個，Q20的比例≥95.2%，低質(zhì)量(Quality<20)的堿基比例較低，說明測序質(zhì)量較好，可用于后續(xù)分析，結果見表1。利用Trinity軟件對clean reads進行合并組裝，并得到轉(zhuǎn)錄本長度分布圖(圖1)，其中最小長度為143，最大長度為1 321，組裝轉(zhuǎn)錄本N50數(shù)值為1 278，數(shù)值較大，說明組裝效果較好，結果見表2。

表1 過濾后的reads質(zhì)量統(tǒng)計

表2 組裝長度分布統(tǒng)計

2.2 轉(zhuǎn)錄本表達水平分析

利用bowtie2[16]，將每個樣本質(zhì)控后的二代序列比對至參考轉(zhuǎn)錄本序列，多重比對序列的平均比對率則為54.35%，唯一比對上參考基因組的平均比對率為29.01%，結果見表3。利用RSEM，調(diào)用bowtie2的比對結果進行統(tǒng)計，得到每個樣本比對到每個轉(zhuǎn)錄本上的reads數(shù)目，并對其進行FPKM分析，結果見圖2。6個樣品中，轉(zhuǎn)錄本的表達豐度主要集中在FPKM 1～15(低水平表達)、FPKM≤1(極低表達水平)這2個區(qū)域，而FPKM≥15(高水平表達)的基因數(shù)目則較少。

表3 reads與組裝轉(zhuǎn)錄本比對結果

2.3 測序基因(或轉(zhuǎn)錄本)功能注釋

對6個樣本的83 234個轉(zhuǎn)錄本做功能注釋，通過與五大數(shù)據(jù)庫進行比對，結果表明，注釋到最多轉(zhuǎn)錄本數(shù)目的是NR 65 371(78.54%)數(shù)據(jù)庫；其他數(shù)據(jù)庫中注釋到的轉(zhuǎn)錄本數(shù)目分別為21 746(26.13%)、28 399(34.12%)、50 980(61.25%)和48 666(58.47%)，其中未知數(shù)據(jù)庫注釋的轉(zhuǎn)錄本數(shù)目為15 248，占總數(shù)的18.32%，結果見表4。

表4 注釋結果統(tǒng)計

2.3 差異表達基因統(tǒng)計及GO和KEGG分析

2.3.1 差異表達基因統(tǒng)計 四倍體甜葉菊中上調(diào)表達量高于二倍體甜葉菊，基因數(shù)目為1 832，二倍體甜葉菊中下調(diào)表達量高于四倍體甜葉菊，基因數(shù)目為3 948，結果見圖3。將二倍體甜葉菊與其同源四倍體中表達的差異基因做層次聚類分析，其中X軸代表進行聚類分析的差異比對，Y軸代表差異基因；每行表示1個基因，每列表示1個樣品。顏色代表差異倍數(shù)，顏色越紅代表上調(diào)倍數(shù)越大，越藍代表下調(diào)倍數(shù)越大。結果表明，甜葉菊植株加倍后其上調(diào)差異基因數(shù)目顯著大于下調(diào)數(shù)目；從顏色深淺可以看出，上調(diào)差異基因倍數(shù)較下調(diào)倍數(shù)大。說明在甜葉菊植株加倍過程中有相關差異基因參與表達，其中上調(diào)的基因調(diào)控作用更加顯著，其差異基因的層次聚類熱圖見圖4。

2.3.2 差異表達基因GO富集分析 為進一步了解加倍后甜葉菊差異基因的表達情況，對DEGs進行GO富集分析。GO分析結果表明，這些差異表達基因分布在分子功能(molecular function)、細胞組分(cellular component)和生物過程(biological process)三大類中的43個類別，包括代謝過程、細胞過程、細胞組分、結合和催化活性等。其中，差異基因大部分富集在生物過程中，而參與分子功能的差異基因相對較少。在生物學進程中，以參與代謝過程、細胞過程、單有機體過程、對刺激的反應、生物調(diào)控和生物過程調(diào)控的差異基因富集最明顯，參與細胞成分的差異基因主要富集在細胞、細胞組分、細胞器、細胞膜等類別；分子功能中差異基因顯著富集在結合與催化活性類別，結果見圖5。

2.3.3 差異表達基因KEGG分析 通過對甜葉菊加倍后的差異表達基因進行KEGG代謝路徑富集分析，結果顯示，差異表達基因被富集到pathway代謝路徑的6個大類、21個亞類。靠前的代謝途徑分別為翻譯(translation)、折疊(folding)、分類和降解(sorting and degradation)、運輸和分解代謝(transport and catabolism)、碳水化合物代謝產(chǎn)物(carbohydrates metabolites)、氨基酸代謝(amino acid metabolism)、脂質(zhì)代謝(lipid metabolism)和能量代謝(energy metabolism)，結果見圖6。

2.4 簡單重復序列(SSR)信息分析

根據(jù)denovo組裝結果，使用MISA軟件對甜葉菊植株加倍過程中獲得的Unigene進行SSR檢測。SSR檢測結果表明，共檢測出的SSR的位點數(shù)目為10 777，最多的是二核苷酸SSR，數(shù)目為4 386，占SSR總數(shù)的40.70%；三核苷酸SSR次之，數(shù)目為 3 412，占31.66%，五核苷酸最少，占0.78%，數(shù)目為84。具體結果見表5。

表5 SSR的類型、數(shù)量及分布比率

3 討論

多倍體具有器官增大、營養(yǎng)成分增多和抗逆性增強等多個優(yōu)點，使得多倍體成為農(nóng)學和園藝學的研究熱點。目前，多倍體器官巨大性的原因是由于細胞增大而不是細胞數(shù)目增多的說法已被證實，但是細胞體積受倍性調(diào)控的分子機制仍不清楚。本研究以甜葉菊二倍體及其同源四倍體植株為試驗材料，利用BGI PE150測序平臺和生物信息學分析方法，分別取四倍體和二倍體的葉片進行轉(zhuǎn)錄組測序并獲得2組差異表達基因，通過對這些DEGs進行功能富集分析，為甜葉菊同源四倍體表型變異形成機制和代謝途徑研究提供新的線索。

對轉(zhuǎn)錄組數(shù)據(jù)差異表達基因進行分析，篩選得到5 780個差異表達基因，其中1 832個基因上調(diào)表達，3 948個下調(diào)。GO分析表明，四倍體與二倍體的差異表達基因主要集中在生物學過程(biological process)和細胞成分(cellular compenent)2個功能區(qū)，對分子功能(molecular function)的調(diào)控作用不大。在細胞成分功能區(qū)發(fā)現(xiàn)差異基因在細胞、細胞組分和細胞器類別顯著富集，表明倍性水平會影響細胞體積，在此過程中大量差異基因在運輸和分解代謝、翻譯、折疊、分類和降解等代謝途徑顯著富集促進相關蛋白質(zhì)的合成，蛋白質(zhì)行使調(diào)控功能從而增加細胞體積，這也是四倍體甜葉菊葉片變大變厚、營養(yǎng)器官增大的原因。同源多倍體藥用植物育種的一個主要目的是增加次生代謝產(chǎn)物的含量，尤其是藥用成分的增加。研究報道，同源四倍體桔梗葉片的葉綠素a、葉綠素b、胡蘿卜素、總葉綠素的含量比二倍體分別增加36.11%、69.23%、31.71%、38.85%[17]，同源四倍體青蒿的青蒿素含量是二倍體的1.5倍[18]。甜葉菊四倍體葉中的活性成分均高于其二倍體，據(jù)文獻報道，甜葉菊四倍體的甜菊糖苷總含量(13.212%)比二倍體的甜菊糖苷總含量(9.998%)提高了32.15%[5]。甜葉菊多倍體的轉(zhuǎn)錄組及SSR分析結果表明，甜葉菊植株面對一些生物和非生物脅迫時，會出現(xiàn)一系列應激反應，形成多種信號轉(zhuǎn)導途徑。同源四倍體甜葉菊具有更強的抗逆性，在面對寒冷刺激時，甜葉菊四倍體會調(diào)控相關差異基因上調(diào)表達來表現(xiàn)出更強的適應性，本研究中，一些與光合作用相關的差異基因在光合作用、光合生物的固碳作用等路徑富集，顯著增強了甜葉菊四倍體株系的光合效率。本研究共發(fā)現(xiàn)10 777個SSR位點，優(yōu)勢重復序列為二核苷酸，占總SSR位點的40.70%，這一研究結果與李俊仁等的結果[19-20]一致。

本研究以甜葉菊二倍體及其同源四倍體為材料，利用轉(zhuǎn)錄組測序技術并結合二者的生理特征，分析同源四倍體產(chǎn)生的基因表達變化，發(fā)現(xiàn)二者在外觀性狀、生長發(fā)育以及抗逆抗脅迫等方面存在差異，以上結果可為進一步定位關鍵調(diào)控基因及基因工程育種等奠定基礎。