鄭 云，崔芳芳，鄭九洲，楊祥飛，孟林峰，王建革，劉齊元

(1. 江西農(nóng)業(yè)大學農(nóng)學院作物生理生態(tài)與遺傳育種教育部重點實驗室，南昌 330045；2. 江西農(nóng)業(yè)大學林學院，南昌 330045)

【研究意義】植物雄性不育(Male sterility,MS)是植物中普遍存在的一種雄蕊不能形成正常花粉的現(xiàn)象[1]。煙草(NicotianatabacumL.)是重要的經(jīng)濟作物，在實驗室中被廣泛應用研究，在生產(chǎn)上迫切要求我們提供更多的優(yōu)質(zhì)抗病、高產(chǎn)品種，因此煙草雄性不育的利用正備受重視。因此利用 RNA-Seq 技術(shù)分析了K326煙草雄性不育系和保持系不同大小花蕾相關(guān)基因的差異表達情況，有助于進一步揭示其雄性不育形成的分子機理?！厩叭搜芯窟M展】我國自1973年開始就對煙草雄性不育系進行了研究[2]。早期的大量研究[3-6]相對集中于煙草雄性不育形態(tài)和生理方面，而國內(nèi)外對煙草雄性不育分子基礎(chǔ)的研究則不多[7-8]?！颈狙芯壳腥朦c】目前，轉(zhuǎn)錄組測序技術(shù)是一種基于Illumina HiSeq等測序平臺，對樣本所有mRNA進行轉(zhuǎn)錄組測序的技術(shù)[9]。不但擁有通量大、檢測范圍廣等優(yōu)點，而且操作起來比較簡單[10]。應用高通量測序技術(shù)的報道在不同物種、不同敗育類型的植物中已有不少。如譚翀[11]利用RNA-Seq技術(shù)，檢測到超過2300個差異表達基因在蔬菜作物大白菜突變體ftms與野生型‘FT’的花蕾中表達，并推測突變體ftms中β-(1,3)-半乳糖基轉(zhuǎn)移酶的突變，最終造成了雄蕊敗育。戴冬洋等[12]對水果作物甜瓜的不育系花蕾中的雄蕊進行轉(zhuǎn)錄組測序，發(fā)現(xiàn)甜瓜雄蕊中相關(guān)基因的表達量減少與雄性不育之間存在一定的相關(guān)性。陶興林等[13]研究了花椰菜溫敏雄性不育系，獲得67 930個Unigene和2170個差異表達基因。毛偉兵等[14]對楸樹不同育性花芽轉(zhuǎn)錄本進行比較分析，發(fā)現(xiàn)丙酮酸代謝異常后進而抑制油菜素內(nèi)酯的合成，最終導致其育性消失。由此可以看出該技術(shù)已經(jīng)應用在甜瓜、花椰菜和楸樹等高等植物中，但在煙草方面的報道卻很少見?！緮M解決的關(guān)鍵問題】本試驗以煙草K326雄性不育系和保持系的花蕾為研究對象，利用轉(zhuǎn)錄組測序技術(shù)對其測序后的數(shù)據(jù)進行組裝、注釋，篩選調(diào)控煙草雄性不育相關(guān)基因，以期揭示煙草雄性不育的形成機制，為研究這些基因具體的生物學功能提供參考和依據(jù)。

1 材料與方法

1.1 材料

供試煙草品種為煙草雄性不育系MSK326及其保持系K326。供試材料于2020年1月播種，2020年4月移栽，盛花期參照劉齊元等[15]的方法取樣后迅速保存于-80 ℃冰箱，將不育系花蕾按大小分為小蕾(sK326-1)、中蕾(sK326-2)和大蕾(sK326-3)3組，保持系花蕾按大小分為小蕾(CK-K326-1)、中蕾(CK-K326-2)和大蕾(CK-K326-3)3組，不育系和保持系相互對照，試驗設(shè)3個生物學重復。

1.2 方法

1.2.1 總RNA文庫構(gòu)建和測序 文庫構(gòu)建和轉(zhuǎn)錄組測序由上海中科新生命有限公司完成。

1.2.2 測序數(shù)據(jù)分析 測序得到的原始數(shù)據(jù)通過trimmomatic去除原始數(shù)據(jù)(raw data)中包含測序接頭、未知堿基以及低質(zhì)量的reads，最后得到后續(xù)轉(zhuǎn)錄組分析使用高質(zhì)量序列(clean reads)。clean reads與參考基因組進行序列比對利用HISAT2軟件，只有完全匹配的reads才可以用于進一步分析和注釋。

1.2.3 差異表達分析 通過DEGSeq R package(1.12.0)軟件分析不同大小花蕾的差異表達，將3個時期進行比對，獲得CK-K326-1-VS-sK326-1、CK-K326-2-VS-sK326-2、CK-K326-3-VS-sK326-3共3個比較組，根據(jù)P<0.05且|log2FoldChange|>1 篩選其中的基因并將其指定為差異表達基因。

1.2.4 富集分析 為研究煙草雄性不育發(fā)育過程中差異表達基因的功能，分別利用采用 Goseq R package (1.10.0)和 KOBAS(2.0)軟件對差異基因集進行GO功能富集分析和 KEGG 通路富集分析，利用iTAK軟件對差異表達基因進行轉(zhuǎn)錄因子分析。

2 結(jié)果與分析

2.1 RNA-Seq測序數(shù)據(jù)統(tǒng)計分析

測序原始數(shù)據(jù)在經(jīng)過過濾后如表1所示，總共84 101 954～134 977 386個reads用于比對分析，堿基測序質(zhì)量值Q30均大于94%，GC含量在41.84%～42.50%，能夠比對到參考基因組的reads比例在95.22%～96.35%；不育系和保持系小、中、大蕾具有唯一匹配的reads比例分別為89.66%、89.87%、90.71%和86.85%、87.30%、86.46%。利用FPKM值進行定量分析，結(jié)果基因表達水平分布范圍在各樣品間較為一致(圖1-A)。為了進一步對獲得轉(zhuǎn)錄組數(shù)據(jù)的可靠性、合理性以及可重復性進行驗證，基于各重復所有基因的FPKM值計算組內(nèi)和組間相關(guān)性系數(shù)的平方(R2)。如圖1-B所示，煙草雄性不育花樣本各重復間的R2分別為0.955(sK326-1/sK326-2)、0.929(sK326-1/sK326-3)、0.922 (sK326-2/sK326-3)，保持系花蕾樣本各重復間的R2分別為0.913 (CK-K326-1/CK-K326-2)、0.881 (CK-K326-1/CK-K326-3)、0.882(CK-K326-2/CK-K326-3)，各樣本重復間相關(guān)系數(shù)的平方均大于0.88，表現(xiàn)出較好的一致性，達到生物學分析的要求。基于相似性和主成分分析(圖1-C)對每個樣本進行區(qū)分和聚類，結(jié)果表明重復樣品間相似度高，重復性好，達到了試驗要求。綜合上述分析，說明轉(zhuǎn)錄組數(shù)據(jù)的質(zhì)量和可靠性較高，可以進行下一步的分析。

表1 不同大小花蕾煙草雄性不育系和保持系的轉(zhuǎn)錄組測序數(shù)據(jù)

2.2 不同時期的煙草雄性不育差異表達基因

由圖2可知，藍點的數(shù)量顯著多于紅點，說明這些下調(diào)基因的表達量顯著高于上調(diào)表達基因。從圖3中可以看出煙草雄性不育系和保持系在中蕾時期的差異表達基因數(shù)目最多，為19 931個；大蕾次之，為10 981個；小蕾最少只有7881個。在小蕾、中蕾和大蕾的差異表達基因中，其中上調(diào)表達的基因數(shù)分別為1718、4992和3925，下調(diào)表達的基因數(shù)目分別為6163、8399和7059個，下調(diào)表達基因數(shù)約為上調(diào)表達基因數(shù)的3.6、1.6和1.8倍，說明煙草K326雄性不育系中的基因通過下調(diào)表達來控制其育性。采用歐氏距離計算基因間相對表達量距離，基于富集分析篩選的煙草雄性不育相關(guān)差異表達基因進行層次聚類，通過聚類圖譜(圖4)可以初步了解與煙草雄性不育相關(guān)的差異表達基因大致可分為5類。

2.3 差異表達基因的GO功能富集分析

通過GO功能富集，煙草雄性不育系中的11 062個差異表達基因被富集于生物過程、細胞組分和分子功能3個部分。其中在生物過程中對刺激的反應term中富集到的差異表達基因最多，有7個；生長素的反應、生長素激活的信號通路、細胞對生長素刺激的反應和信號轉(zhuǎn)導系統(tǒng)等term富集到差異表達基因最少，只有2個；在細胞組分中富集到分子功能term的差異表達基因最多，有18個；富集到轉(zhuǎn)移酶活性term的差異表達基因次之，有5個；在分子功能中富集到差異表達基因都是2個(圖5)，此處只列出最顯著富集的35個功能條目。

2.4 煙草雄性不育差異基因的KEGG通路富集分析

利用KEGG數(shù)據(jù)庫富集出煙草雄性不育發(fā)育過程中表現(xiàn)出顯著性差異的通路，探尋主要代謝和信號轉(zhuǎn)導途徑中差異表達基因所參與的情況。結(jié)果(圖6)顯示最顯著富集的20個通路，其中，代謝途徑(ko01100)差異表達基因數(shù)量最多有717條(表2)；次生代謝物生物合成代謝通路(ko01110)次之，有408條；亞油酸代謝(ko00591)和甲狀腺激素合成(ko04918)相關(guān)的差異表達基因最少，分別為14和13條。并且除玉米素生物合成(ko00908)外，其于代謝通路途徑中的差異表達基因上調(diào)都比下調(diào)少。

A. 樣品FPKM值；B. 樣品相似性聚類分析；C. 樣品 PCA 主成分分析A. FPKM values of samples; B. Clustering of samples based on sample-to-sample distance; C. Principal component analysis of samples圖1 基于RNA 測序分析的基因表達概況和聚類分析Fig.1 Gene expression profile and cluster analysis based on RNA sequencing analysis

Up 表示上調(diào)，Down表示下調(diào)，Not表示既不上調(diào)也不下調(diào)，下同‘Up’ means up-regulation,‘Down’ means down-regulation and ‘Not’ means neither up-regulation nor dowh-regulation.The same as below圖2 煙草雄性不育小蕾、中蕾和大蕾時期差異表達基因Fig.2 DEGs in the periods of small bud, middle bud and large bud in tobacco male sterility

2.5 轉(zhuǎn)錄因子分析

從圖7可以看出，在小蕾、中蕾、大蕾時期較多的差異基因都屬于ERF、MYB、C2H2、NAC、ERKY、bHLH等轉(zhuǎn)錄因子家族，而在BBR-BPC、GRF、ARR-B、CAMTA、EIL、GeBP等轉(zhuǎn)錄因子家族中的差異基因較少，只有1～2個，編碼MYB和bHLH轉(zhuǎn)錄因子的基因在植物花藥中表達，參與調(diào)控花藥形成和花粉發(fā)育以及植物育性。在小蕾、中蕾、大蕾時期分別有38、43、44個轉(zhuǎn)錄因子家族。轉(zhuǎn)錄因子通過激活或者抑制基因表達在植物次生代謝產(chǎn)物的生物合成過程中發(fā)揮至關(guān)重要的作用。

圖3 煙草雄性不育不同時期的差異基因數(shù)目Fig.3 Numbers of differentially expressed genes in different periods of tobacco male sterility

3 討 論

植物雄性不育與制種之間息息相關(guān)，想要充分利用煙草雄性不育培育優(yōu)質(zhì)品種，那么挖掘新的煙草雄性不育基因，拓寬其來源具有重要的實踐意義。要想對某種動植物在某一狀態(tài)下或者某一時期的相關(guān)基因進行篩選，并且獲得相關(guān)基因表達豐度和剪接方式等信息，那么可以選擇對其進行轉(zhuǎn)錄組測序[16]。本試驗對煙草雄性不育花蕾中的高表達基因進行KEGG通路富集分析表明，在次生代謝物的生物合成、代謝途徑、內(nèi)質(zhì)網(wǎng)中的蛋白質(zhì)加工等途徑中富集，這與王志敏等[17]研究蔬菜作物中茄子的功能雄性不育類似，這些通路通過參與植物能量的供應和營養(yǎng)物質(zhì)的分配過程進而造成植物的育性的改變。而馮小雨等[18]在雄性不育小麥的KEGG通路富集中也發(fā)現(xiàn)在苯丙烷生物合成通路中所富集到的差異基因最多。然而，植物苯丙烷合成途徑的次生代謝產(chǎn)物，在參與植物抵抗紫外線輻射和正常的生長發(fā)育過程中同樣起著至關(guān)重要的作用，這也暗示了煙草雄性不育形成機制的復雜性。

紅色表示上調(diào)，藍色表示下調(diào)；橫表示同一差異表達基因在不同樣本中的表達，列表示同一樣本中不同差異表達基因的表達Red means up and blue means down；The horizontal line represents the expression of the same differentially expressed gene in different samples and the column represents the expression of different differentially expressed genes in the same sample圖4 煙草雄性不育相關(guān)差異表達基因聚類分析Fig.4 Cluster analysis of differentially expressed genes rclated to male sterility in toacco

圖5 煙草雄性不育差異表達基因的 GO 功能分類Fig.5 GO functional classification of differentially expressed genes in tobacco male sterility

圖6 KEGG代謝通路分析Fig.6 KEGG mmetabolic pathway analysis

表2 煙草雄性不育轉(zhuǎn)錄組差異表達基因KEGG通路注釋

續(xù)表2 Continued table 2

圖7 煙草雄性不育不同時期轉(zhuǎn)錄組中差異基因?qū)霓D(zhuǎn)錄因子Fig.7 Transcription factors corresponding to differential genes in the transcriptome of different periods of tobacco male sterility

轉(zhuǎn)錄因子在植物生長發(fā)育過程有很重要的作用[19]。MYB和bHLH這兩類轉(zhuǎn)錄因子的基因參與調(diào)控花藥形成和花粉發(fā)育，在植物花藥中表達。因此，對植物育性具有一定的調(diào)控作用。bHLH和MYB轉(zhuǎn)錄因子都是與雄性不育發(fā)育相關(guān)的轉(zhuǎn)錄因子之一。在干燥低溫和高鹽分的環(huán)境中調(diào)控相關(guān)基因的表達，bHLH轉(zhuǎn)錄因子主要參與調(diào)控植物發(fā)育進程和生化過程[20-22]。從煙草雄性不育系和保持系的相互對比中我們發(fā)現(xiàn)此類轉(zhuǎn)錄因子。bHLH這類轉(zhuǎn)錄因子表達異常會使植物絨氈層出現(xiàn)高度液泡化[23]。MYB作為最大的轉(zhuǎn)錄因子之一，參與著植物的生長、代謝、脅迫等多種生理生化過程，在真菌和動物中的數(shù)量相對較少[24]。有在水稻中發(fā)現(xiàn)MYB轉(zhuǎn)錄因子是雄性不育的關(guān)鍵調(diào)控因子[25]。BBR-BPC[26]、GRF[27-28]這兩個轉(zhuǎn)錄因子家族中的差異基因也調(diào)控著植物生長發(fā)育。因此，我們后續(xù)在對煙草育性相關(guān)差異基因做進一步的分析驗證時, 應該尤其重視這些轉(zhuǎn)錄因子的表達情況。

本試驗從分子水平上研究了雄性不育基因及其差異表達，為篩選煙草雄性不育關(guān)鍵基因和功能驗證提供重要的分子基礎(chǔ)，同時對煙草育種具有重要的參考意義。

4 結(jié) 論

本研究通過對煙草K326雄性不育系及其保持系小、中、大蕾時期的花蕾進行轉(zhuǎn)錄組測序，結(jié)果表明測序原始數(shù)據(jù)過濾后，總共得到了84 101 954～134 977 386個clean reads用于比對分析，堿基測序質(zhì)量值Q30均大于94%，GC含量在41.84%～42.50%。小蕾、中蕾和大蕾中分別得到7881、19 931和10 984個差異表達基因。GO功能分析發(fā)現(xiàn)差異基因在刺激的反應和分子功能上被富集。KEGG通路富集分析表明，次生代謝物的生物合成、代謝途徑、內(nèi)質(zhì)網(wǎng)中的蛋白質(zhì)加工這幾個代謝途徑差異基因較多。對這些基因進行轉(zhuǎn)錄因子分析得知在小蕾、中蕾、大蕾時期較多的差異基因?qū)儆贓RF、MYB、C2H2、NAC、ERKY、bHLH等轉(zhuǎn)錄因子家族；BBR-BPC、GRF、ARR-B、CAMTA、EIL、GeBP等轉(zhuǎn)錄因子家族中的差異基因在少數(shù)。