毛偉兵 陳發(fā)菊 王長(zhǎng)蘭 梁宏偉

(三峽大學(xué)生物技術(shù)研究中心 宜昌 443002)

?

楸樹雄性不育花芽轉(zhuǎn)錄組測(cè)序及分析*

毛偉兵 陳發(fā)菊 王長(zhǎng)蘭 梁宏偉

(三峽大學(xué)生物技術(shù)研究中心 宜昌 443002)

【目的】 以自然突變的楸樹雄性不育花芽為研究材料，分析與楸樹雄性不育相關(guān)基因的表達(dá)模式，從基因的表達(dá)水平上揭示楸樹雄性不育的分子機(jī)制，為楸樹及木本植物的雄性不育研究提供有價(jià)值的參考?！痉椒ā?采用高通量測(cè)序技術(shù)對(duì)自然狀態(tài)下楸樹的雄性不育的花芽以及可育的花芽進(jìn)行轉(zhuǎn)錄組測(cè)序，通過(guò)生物信息學(xué)對(duì)可育花芽與不育花芽的轉(zhuǎn)錄本進(jìn)行比較分析，預(yù)測(cè)和篩選出與楸樹雄性不育有關(guān)的基因?！窘Y(jié)果】 轉(zhuǎn)錄組測(cè)序共產(chǎn)生27.18 Gb數(shù)據(jù)，組裝并去冗余后得到86 076個(gè)Unigene，然后將Unigene比對(duì)到7大功能數(shù)據(jù)庫(kù)(NR, NT, GO, COG, KEGG, Swissprot，Interpro)進(jìn)行功能注釋，最終被7大數(shù)據(jù)庫(kù)中任意一個(gè)數(shù)據(jù)庫(kù)注釋上的Unigene總數(shù)為64 600(75.05%)。將試驗(yàn)組(雄性不育花芽, SL)與對(duì)照組(可育花芽, FL)所測(cè)得的Unigene進(jìn)行表達(dá)量的分析后，篩選出表達(dá)量有差異且可信度高的差異表達(dá)基因。在噪音分布中(差異表達(dá)倍數(shù)在2以上，可信度在0.8以上)共篩選出試驗(yàn)組中有6 915個(gè)基因上調(diào)表達(dá)，3 504個(gè)基因下調(diào)表達(dá)。在泊松分布中(差異表達(dá)倍數(shù)在2以上，錯(cuò)誤發(fā)生率在0.001以下)，SL-1 vs FL-1、SL-2 vs FL-2、SL-3 vs FL-3三個(gè)生物學(xué)重復(fù)中得出差異上調(diào)表達(dá)基因分別有13 979，13 513，13 055個(gè)，差異下調(diào)表達(dá)基因分別有12 170，13 807，10 411個(gè)。差異基因的GO功能分析表明，在生物過(guò)程中顯著性富集的條目集群頻率較高的是生殖過(guò)程、生殖發(fā)育過(guò)程、生殖系統(tǒng)發(fā)育、生殖結(jié)構(gòu)發(fā)育，在分子功能中只有生長(zhǎng)素流出跨膜轉(zhuǎn)運(yùn)蛋白活性顯著性富集。差異基因KEGG通路分析中，差異基因映射到127個(gè)不同的生物途徑，顯著富集代謝通路中差異基因的生物途徑主要為： 代謝途徑、次生代謝產(chǎn)物的生物合成、剪接體、RNA轉(zhuǎn)運(yùn)以及甘油磷脂及淀粉和蔗糖的代謝。在用差異基因與同源基因比對(duì)分析中，共比對(duì)出246個(gè)同源性高的Unigene，其COG功能分類多聚集在RNA加工和修飾，細(xì)胞周期控制、細(xì)胞分裂、染色體分區(qū)，轉(zhuǎn)錄等，同時(shí)將這些同源性高的差異表達(dá)基因映射到了丙酮酸代謝、植物激素信號(hào)轉(zhuǎn)導(dǎo)途徑中?！窘Y(jié)論】 楸樹雄性不育的形成與生殖發(fā)育的多個(gè)過(guò)程、丙酮酸代謝途徑、生長(zhǎng)素流出跨膜轉(zhuǎn)運(yùn)蛋白活性以及油菜素內(nèi)酯介導(dǎo)的信號(hào)轉(zhuǎn)導(dǎo)通路密切相關(guān)。根據(jù)分析的結(jié)果和已經(jīng)完成的相關(guān)細(xì)胞學(xué)觀察，推測(cè)楸樹雄性不育的形成可能是楸樹體內(nèi)丙酮酸代謝過(guò)程異常，抑制油菜素內(nèi)酯的合成，使絨氈層發(fā)育異常并進(jìn)一步影響到小孢子減數(shù)分裂，最終導(dǎo)致不育花粉的形成。

楸樹； 雄性不育； 轉(zhuǎn)錄組測(cè)序； 差異表達(dá)基因

楸樹(Catalpabungei)是紫葳科(Bignoniaceae)梓屬(Catalpa)落葉喬木，是我國(guó)所特有的珍貴優(yōu)質(zhì)用材樹種和觀賞樹種，素有“木王”之稱。其樹干通直，花形若鐘，花朵盛開時(shí)極其優(yōu)美，自古人們就把楸樹作為觀賞樹種廣泛栽培(王新建等， 2004)。自然條件下，楸樹個(gè)體存在雄性不育的現(xiàn)象，小孢子在四分體時(shí)期由于絨氈層結(jié)構(gòu)和功能異常而不能正常發(fā)育，造成花藥敗育而不能形成正常的花粉(張博等， 2015)。楸樹通常采用嫁接和扦插等方式進(jìn)行無(wú)性繁殖，長(zhǎng)期的無(wú)性繁殖導(dǎo)致品種、類型和無(wú)性系單一化，且由于其自花不孕，常常使得楸樹只開花不結(jié)實(shí)，加之對(duì)其過(guò)度的開發(fā)利用，造成楸樹的資源缺乏。

近年來(lái)，轉(zhuǎn)錄組測(cè)序技術(shù)迅速發(fā)展，已被廣泛應(yīng)用到了SSR標(biāo)記規(guī)模化開發(fā)(文亞峰等， 2015； 張振等， 2015)、代謝途徑的分析(陳昊等， 2015； 馬婧等， 2016)、基因挖掘等方面，同時(shí)也在植物雄性不育研究上應(yīng)用并取得進(jìn)展。根據(jù)已報(bào)道的雄性不育轉(zhuǎn)錄組數(shù)據(jù)分析表明，蔥(Alliumfistulosum)細(xì)胞質(zhì)雄性不育與線粒體的氧化磷酸化有關(guān)(Liuetal.， 2016a)，芝麻(Sesamumindicum)雄性不育的形成是在乙烯和茉莉酸介導(dǎo)的信號(hào)通路以及NAC、WRKY 2類轉(zhuǎn)錄因子共同作用下形成的(Liuetal.， 2016b)，而卷心菜(Brassicaoleraceavar.capitata)細(xì)胞核雄性不育則主要涉及ATP合成酶(Guoetal.， 2016)。西瓜(Citrulluslanatus)(Rheeetal.， 2015)、辣椒(Capsicumannuum)(Liuetal.， 2013)等植物也有報(bào)道，但在木本植物中通過(guò)轉(zhuǎn)錄組數(shù)據(jù)分析對(duì)雄性不育相關(guān)基因進(jìn)行研究尚未見報(bào)道。

本研究通過(guò)轉(zhuǎn)錄組測(cè)序技術(shù)，對(duì)楸樹雄性可育與不育花芽進(jìn)行差異表達(dá)基因的分析，在基因轉(zhuǎn)錄水平上對(duì)楸樹雄性不育的機(jī)制進(jìn)行探索，分析楸樹雄性不育發(fā)育過(guò)程中的特定基因的轉(zhuǎn)錄表達(dá)信息，以期為楸樹雄性不育機(jī)制的研究提供理論參考，同時(shí)也為重要性狀相關(guān)基因的克隆以及功能分析等奠定基礎(chǔ)。

1 材料與方法

1.1 試驗(yàn)材料 楸樹花芽采自河南洛陽(yáng)虎頭山基地的楸樹植株，以自然突變的楸樹雄性不育花芽為試驗(yàn)組，記為SL，以雄性可育花芽作為對(duì)照組，記為FL。試驗(yàn)材料采集完成后立即用錫箔紙包好后做好標(biāo)記用液氮處理，保存于-80 ℃冰箱中備用，轉(zhuǎn)錄組測(cè)序樣本每一組各取3個(gè)生物學(xué)重復(fù)，分別為SL-1，SL-2，SL-3和FL-1，F(xiàn)L-2，F(xiàn)L-3。

1.2 RNA提取 采用多糖多酚植物總RNA提取試劑盒(北京天根生化科技有限公司)提取楸樹雄性不育花芽和雄性可育花芽總RNA。利用NanoDrop檢測(cè)RNA濃度，用Agilent 2100檢測(cè)28S/18S以及RIN值，同時(shí)使用1%瓊脂糖凝膠電泳檢測(cè)提取RNA的質(zhì)量以及完整性。

1.3 轉(zhuǎn)錄組建庫(kù)及測(cè)序 提取樣品總RNA并使用DNaseⅠ消化DNA后，用帶有Oligo(dT)的磁珠富集真核生物mRNA； 加入打斷試劑在Thermomixer中適溫將mRNA打斷成短片段，以打斷后的mRNA為模板合成一鏈cDNA，再通過(guò)二鏈合成反應(yīng)體系合成二鏈cDNA，并使用試劑盒純化回收、粘性末端修復(fù)、cDNA的3′末端加上堿基“A”并連接接頭，然后進(jìn)行片段大小選擇，最后進(jìn)行PCR擴(kuò)增； 構(gòu)建好的文庫(kù)用Agilent 2100 Bioanalyzer和ABI StepOnePlus Real-Time PCR System質(zhì)檢，合格后使用Illumina HiSeq 4000或其他平臺(tái)進(jìn)行測(cè)序。

1.4 測(cè)序讀長(zhǎng)分析及拼接 測(cè)序所得數(shù)據(jù)稱為raw reads。首先，過(guò)濾掉低質(zhì)量、接頭污染以及未知堿基N含量過(guò)高的reads，過(guò)濾后的數(shù)據(jù)稱為clean reads。對(duì)過(guò)濾后的數(shù)據(jù)使用Trinity進(jìn)行De novo組裝得到Unigene，接下來(lái)使用Tgicl對(duì)轉(zhuǎn)錄本進(jìn)行聚類去冗余得到Unigene。

1.5 Unigene的功能注釋、分類、代謝途徑分析以及表達(dá)量計(jì)算 使用Blast對(duì)Unigene進(jìn)行NT、NR、COG、KEGG以及SwissProt注釋，使用Blast2GO以及NR注釋結(jié)果進(jìn)行GO注釋，使用InterProScan 5進(jìn)行InterPro注釋。在NT數(shù)據(jù)庫(kù)中比對(duì)分析非冗余核酸序列，利用NR、SwissProt數(shù)據(jù)庫(kù)分析蛋白序列，然后利用COG數(shù)據(jù)庫(kù)根據(jù)直系同源家族蛋白進(jìn)行功能分類以及通過(guò)InterPro數(shù)據(jù)庫(kù)根據(jù)蛋白結(jié)構(gòu)域進(jìn)行家族分類以提供功能注釋。對(duì)所有Unigene基因在分子功能、細(xì)胞組分和生物過(guò)程進(jìn)行GO功能分析以及在KEGG數(shù)據(jù)庫(kù)中進(jìn)行代謝途徑的分析。然后使用RSEM計(jì)算各個(gè)樣品的基因表達(dá)水平。然后使用R軟件里的princomp函數(shù)進(jìn)行PCA分析。

1.6 SL/FL庫(kù)中差異表達(dá)基因分析 根據(jù)需求，使用NOIseq和PossionDis方法進(jìn)行差異基因檢測(cè)。NOIseq(參數(shù)： Fold Change ≥ 2.00； Probability ≥ 0.8)方法基于噪音分布原理，本試驗(yàn)根據(jù)Tarazona等(2011)描述的方法進(jìn)行差異表達(dá)基因(DEG)檢測(cè)。PossionDis(參數(shù)： Fold Change ≥ 2.00； FDR ≥ 0.001)方法基于泊松分布原理，本試驗(yàn)根據(jù)Audic等(1997)中描述的方法進(jìn)行DEG檢測(cè)。將得到的差異表達(dá)基因再進(jìn)行單獨(dú)的聚類分析、GO功能分析、Pathway功能分析，并與已知雄性不育同源基因進(jìn)行比對(duì)分析。

2 結(jié)果與分析

2.1 轉(zhuǎn)錄組測(cè)序以及功能注釋 使用Illumina Hiseq平臺(tái)對(duì)楸樹雄性可育花芽材料FL-1、FL-2、FL-3以及雄性不育材料SL-1、SL-2、SL-3進(jìn)行轉(zhuǎn)錄組測(cè)序，共測(cè)得27.18 Gb原始數(shù)據(jù)。過(guò)濾掉原始數(shù)據(jù)中不合格的reads后，過(guò)濾后的clean reads Q20都達(dá)到了96.84%以上，每組數(shù)據(jù)質(zhì)量見表1。對(duì)過(guò)濾后的數(shù)據(jù)使用Trinity進(jìn)行de novo組裝，聚類去冗余得到Unigene，共得到86 076個(gè)Unigene。將得到的86 076個(gè)Unigene在7大功能數(shù)據(jù)庫(kù)進(jìn)行注釋，注釋結(jié)果分布見表2，最終總共有64 600(75.05%)個(gè)Unigene得到注釋。

表1 過(guò)濾后的reads質(zhì)量統(tǒng)計(jì)Tab.1 The filtered reads quality statistics

表1 過(guò)濾后的reads質(zhì)量統(tǒng)計(jì)Tab.1 The filtered reads quality statistics

樣品Sample原始序列Totalrawreads/Mb有效序列Totalcleanreads/Mb有效堿基數(shù)Totalcleanbases/Gb有效讀數(shù)Q20CleanreadsQ20(%)有效讀數(shù)Q30CleanreadsQ30(%)有效讀數(shù)比例Cleanreadsratio(%)FL?145424529453971293139971FL?245414528453968492569972FL?345424530453975093869975SL?145424530453975294019973SL?245424530453974393779974SL?345364525453972493419975

① FL: 雄性可育花芽(對(duì)照組)； SL： 雄性不育花芽(試驗(yàn)組)。Q20： 質(zhì)量值大于20的堿基數(shù)目占總堿基數(shù)目的比例; Q30： 質(zhì)量值大于30的堿基數(shù)目占總堿基數(shù)目的比例。FL： Male fertile flower buds(control group); SL: Male sterile flower buds(experimental group). Q20: The ratio of the number of bases having a quality value greater than 20 to the total number of bases; Q30:The ratio of the number of bases having a quality value greater than 30 to the total number of bases.

表2 Unigene功能注釋結(jié)果統(tǒng)計(jì)Tab.2 Unigene functional annotation statistical results

表2 Unigene功能注釋結(jié)果統(tǒng)計(jì)Tab.2 Unigene functional annotation statistical results

被注釋unigeneTheannotatedunigene總數(shù)Total數(shù)據(jù)庫(kù)DatabasesNrNtSwissprotKEGGCOGInterproGO總體Overall數(shù)目Number86076604695411242547372542578447203688364600百分比Percentage100%7025%6287%4943%4328%2995%5484%800%7505%

① Nr: Nr庫(kù)屬于非冗余蛋白序列數(shù)據(jù)庫(kù)，數(shù)據(jù)來(lái)源于GenPept, Swissprot, PIR, PDF, PDB以及NCBI RefSeq(ftp://ftp.ncbi.nlm.nih.gov/blast/db); Nt: Nt庫(kù)屬于非冗余核酸序列數(shù)據(jù)庫(kù)，數(shù)據(jù)來(lái)源于GenBank, EMBL以及DDBJ(ftp://ftp.ncbi.nlm.nih.gov/blast/db); SwissProt: UniProtKB/Swiss-Prot利用人工注釋，是一個(gè)高質(zhì)量非冗余的蛋白數(shù)據(jù)庫(kù)(http://ftp.ebi.ac.uk/pub/databases/swissprot); KEGG: KEGG (Kyoto Encyclopedia of Genes and Genomes)用于生物信息研究等，包括基因組、代謝組學(xué)等其他組學(xué)的數(shù)據(jù)分析(http://www.genome.jp/kegg); COG: 直系同源家族蛋白數(shù)據(jù)庫(kù),根據(jù)系統(tǒng)發(fā)生進(jìn)行分類(http://www.ncbi.nlm.nih.gov/COG); InterPro: InterPro根據(jù)蛋白結(jié)構(gòu)域進(jìn)行家族分類以提供功能注釋(http://www.ebi.ac.uk/interpro); GO: Gene Ontology(GO)旨在開發(fā)一個(gè)計(jì)算方法來(lái)描述基因在分子、細(xì)胞和組織水平的功能體現(xiàn)(http://geneontology.org).

Nr: Nr library belongs to non redundant protein sequence databases, data from GenPept, Swissprot, PIR, PDF, PDB and NCBI RefSeq (ftp://ftp.ncbi.nlm.nih.gov/blast/db); Nt: Nt library belongs to non redundant nucleic acid sequence database, data from GenBank, EMBL and DDBJ(ftp://ftp.ncbi.nlm.nih.gov/blast/db); SwissProt: UniProtKB/SwissProt use manual annotation, is a high quality of non redundant protein database (http://ftp.ebi.ac.uk/pub/databases/swissprot); KEGG: KEGG (Kyoto Encyclopedia of Genes and Genomes) is used for the study of biological information, etc., including genome, metabolomics, and other omics data analysis (http://www.genome.jp/kegg); GOG: A family of orthologous protein database, classified according to the system (http://www.ncbi.nlm.nih.gov/COG); InterPro: According to the structure of protein domain family classification to provide functional annotation(http://www.ebi.ac.uk/interpro); GO: The Gene Ontology (GO) attempt to exploit a method to describe Gene in the function of the level of molecules, cells and tissues (http://geneontology.org).

圖1 差異表達(dá)基因統(tǒng)計(jì)結(jié)果Fig.1 The statistical results of differentially expressed genes

2.2 差異表達(dá)基因分析 采用Audic等(1997)的方法對(duì)3組測(cè)序的差異表達(dá)基因進(jìn)行篩選，每組差異表達(dá)基因分布見圖1。差異表達(dá)基因定義為FDR(fail discovery rate)<0.001并且倍數(shù)差異在2倍以上的基因。楸樹雄性不育與雄性可育3組轉(zhuǎn)錄組測(cè)序結(jié)果分析獲得的差異表達(dá)基因數(shù)分別為26 149，27 320，23 466，差異上調(diào)表達(dá)基因分別有13 979，13 513，130 55個(gè)，差異下調(diào)表達(dá)基因分別有12 170，13 807，10 411個(gè)，3組測(cè)序篩選的結(jié)果相似。用噪音分布原理對(duì)楸樹雄性不育與雄性可育進(jìn)行差異表達(dá)基因的篩選，差異表達(dá)基因?yàn)?0 419個(gè)，表明這些差異基因與楸樹雄性不育的關(guān)系密切。在楸樹雄性不育花芽中，有3 504個(gè)差異基因(33.63%)在花芽組織中下調(diào)表達(dá)，有6 915個(gè)差異基因在花芽組織中上調(diào)表達(dá)(66.37%)。

圖2 差異表達(dá)基因GO聚類分析Fig.2 GO clustering analysis of Differentially expressed genes 1-20： 生物過(guò)程； 21-31:細(xì)胞組分； 32-43： 分子功能。 1.生物黏附； 2. 生物調(diào)節(jié)； 3. 組織或生物發(fā)生細(xì)胞組件； 4. 細(xì)胞過(guò)程； 5. 發(fā)育過(guò)程； 6. 生長(zhǎng)； 7. 定位； 8. 細(xì)胞活動(dòng)； 9. 代謝過(guò)程； 10. 多有機(jī)體進(jìn)程； 11. 有機(jī)體多細(xì)胞進(jìn)程； 12. 負(fù)調(diào)控生物過(guò)程； 13. 正調(diào)控生物過(guò)程； 14. 生物過(guò)程調(diào)節(jié)； 15. 生殖； 16. 生殖過(guò)程； 17. 應(yīng)激反應(yīng)； 18. 節(jié)律進(jìn)程； 19. 信號(hào)傳導(dǎo)； 20. 單有機(jī)體進(jìn)程； 21. 細(xì)胞； 22. 細(xì)胞連接； 23. 細(xì)胞成分； 24.細(xì)胞外區(qū)域； 25. 復(fù)雜大分子； 26. 膜； 27. 膜組件； 28. 附膜腔； 29. 細(xì)胞器； 30. 細(xì)胞器組分； 31. 共質(zhì)體； 32. 抗氧化活性； 33. 結(jié)合活性； 34. 催化活性； 35. 電子載體活性； 36. 酶調(diào)節(jié)活性； 37. 分子轉(zhuǎn)導(dǎo)活性； 38. 核酸結(jié)合轉(zhuǎn)錄因子活性； 39. 養(yǎng)分貯液囊活性； 40. 蛋白結(jié)合轉(zhuǎn)錄因子的活性； 41. 受體活性； 42. 分子結(jié)構(gòu)活性； 43. 轉(zhuǎn)錄活性。 1-20： Biological process; 21-31:Cellular component; 32-43： Molecular function. 1.Biological adhesion; 2. Biological regulation; 3. Cellular component organization or biogenesis; 4. Cellular process; 5. Developmental process; 6. Growth; 7. Localization; 8. Locomotion; 9. Metabolic process; 10.Multi-organism process; 11. Multicellular organismal process; 12. Negative regulation of biological process; 13. Positive regulation of biological process; 14. Regulation of biological process; 15. Reproduction; 16. Reproductive process; 17. Response to stimulus; 18. Rhythmic process; 19. Signaling; 20. Single-organism process; 21. Cell; 22. Cell junction; 23. Cell part; 24. Extracellular region; 25. Macromolecular complex; 26. Membrane; 27. Membrane part; 28. Membrane-enclosed lumen; 29. Organelle; 30. Organelle part; 31. Symplast; 32. Antioxidant activity; 33. Binding; 34. Catalytic activity; 35. Electron carrier activity; 36. Enzyme regulator activity; 37. Molecular transducer activity; 38. Nucleic acid binding transcription factor activity; 39. Nutrient reservoir activity; 40. Protein binding transcription factor activity; 41. Receptor activity; 42. Structural molecule activity; 43. Transporter activity.

2.3 差異表達(dá)基因的GO功能分析 對(duì)楸樹雄性不育以及可育的差異表達(dá)基因進(jìn)行GO功能富集，圖2為差異表達(dá)基因GO聚類分析圖。發(fā)現(xiàn)富集在生物過(guò)程中的Unigene有2 329個(gè)，且富集較多的是代謝過(guò)程和細(xì)胞過(guò)程亞類； 富集在細(xì)胞組分中Unigene有1 538個(gè)，富集較多的2個(gè)亞類是細(xì)胞和細(xì)胞部分； 分子功能富集到的Unigene數(shù)為1 361，在該類中結(jié)合作用和催化活性2個(gè)亞類富集最多。顯著富集的(correctedP< 0.05) GO功能亞類9個(gè)(表3)，分子功能中的生長(zhǎng)素外排膜轉(zhuǎn)運(yùn)蛋白活性，生物過(guò)程中的種子休眠期的胚胎發(fā)育、種子發(fā)育、果實(shí)發(fā)育、繁殖器官結(jié)構(gòu)的發(fā)育、繁殖系統(tǒng)的發(fā)育、體軸建成、生殖發(fā)育過(guò)程、繁殖過(guò)程。這些類別在楸樹雄性不育花芽組織中發(fā)揮了重要作用。

2.4 差異表達(dá)基因的KEGG通路分析 以楸樹雄性不育與雄性可育花芽組織的差異表達(dá)基因KEGG pathway聚類分析，差異表達(dá)基因映射到127個(gè)不同的生物途徑，顯著富集的Pathway (Q< 0.05)有42條(表4)，主要為： 代謝途徑、生物合成的次生代謝產(chǎn)物、剪接體、RNA轉(zhuǎn)運(yùn)、內(nèi)吞作用、甘油磷酸酯代謝、淀粉和蔗糖代謝、醚脂類代謝、苯丙素類生物合成等。這些途徑在楸樹雄性不育花芽組織中起重要調(diào)節(jié)作用。

2.5 雄性不育同源基因比對(duì)結(jié)果分析 根據(jù)同源性比對(duì)分析，最終共比對(duì)出246個(gè)Unigene。將得到的Unigene進(jìn)行COG功能分類(圖3)，大多數(shù)都能夠在RNA加工和修飾，細(xì)胞周期控制、細(xì)胞分裂、染色體分區(qū)，轉(zhuǎn)錄中得到注釋，在信號(hào)轉(zhuǎn)導(dǎo)機(jī)制與復(fù)制、重組和修復(fù)中也有Unigene得到注釋。將比較分析后篩選得到的Unigene映射到代謝通路當(dāng)中，最終多集中在丙酮酸代謝通路(圖4)以及與油菜素內(nèi)酯(圖5)相關(guān)的代謝通路當(dāng)中。在丙酮酸的代謝途徑中有丙酮酸轉(zhuǎn)化為草酰乙酸的途徑中(圖4加粗線條標(biāo)示)富集結(jié)果最為顯著； 而在油菜素內(nèi)酯的信號(hào)轉(zhuǎn)導(dǎo)途徑中，合成的油菜素內(nèi)酯在對(duì)油菜素內(nèi)酯不敏感相關(guān)受體激酶及油菜素內(nèi)酯不敏感蛋白復(fù)合物(BRI1)的作用(圖5加粗框標(biāo)示)過(guò)程中富集結(jié)果顯著。因此根據(jù)富集結(jié)果可以分析出丙酮酸代謝途徑和油菜素內(nèi)酯信號(hào)轉(zhuǎn)導(dǎo)途徑對(duì)楸樹雄性不育的形成起著重要的作用。

3 討論

植物雄性不育的性狀是在雄蕊發(fā)育過(guò)程中由一系列極其復(fù)雜的與調(diào)控育性有關(guān)的基因在時(shí)間和空間上共同作用調(diào)控的。引起植物雄性不育主要與花粉形成過(guò)程中的調(diào)節(jié)基因有關(guān)，這些基因在植物花粉發(fā)育的各個(gè)過(guò)程，包括減數(shù)分裂異常(Zhouetal.， 2011； Nonomuraetal.， 2004)、胼胝質(zhì)代謝異常(Wanetal.， 2011)、絨氈層發(fā)育異常(Jungetal.， 2005； Lietal.， 2006)、花粉壁發(fā)育異常(Shietal.， 2011)以及花藥開裂異常(Stintzi， 2000； Steiner-Langeetal.， 2003)等其他過(guò)程中都起著極其重要的作用。

表3 差異表達(dá)基因的GO功能顯著性富集可視化分析Tab.3 Visual analysis of differentially expressed genes which are significant enrichment by GO function

表4 差異表達(dá)基因在KEGG pathway顯著性富集可視化分析Tab.4 Visual analysis of differentially expressed genes which are significant enrichment by KEGG pathway

續(xù)表Continued

SL/FL通路項(xiàng)Pathwayterm差異基因通路注釋DEGswithpathwayannotation(5886)PQ通路位置PathwayID31丙酸代謝Propanoatemetabolism41(070%)389×10-3159×10-2ko0064032RNA轉(zhuǎn)運(yùn)RNAtransport280(476%)405×10-3161×10-2ko0301333苯丙氨酸代謝Phenylalaninemetabolism56(095%)448×10-3173×10-2ko0036034內(nèi)吞作用Endocytosis201(341%)651×10-3239×10-2ko0414435丙氨酸、天冬氨酸和谷氨酸代謝Alanine，aspartateandglutamatemetabolism35(059%)658×10-3239×10-2ko0025036糖胺聚糖降解Glycosaminoglycandegradation24(041%)879×10-3310×10-2ko0053137抗壞血酸代謝Ascorbatemetabolism46(078%)926×10-3318×10-2ko0005338糖酵解和糖異生Glycolysis/Gluconeogenesis95(161%)103×10-2344×10-2ko0001039賴氨酸降解Lysinedegradation32(054%)107×10-2349×10-2ko0031040β-丙氨酸代謝Beta?alaninemetabolism34(058%)141×10-2440×10-2ko0041041甘油磷脂代謝Glycerophospholipidmetabolism189(321%)145×10-2440×10-2ko0056442谷胱甘肽代謝Glutathionemetabolism50(085%)146×10-2440×10-2ko00480

圖3 同源基因COG功能分類Fig.3 COG functional classification of homologous genesA. RNA加工和修飾； D. 細(xì)胞周期控制、細(xì)胞分裂、染色體分區(qū)； K. 轉(zhuǎn)錄； L. 復(fù)制、重組和修復(fù)； R. 一般功能預(yù)測(cè)； T. 信號(hào)傳導(dǎo)機(jī)制。A. RNA processing and modification; D. Cell cycle control, cell division, chromosome partitioning; K. Transcription; L. Replication, recombination and repair; R. General function prediction only; T. Signal transduction mechanisms.

眾所周知，可育花粉的形成與絨氈層的正常發(fā)育密切相關(guān)。絨氈層細(xì)胞發(fā)育異常將直接或間接導(dǎo)致花粉敗育。張博等(2015)發(fā)現(xiàn)，楸樹雄性不育的形成，是由于絨氈層出現(xiàn)高度液泡化，導(dǎo)致其提前解體。本研究結(jié)合差異表達(dá)基因的GO分析和Pathway分析以及同源基因的分析，尋找出與楸樹雄性不育有關(guān)的差異基因顯著富集通路。通過(guò)對(duì)差異表達(dá)基因的GO功能分析，楸樹的生殖過(guò)程、生殖發(fā)育過(guò)程、生殖系統(tǒng)發(fā)育、生殖結(jié)構(gòu)發(fā)育均有顯著性富集，生長(zhǎng)素流出跨膜轉(zhuǎn)運(yùn)蛋白活性也顯著性富集。差異表達(dá)基因的Pathway顯著性分析表明，代謝途徑、次生代謝產(chǎn)物的生物合成、剪接體、RNA轉(zhuǎn)運(yùn)及甘油磷脂、淀粉和蔗糖的代謝富集的差異基因數(shù)目較多。在同源基因的比對(duì)分析基礎(chǔ)上，將與已知相似度高的差異基因再次進(jìn)行功能富集分析。最終顯著性富集結(jié)果主要映射在丙酮酸代謝以及植物激素信號(hào)轉(zhuǎn)導(dǎo)途徑中，這些分析為闡明楸樹雄性不育的調(diào)控途徑提供了線索。

GO分析結(jié)果表明楸樹雄性不育的形成與花發(fā)育過(guò)程中花繁殖器官的形成過(guò)程相關(guān)，才會(huì)形成不育的花粉。這個(gè)結(jié)果與本實(shí)驗(yàn)室在細(xì)胞學(xué)上觀察的楸樹不育花粉的形成是絨氈層的發(fā)育異常相一致。Pathway通路分析表明楸樹雄性不育的形成與生化過(guò)程有著較為緊密的聯(lián)系。在已經(jīng)報(bào)道的雄性不育與絨氈層發(fā)育相關(guān)的基因中，有一類基因?qū)儆赽HLH轉(zhuǎn)錄因子。這一類轉(zhuǎn)錄因子主要參與調(diào)控植物發(fā)育進(jìn)程以及生化過(guò)程(Andradeetal.， 2014； Andriankajaetal.， 2014)； 而且這類轉(zhuǎn)錄因子表達(dá)異常后絨氈層也會(huì)出現(xiàn)高度液泡化(Zhangetal.， 2006； Jungetal.， 2005； Moonetal.， 2013； Liuetal.， 2014)。因此，在對(duì)差異基因Pathway分析結(jié)果做進(jìn)一步的驗(yàn)證時(shí)，應(yīng)將潛在的bHLH轉(zhuǎn)錄因子當(dāng)作進(jìn)一步分析的對(duì)象。

圖4 丙酮酸參與相關(guān)代謝途徑Fig.4 Pyruvate involved in the metabolic pathways粗線條標(biāo)示丙酮酸轉(zhuǎn)化為草酰乙酸的代謝途徑，丙酮酸→磷酸烯醇丙酮酸→草酰乙酸。The thick lines indicate the metabolic pathways of pyruvate conversion to oxaloacetate： Pyruvate→Phosphoenol-pyruvate→Oxaloacetate.

圖5 油菜素內(nèi)酯相關(guān)信號(hào)轉(zhuǎn)導(dǎo)途徑Fig.5 Brassinosteroid related signal transduction pathways加粗框標(biāo)示的BRI1為油菜素內(nèi)酯不敏感蛋白1。BAK1: 油菜素內(nèi)酯不敏感1受體激酶1; BKI1: 油菜素內(nèi)酯不敏感蛋白1激酶抑制劑1； BSK: 油菜素內(nèi)酯信號(hào)激酶； BSU1: 蘇氨酸蛋白磷酸酶激酶； BIN2: 油菜素內(nèi)酯不敏感蛋白2； BZR1/2: 油菜素內(nèi)酯耐受蛋白1/2； TCH4: 低聚木糖葡糖基轉(zhuǎn)移酶； CYCD3: 細(xì)胞周期蛋白D3。The BRI1 marked with the bold box is brassinosteroid insensitive protein 1. BAK1: Brassinosteroid insensitive 1 associated receptor kinase 1; BKI1:Brassinosteroid insensitive protein 1 kinase inhibitor 1; BSK: Brassinosteroid-signaling kinase; BSU1: Threonine-protein phosphatase; BIN2: Brassinosteroid insensitive protein 2; BZR1/2: Brassinosteroid resistant protein1/2; TCH4: Xyloglucosyl transferase; CYCD3: Cyclin D3.

在與雄性不育有關(guān)的同源基因的分析中，得到Unigene集中映射在丙酮酸代謝和油菜素內(nèi)酯參與的信號(hào)轉(zhuǎn)導(dǎo)途徑中。丙酮酸在生物體的3大物質(zhì)代謝過(guò)程中起著極其重要的作用。在丙酮酸的代謝途徑中，它可以通過(guò)乙酰CoA和重要的代謝途徑三羧酸循環(huán)實(shí)現(xiàn)3大物質(zhì)的相互轉(zhuǎn)換。而一旦丙酮酸代謝異常就有可能對(duì)與生殖過(guò)程有關(guān)的包括減數(shù)分裂、絨氈層發(fā)育等途徑造成影響。丙酮酸參與的代謝網(wǎng)絡(luò)途徑復(fù)雜，但由丙酮酸到草酰乙酸的轉(zhuǎn)化這一過(guò)程得到顯著的富集。該過(guò)程也是生物體內(nèi)糖代謝的關(guān)鍵步驟。植物激素轉(zhuǎn)導(dǎo)過(guò)程中，油菜素內(nèi)酯的合成至關(guān)重要。油菜素內(nèi)酯廣泛存在于包括花粉在內(nèi)的許多花器官當(dāng)中。在油菜素內(nèi)酯參與調(diào)節(jié)的生物過(guò)程中，通過(guò)合成的油菜素內(nèi)酯對(duì)油菜素內(nèi)酯不敏感相關(guān)受體激酶和油菜素內(nèi)酯不敏感蛋白復(fù)合物的作用，最終對(duì)DNA的轉(zhuǎn)錄表達(dá)產(chǎn)生影響，使得細(xì)胞伸長(zhǎng)和細(xì)胞分裂出現(xiàn)異常從而影響花粉的正常形成。

結(jié)合本實(shí)驗(yàn)室對(duì)楸樹雄性不育的細(xì)胞學(xué)觀察(張博等， 2015)以及對(duì)轉(zhuǎn)錄組數(shù)據(jù)的分析，對(duì)楸樹雄性不育得出一個(gè)初步的結(jié)果： 在楸樹花藥發(fā)育過(guò)程中，由于丙酮酸參與的代謝途徑出現(xiàn)異常，導(dǎo)致油菜素內(nèi)酯合成途徑受到影響。合成的油菜素內(nèi)酯是通過(guò)作用于細(xì)胞表面由BAK1和BRI1形成的復(fù)合物來(lái)進(jìn)行信號(hào)轉(zhuǎn)導(dǎo)的。油菜素內(nèi)酯合成異常，使得由BAK1和BRI1形成的復(fù)合物活性受到影響，導(dǎo)致參與調(diào)控細(xì)胞伸長(zhǎng)和分裂的基因表達(dá)異常，從而使絨氈層細(xì)胞出現(xiàn)液泡化提前解體，楸樹小孢子在減數(shù)分裂時(shí)期出現(xiàn)異常，最終形成不育的花粉。楸樹的花粉敗育形成的雄性不育，在本試驗(yàn)中已經(jīng)有了初步的結(jié)論，但絨氈層異常形成的雄性不育，其形成過(guò)程十分復(fù)雜，只有通過(guò)更加深入的探索和研究才能對(duì)楸樹雄性不育的機(jī)制有更加清晰的了解。

4 結(jié)論

為了研究楸樹雄性不育形成的機(jī)制，對(duì)1株雄性不育的楸樹花芽和另1株正常發(fā)育的楸樹花芽進(jìn)行轉(zhuǎn)錄組測(cè)序并對(duì)差異表達(dá)基因進(jìn)行分析，結(jié)果顯示顯著性富集的差異表達(dá)基因多在與花的繁殖器官相關(guān)的途徑、重要的代謝途徑和與油菜素內(nèi)酯相關(guān)的信號(hào)轉(zhuǎn)導(dǎo)途徑中。據(jù)此推斷可能是在花粉發(fā)育過(guò)程中，由于重要代謝途徑異常使得油菜素內(nèi)酯的合成途徑出現(xiàn)變化最終導(dǎo)致絨氈層的不正常發(fā)育，使得小孢子減數(shù)分裂異常形成花粉敗育。本研究從轉(zhuǎn)錄水平初步分析了楸樹雄性不育的相關(guān)代謝途徑，獲得的有關(guān)結(jié)果對(duì)楸樹雄性不育的研究有著重要的參考價(jià)值。

陳 昊, 譚曉風(fēng). 2015. 基于油脂合成期油桐種仁轉(zhuǎn)錄組數(shù)據(jù)的α-亞麻酸代謝途徑解析. 林業(yè)科學(xué), 51(3): 41-48.

(Chen H, Tan X F. 2015. Identification of α-linolenic acid metabolism pathway based on transcriptome data ofVerniciafordiikernels during tung oil synthesis stage. Scientia Silvae Sinicae, 51(3): 41-48. [in Chinese])

馬 婧, 鄧 楠, 褚建民, 等. 2016. 泡泡刺高通量轉(zhuǎn)錄組鑒定及其黃酮類代謝途徑初步分析. 林業(yè)科學(xué)研究, 29(1): 61-66.

(Ma J, Deng N, Chu J M,etal. 2016. High-throughput transcriptome identification and flavonoids metabolic pathways inNitrariasphaerocarpa. Forest Research, 29(1): 61-66. [in Chinese])

王新建, 張秋娟, 祝亞軍, 等. 2004. 楸樹新品種及速生豐產(chǎn)技術(shù)研究的現(xiàn)狀與展望. 河南林業(yè)科技, 24(1): 30-31.

(Wang X J, Zhang Q J, Zhu Y J,etal. 2004.Catalpabungeivarieties and fast-growing productive technology research status quo and prospect. Journal of Henan Forestry Science and Technology, 24(1): 30-31. [in Chinese])

文亞峰, 韓文軍, 周 宏, 等. 2015. 杉木轉(zhuǎn)錄組SSR挖掘及EST-SSR標(biāo)記規(guī)模化開發(fā). 林業(yè)科學(xué), 51(11): 40-49.

(Wen Y F, Han W J, Zhou H,etal. 2015. Transcriptome sequencing analysis and development of EST-SSR markers forPinuskoraiensis. Scientia Silvae Sinicae, 51(11): 40-49. [in Chinese])

張 博, 李利平, 毛偉兵, 等. 2015. 雄性不育與可育楸樹花發(fā)育的細(xì)胞學(xué)比較研究. 植物研究, 35(6): 812-818.

(Zhang B, Li L P, Mao W B,etal. 2015. Development of male and female gametophytes between fertile line and male sterile line’s flowers ofCatalpabungeiC. A. Meyer. Bulletin of Botanical Research, 35(6): 812-818. [in Chinese])

張 振, 張含國(guó), 莫 遲, 等. 2015. 紅松轉(zhuǎn)錄組SSR分析及EST-SSR標(biāo)記開發(fā). 林業(yè)科學(xué), 51(8): 114-120.

(Zhang Z, Zhang H G, Mo C,etal. 2015. Transcriptome sequencing analysis and development of EST-SSR markers forPinuskoraiensis. Scientia Silvae Sinicae, 51(8): 114-120. [in Chinese])

Andrade Z I, Baonza A. 2014. The bHLH factors extramacrochaetae and daughterless control cell cycle inDrosophilaimaginal discs through the transcriptional regulation of the Cdc25 phosphatase string. PLoS Genet, 10(3): e1004233.

Andriankaja M E, Danisman S,Mignolet-Spruyt L F,etal. 2014. Transcriptional coordination between leaf cell differentiation and chloroplast development established by TCP20 and the subgroup IbbHLHtranscription factors. Plant Molecular Biology, 85(3): 233-245.

Audic S, Claverie J M. 1997. The significance of digital gene expression profiles. Genome Research, 7(10): 986-995.

Guo J, Zhang Y, Hui M,etal. 2016. Transcriptome sequencing and de novo analysis of a recessive genic male sterile line in cabbage (BrassicaoleraceaL. var.capitata). Molecular Breeding, 36(8): 117.

Jung K H, Han M J, Lee Y S,etal. 2005. RiceUndevelopedTapetum1 is a major regulator of early tapetum development. The Plant Cell, 17(10): 2705-2722.

Li N, Zhang D S, Liu H S,etal. 2006. The ricetapetumdegenerationretardationgene is required for tapetum degradation and anther development. The Plant Cell, 18(11): 2999-3014.

Liu T, Li Y, Zhang C,etal. 2014. Basic helix-loop-helix transcription factorBcbHLHpolfunctions as a positive regulator of pollen development in non-heading Chinese cabbage. Functional and Integrative Genomics, 14(4): 731-739.

Liu C, Ma N, Wang P Y,etal. 2013. Transcriptome sequencing and de novo analysis of a cytoplasmic male sterile line and its near-isogenic restorer line in chili pepper (CapsicumannuumL.). PloS One, 8(6): e65209.

Liu Q, Lan Y, Wen C,etal. 2016a. Transcriptome sequencing analyses between the cytoplasmic male sterile line and its maintainer line in Welsh onion (AlliumfistulosumL.). International Journal of Molecular Sciences, 17(7): 1058. doi: 10.3390/ijms17071058.

Liu H, Tan M, Yu H,etal. 2016b. Comparative transcriptome profiling of the fertile and sterile flower buds of a dominant genic male sterile line in sesame (SesamumindicumL). BMC Plant Biology, 16(1): 250.

Moon J, Skibbe D, Timofejeva L,etal. 2013. Regulation of cell divisions and differentiation byMALESTERILITY32 is required for anther development in maize. The Plant Journal, 76(4): 592-602.

Nonomura K I, Nakano M, Fukuda T,etal. 2004. The novel geneHOMOLOGOUSPAIRINGABERRATIONINRICEMEIOSIS1 of rice encodes a putative coiled-coil protein required for homologous chromosome pairing in meiosis. The Plant Cell, 16(4): 1008-1020.

Rhee S J, Seo M, Jang Y J,etal. 2015. Transcriptome profiling of differentially expressed genes in floral buds and flowers of male sterile and fertile lines in watermelon. BMC Genomics, 16(1): 914. doi: 10.1186/s12864-015-2186-9.

Shi J, Tan H, Yu X H,etal. 2011. Defective pollen wall is required for anther and microspore development in rice and encodes a fatty acyl carrier protein reductase. The Plant Cell, 23(6): 2225-2246.

Steiner-Lange S, Unte U S, Eckstein L,etal. 2003. Disruption ofArabidopsisthalianaMYB26 results in male sterility due to non-dehiscent anthers. The Plant Journal, 34(4): 519-528.

Stintzi A. 2000. TheArabidopsismale-sterile mutant,opr3, lacks the 12-oxophytodienoic acid reductase required for jasmonate synthesis. Proceedings of the National Academy of Sciences, 97(19): 10625-10630.

Tarazona S, García-Alcalde F, Dopazo J,etal. 2011. Differential expression in RNA-seq: a matter of depth. Genome Research, 21(12): 2213-2223.

Wan L, Zha W, Cheng X,etal. 2011. A rice β-1,3-glucanase geneOsg1 is required for callose degradation in pollen development. Planta, 233(2): 309-323.

Zhang W, Sun Y, Timofejeva L,etal. 2006. Regulation ofArabidopsistapetum development and function byDYSFUNCTIONALTAPETUM1 (DYT1) encoding a putative bHLH transcription factor. Development, 133(16): 3085-3095.

Zhou S, Wang Y, Li W,etal. 2011. Pollensemi-sterility1 encodes a kinesin-1-like protein important for male meiosis, anther dehiscence, and fertility in rice. The Plant Cell, 23(1): 111-129.

(責(zé)任編輯 徐 紅)

Transcriptome Sequencing and Analysis of Male SterileFlower Buds inCatalpabungei

Mao Weibing Chen Faju Wang Changlan Liang Hongwei

(BiotechnologyResearchCenter,ChinaThreeGorgesUniversityYichang443002)

【Objective】In order to reveal molecular mechanism ofCatalpabungeimale sterility from the gene expression levels, we analyze the regulation of gene expression patterns forC.bungeimale sterility by studying its natural mutation male sterile flower buds, in order to understand male sterility ofC.bungeiand other woody plants. 【Method】The transcriptome sequencing was carried out on male sterile and fertile flower buds. We use the method of comparative analysis to deal with the sterile and fertile bud transcript through bioinformatics, and predict and filter genes aboutC.bungeimale sterility. 【Result】Transcriptome sequencing produced a total of 27.18 Gb data. Finally, we obtained 86 076 Unigenes after assembling and removing redundant. Then we make use of the seven function databases (NR, NT, GO, COG, KEGG, Swissprot and Interpro) annotating all the Unigenes, a total of 64 600 (75.05%) Unigenes were eventually annotated by any one of the seven databases. Based on the analysis of expression levels of the Unigenes of the experimental group (male sterile flower buds, SL) and control group (male fertile flower buds, FL). Differently expressed genes with different expression levels and high reliability were screened out. In the noise distribution (expression levels differed over 2 times, reliability over 0.8), 6 915 up-regulated genes and 3 504 down-regulated genes were selected from the experimental group. Poisson distribution (expression levels differed over 2 times, incidence of errors under 0.001), in three biological repeats SL-1 vs FL-1, SL-2 vs FL-2, SL-3 vs FL-3 obtained up-regulated genes respectively of 13 979, 13 513 and 13 055 and down-regulated genes respectively of 12 170, 13 807 and 10 411. The GO functional analysis for differently expressed gene showed that the reproduction process, reproduction development process, reproduction system development, reproductive structure development in biological processes were significantly enriched. In molecular function the auxin efflux transmembrane transporter protein activity was significantly enriched. KEGG pathway analysis indicated that differently expressed genes were mapped to 127 different biological pathways. Those significantly enriched pathways mainly contain metabolic pathways, biosynthesis of secondary metabolites, spliceosome, RNA transport and metabolism of glycerophospholipid and starch and sucrose. By comparing differently expressed genes with these genes that are related to male sterility that have been reported, 246 highly homology Unigenes were distinguished. The COG function classification were more focused on RNA processing and modification, cell cycle control, cell division and chromosome partitions, transcription, etc. The differently expressed genes which were highly homology with the male sterility genes were mapped to the pyruvate metabolism, the plant hormone signal transduction pathway. 【Conclusion】 The formation of the male sterility forC.bungeiis involved in multiple processes of reproductive development, pyruvate metabolism pathway, the auxin efflux transmembrane transporter protein activity and brassinosteroid-mediated signal transduction path. Based on the analyses and relevant cytological observations that have been completed, we assume that the male sterility ofC.bungeiis possibly due to abnormal pyruvate metabolism process, leading to abnormal brassinosteroid synthesis and tapetum dysplasia which further affect the meiosis of microspores. These eventually resulted in the formation of sterile pollen.

Catalpabungei； male sterility； transcriptome sequencing； differentially expressed gene

10.11707/j.1001-7488.20170617

2016-07-25；

2016-11-02。

國(guó)家“十二五”科技支撐計(jì)劃課題(2012BAD21B03)； 三峽區(qū)域珍稀植物遺傳發(fā)育與種質(zhì)創(chuàng)新重點(diǎn)實(shí)驗(yàn)室(培育)開放基金項(xiàng)目(2016KBC07)。

S718.46

A

1001-7488(2017)06-0141-10

*梁宏偉為通訊作者。