杜仲雄花芽2個發育時期轉錄組分析

朱利利 杜慶鑫 何 鳳 慶 軍 杜紅巖*

(1.中國林業科學研究院經濟林研究開發中心，鄭州 450003； 2.國家林業局杜仲工程技術研究中心，鄭州 450003)

杜仲(Eucommiaulmoides)為杜仲科(Eucommiaceae)的單科單屬植物，是第四季冰川侵襲后僅留存于我國的孑遺植物,具有重要的經濟價值、藥用價值和生態價值，廣泛分布于我國27個省(市、區)[1]。杜仲為嚴格的雌雄異株植物，雄株所占比例接近50%，雄花結構簡單，為8～12枚針狀雄蕊組成一個雄蕊群簇生長于苞片腋間，產量高、易采集[2～3]。杜仲雄花富含京尼平苷酸、綠原酸、黃酮等多種活性成分，同時含有多種人體必需氨基酸，具有較高的藥用和營養價值[4～6]。目前，杜仲雄花已作為多種飲品原料進行開發，如杜仲雄花芽、杜仲雄花酒、杜仲雄花功能飲料等，具有廣闊的市場開發前景[7]。

植物花芽分化過程是內外因素共同作用的結果，分化數量和質量直接決定花和果實的數量和品質，影響產品開發和商品價值[8]。因此，杜仲雄花芽分化結果特別是雄蕊原基分化對雄花的產量和質量起決定性作用。目前，關于杜仲雄花芽分化過程和時間已進行了觀察，但關于調控杜仲雄花芽分化的基因研究較為匱乏。利用高通量測序技術進行轉錄組學研究，可了解物種在特定環境條件下基因功能表達、生物過程和分子運行機制。本研究前期觀察了杜仲雄花芽分化過程，可分為花序原基分化期、苞片原基分化期、雄蕊原基分化初期、雄蕊形態建成期4個時期，與前人研究結果較一致[9]。根據前期結果，本研究以苞片原基分化期和雄蕊原基分化期的雄花芽為材料，采用Illumina Hiseq X-ten平臺對杜仲雄花芽進行轉錄組測序，并對測序的數據進行過濾，過濾數據以杜仲基因組為參考進行序列比對，利用StringTie進行組裝。根據基因組注釋信息，查詢未被注釋的轉錄區，挖掘新基因，同時根據常用數據庫信息對新基因進行注釋。根據轉錄組信息，了解杜仲雄花芽在苞葉原基發育期和雄蕊原基分化初期的基因表達情況，挖掘調控雄蕊原基發育的基因，為促進雄花發育及分子育種提供重要的數據參考。

1 材料與方法

1.1 試驗材料

試驗材料采集于中國林業科學研究院經濟林研究中心原陽試驗基地(34°55′18″～34°56′27″N，113°46′14″～113°47′35″E)栽培的杜仲良種雄株“華仲11號”(“Huazhong No.11”)。選擇3株長勢良好、無病蟲害、長勢一致的植株，根據杜仲花芽分化結果分別于5月30日(苞葉原基分化中期)，6月14日(雄蕊原基分化初期)采集花芽(圖1)，每株樹采集30個花芽作為一個生物學重復，共3個生物學重復?；ㄑ坎杉罅⒓捶湃胍旱校缓笥?80℃超低溫冰箱保存備用。樣品命名方式：M11-1-1表示苞葉原基分化期雄花芽，M11-2-1表示雄蕊原基分化初期雄花芽，最后一個數字表示生物學重復。

1.2 cDNA文庫構建和轉錄組測序

杜仲雄花芽在液氮中充分研磨后，根據TRIzol試劑盒(Invitrogen)操作步驟提取RNA，其中RNA的濃度、純度和完整性分別用NanoDrop和Agilent 2100檢測。質檢合格的RNA，根據cDNA文庫構建試劑盒(NEB)說明書構建cDNA文庫。利用lllumina Hiseq X-10平臺對cDNA文庫進行雙端測序，測序讀長為150 bp，獲得原始數據。

1.3 數據組裝和新轉錄本預測

根據轉錄組分析標準，測序得到的Raw Data去除質量不合格序列后得到高質量的clean reads，用于后續分析。使用杜仲基因組作為參考和分析，利用軟件HISAT2將clean reads與基因組進行比對[10]。比對分析完成后利用StringTie將比對上的Reads進行組裝，且與杜仲基因組注釋信息進行比較(NCBI登錄號PRJNA357336、SRP095726及SRS2666014)， 查詢未被注釋的轉錄區，挖掘杜仲新轉錄本，補充完善基因組注釋信息[11]。使用BLAST[12]軟件將挖掘的新基因與NR[13]，Swiss-Prot[14]，GO[15]，COG[16]，KOG[17]，Pfam[18]和KEGG[19]數據庫進行序列比對，使用KOBAS2.0[20]得到新轉錄本的KEGG Orthology結果，完成新轉錄本的氨基酸序列預測后使用HMMER[21]軟件與Pfam數據庫比對，獲得新基因的注釋信息。

圖1 杜仲雄花芽2個花芽分化時期形態特征 Br.苞片；SAM.莖尖分生組織；Sq.鱗片；Sta.雄蕊Fig.1 Morphological characteristics of male floral bud at two differentiated stages Br.Bract; SAM.Stem apical meristem; Sq.Squama; Sta.Stamen

1.4 基因表達量定量和差異基因篩選

轉錄本組裝完成采用FPKM(Fragments Per Kilobase of transcript per Million fragments mapped)計算基因表達水平[22]。使用Ballgown對不同分化階段的花芽差異表達基因進行計算和統計[22]。差異基因篩選條件設置為|log2foldchange|>1且校正P-value(FDR)<0.01。

2 結果與分析

2.1 杜仲雄花芽轉錄組測序數據統計

杜仲6個雄花芽轉錄組測序(表1)，原始數據經過濾后共獲得40.48 Gb clean reads,各樣品clean reads均達到6.09 Gb及以上，Q30堿基百分比在90.75%及以上，GC平均含量為46.86%，數據質量較高可用于后續分析。分別將各樣品的clean reads與杜仲基因組進行序列比對，比對效率為91.59%～92.05%。

表1 測序數據統計

圖2 新轉錄本功能注釋結果統計Fig.2 Summary of annotated new transcripts in seven databases

圖3 杜仲雄花芽不同發育時期表達基因和差異表達基因 a.2個發育時期表達基因韋恩圖；b.2個發育時期差異表達基因火山圖Fig.3 Identified and differentially expressed genes at two differentiated stages of E.ulmoides a. Expressed genes at two differentiated stages; b. Differentially expressed genes at two differentiated stages

圖4 差異基因GO功能富集 生物學過程(A.代謝過程；B.細胞過程；C.單個有機體過程；D.生物調節；E.定位；F.刺激應答；G.細胞成分組織或生物發生；H.發育過程；I.信號；J.多細胞生物的過程；K.生殖；L.生殖過程；M.多有機體過程；N.解毒；O.增長；P.免疫系統過程；Q.運動；R.細胞殺傷；S.生物節奏的過程；T.生物附著)；細胞組分(A.細胞；B.細胞部分；C.膜；D.細胞器；E.膜部分；F.細胞器部分；G.大分子復合體；H.胞外區；I.膜密封腔；J.細胞連接；K.超分子復合物；L.細胞外區域部分；M.病毒體；N.病毒部分；O.類核)；分子功能(A.催化活性；B.綁定；C.載體活性；D.結構分子活性；E.核酸結合轉錄因子活性；F.分子功能調節；G.信號傳感器活性；H.電子載體活性；I.抗氧化活性；J.分子傳感器活性；K.營養庫活性；L.轉錄因子活性、蛋白質綁定；M.金屬伴侶活性；N.蛋白質標記)Fig.4 GO analyses of differentially expressed genesBiological process(A.Metabolic process; B.Cellular process; C.Single-organism process; D.Biological regulation; E.Localization; F.Response to stimulus; G.Cellular component organization or biogenesis; H.Developmental process; I.Signaling; J.Multicellular organismal process; K.Reproduction; L.Reproductive process; M.Multi-organism process; N.Detoxification; O.Growth; P.Immune system process; Q.Locomotion; R.Cell killing; S.Rhythmic process; T.Biological adhesion); Cellular component(A.Cell; B.Cell part; C.Membrane; D.Organelle; E.Membrane part; F.Organelle part; G.Macromolecular complex; H.Extracellular region; I.Membrane-enclosed lumen; J.Cell junction; K.Supramolecular complex; L.Extracellular region part; M.Vvirion; N.Virion part; O.Nucleoid; P.Catalytic activity); Molecular function(A.Catalytic activity; B.Binding; C.Transporter activity; D.Structural molecule activity; E.Nucleic acid binding transcription factor activity; F.Molecular function regulator; G.Signal transducer activity; H.Electron carrier activity; I.Antioxidant activity; J.Molecular transducer activity; K.Nutrient reservoir activity; L.Transcription factor activity,protein binding; M.Metallochaperone activity; N.Protein tag)

2.2 新轉錄本發掘和注釋

基于clean reads在基因組比對結果，進行新轉錄本的發掘，共發掘新轉錄本2 929個,將發掘的新轉錄本與NR，Swiss-Prot，GO，COG，KOG，Pfam，KEGG 7大數據庫進行序列比對，共2 166個新轉錄本得到功能注釋，占總得到新轉錄本的73.95%(圖2)。其中，注釋到NR蛋白數據庫中的新轉本數量最多為2 136個占總新轉錄本數72.93%，然后依次注釋到Swiss-Prot數據庫數量為1 509個占51.52%、GO數據庫1 336個占45.61%，Pfam數據庫1 235個占42.16%，KOG數據庫為1 174個占40.08%，KEGG數據庫865個占29.53%，而在COG數據庫注釋數量最少為571個占19.49%。

2.3 杜仲雄花芽2個發育時期差異基因表達分析

杜仲雄花芽2個發育時期共鑒定獲得20 226個表達基因，2個發育時期共表達基因18 094個，苞葉基分化期特異表達778個，雄蕊原基分化初期特異表達1 354個(圖3a)。結果顯示約89.5%基因在杜仲雄花芽2個發育階段均表達，僅約3.8%和6.7%基因在苞葉原基和雄蕊原基發育階段特定表達。根據差異基因篩選條件，在2個發育時期篩選出583個差異表達基因，其中在雄蕊原基發育初期上調基因315個，下調基因267個(圖3b)，可能是這些特異表達基因特別是差異表達基因促進了杜仲雄花芽的雄蕊原基萌發。

2.4 杜仲雄花芽差異基因GO功能分析

為了進一步挖掘參與杜仲雄蕊發育的基因，對雄花芽苞葉原基發育時期和雄蕊原基發育初期的差異基因進行GO功能富集分析，篩選出雄花芽中與成花誘導緊密相關的GO條目(圖4)。雄花芽在苞葉原基發育期和雌蕊原基發育初期之間共有359個差異基因富集到個1 388個GO條目(GO terms)，其中顯著性富集GO條目(P-value<0.05)有91個，其中有52個GO條目顯著富集到生物學過程，11個顯著富集到細胞組分，28個顯著性富集到分子功能。雄花芽顯著性富集的GO條目包括與分生組織發育相關的過程如分生組織從營養生長向生殖生長的轉變(vegetative to reproductive phase transition of meristem)、生長調節(regulation of growth)、生長負向調控(negative regulation of growth)，與激素相關的激素生物合成過程如吲哚硫代葡萄糖苷代謝過程(indole glucosinolate metabolic process)、乙烯代謝過程(ethylene metabolic process)和乙烯生物合成過程(ethylene biosynthetic process)，與光周期相關的遠紅光反應(response to far red light)，光周期和開花(photoperiodism,flowering)，光周期現象(photoperiodism)。同時差異基因被富集到核酸結合轉錄因子活性(nucleic acid binding transcription factor activity)、發育過程(developmental process)，生殖(reproduction)和生殖過程(reproductive process)等與開花緊密相關的生物過程。說明杜仲雄花芽中參與營養生長向生殖生長的轉變、激素合成與代謝、光周期過程、核酸結合轉錄因子活性、生殖過程等基因的差異表達可能是促進杜仲雄花芽由苞葉原基發育向雌、雄蕊原基發育轉變的直接促進因子。

2.5 杜仲雄花芽差異基因KEGG功能分析

杜仲雄花芽2個發育時期有98個差異表達基因被富集到70條KEGG pathway。其中富集比例前5的pathway分別是苯丙素生物合成富集比例11.22%，其次是淀粉和蔗糖代謝富集比例11.22%，碳代謝富集比例10.20%，內質網蛋白質加工富集比例10.20%，氰基氨基酸代謝富集比例8.16%。同時差異基因被富集到植物激素信號轉導富集比例為3.06%，以及顯著富集(P-value<0.05)到植物生物節律富集比例6.12%(表2)。代謝通路富集結果說明杜仲雄花芽在雄蕊發育過程中需要碳水化合物作為基礎，且需要激素和一些次生代謝物調控，同時杜仲雄花芽形態分化與植物節律密切相關。

表2 差異表達基因代謝通路注釋

***P<0.001

2.6 開花相關差異基因分析結果

根據GO和KEGG富集分析，營養生長向生殖生長的轉變、激素合成與代謝、激素信號轉導、光周期過程和植物生物節律，核酸結合轉錄因子活性等生物過程與杜仲雄花芽發育密切相關，對調節這些生物途徑的差異基因進行詳細分析。分析結果顯示，與激素相關的差異表達基因有16個，其中與生長素(IAA)相關的有6個差異基因，有4個上調差異基因分別是ARF9(EUC06536-RA)、GH3.10(EUC06844-RA)、ABCB19(EUC13060-RA)、PILS5(EUC02516-RA)和下調基因PCNT115(EUC23432-RA和EUC25737-RA)；與乙烯(ETH)代謝和傳導途徑相關的有9個差異基因，包含上調因子RAP2-7(EUC00521-RA)、AFP3(EUC_newGene_23655)、AHK1(EUC08386-RA)、EREBP-like(EUC00228-RA)，FER(EUC20697-RA)、AFP2(EUC15307-RA)、ACO(EUC22557-RA)、CTR1(EUC03132-RA)和一個下調因子ERF13(EUC_newGene_3252)；而與脫落酸(ABA)相關的僅有下調因子PP2C基因(EUC22646-RA)。

杜仲雄花芽2個發育時期發現了調控光周期開花途徑重要轉錄因子COL差異表達，包括上調轉錄因子COL4(EUC12564-RA)和COL2(EUC15205-RA)，下調因子COL9(EUC17024-RA)，同時還發現與生物節律相關2個上調基因CCA1(EUC05190-RA)和光敏色素PHYA(EUC15838-RA)以及五個下調基因PRR73(EUC00514-RA)，APRR5(EUC16182-RA)，ADO3(EUC22871-RA)，GI(EUC19998-RA)和ELF4(EUC18826-RA)。此外，結果顯示有多種轉錄因子家族成員差異表達調控杜仲雄花芽從營養生長向生殖生長轉變，包括MADS-box中MIKC亞家族基因上調轉錄因子AGL3(EUC07615-RA)和AGL8(EUC07648-RA)，以及開花整合子FLC(EUC06393-RA)和SOC1(EUC03328-RA)；WRKY家族上調轉錄因子WRKY41(EUC07015-RA、EUC22632-RA和EUC07417-RA)、WRKY24(EUC00375-RA)；MYB家族上調轉錄因子MYB61(EUC06946-RA)、CCA1(EUC05190-RA)和RVE8(EUC03052-RA)，以及下調轉錄因子LUX(EUC07116-RA和EUC11242-RA)；GATA家族上調因子GATA18(EUC11132-RA)和下調因子GATA9(EUC18806-RA)；同時還包括bHLH、TCP、PIF1以及新轉錄因子家族PLATZ參與雄花芽的雄蕊發育過程。在本研究中，杜仲雄花芽轉錄組差異分析結果還顯示與植物分生組織相關的基因AGO2和AGO7在雄花芽苞葉原基分化期和雄蕊原基分化初期差異表達。

3 討論

本研究以杜仲雄花良種“華仲11號”為材料，對處于苞葉原基發育時期和雄蕊原基發育時期的花芽進行轉錄組測序，功能富集結果顯示差異表達基因顯著富集在淀粉和蔗糖代謝途徑和碳代謝途徑。研究表明糖類在成花過程中可作為能量物質和結構物質促進花芽分化，而且可以作為成花信號物質直接作用于莖尖調控花芽分化[23]。植物成花過程中除了需要以碳水化合物作為基礎外，而且與激素調控密切相關。大量研究證明激素參與植物的開花過程，對植物花芽分化有顯著調控作用。生長素極性運輸對花的形成具有重要作用，其中擬南芥(Arabidopsisthaliana)pin-1突變體的花序軸的內源生長素和生長素的極性運輸能力比起野生型明顯降低，其花芽不能正常發育成完整的花[24]。本文結果顯示有15個差異基因與生長素、乙烯和ABA的生物代謝和傳導途徑相關，其中與乙烯和生長素相關的基因除ERF13外其他都在雄蕊原基發育階段表現為上調，而與ABA相關的基因則表現為下調。其中與生長素響應元件結合的生長素響應因子(ARF9)，具有維持植物生長素動態平衡的GH3基因，生長素轉運蛋白基因ABCB19和PILS5在杜仲雄蕊原基發育時期表達上調，而生長素誘導蛋白基因PCNT115則表達下調。乙烯具有刺激開花及花的衰老，其中參與乙烯信號傳導的RAP2，AFP3，CTR1等是乙烯信號途徑中一個重要的調節因子，在杜仲雄蕊原基發育時期表達上調。乙烯響應因子ERF在玉米中可以調控胚乳發育[25]，本研究顯示ERF在杜仲雄蕊原基發育時期表達下調。對ABA調控花芽分化研究顯示，低含量的ABA有利于板栗(Castaneamollissima)花簇原基的分化和生長[26]。ABA受體PYL蛋白是ABA信號轉導通路中重要成員，可與ABA結合抑制PP2C表達[23]；在本研究中，PP2C在杜仲雄花芽苞葉原基時期表達上調，而在雄蕊原基發育期可能受到抑制表達下調。綜合以上結果，生長素、乙烯可能正調控杜仲雄蕊原基的形態發育，而ABA可能負調控杜仲雄蕊原基發育。

光周期途徑是植物成花誘導的重要途徑，植物在感受適宜的光周期后，通過生物鐘調控開花基因表達。本研究顯示與光周期和物種生物鐘密切相關的基因COL成員(COL4、COL2、COL9)、CCA1、PHYA、PRR73、APRR5、ADO3、GI和ELF4在杜仲雄花芽發育不同時期差異表達。其中，CO同源基因COL是光周期途徑成花的關鍵基因，擬南芥ATCOL2表達受光照調控，在黎明時表達最高，傍晚時表達急劇下降，但在轉基因植物中其表達量對開花時間無太大影響[27]；而從小麥(Triticumaestivum)中克隆的TACOL4轉化到水稻(Oryzasativa)中可延遲水稻抽穗[28]；在玉米(Zeamays)中COL9可通過調控Ehd1途徑延長玉米開花時間[29]。GI是擬南芥核心生物鐘基因之一，調控下游CO基因的表達來影響植物開花時間和進程[30]。光敏色素基因PHYA調控CO表達，生物鐘蛋白基因CCA1和LHY通過降低SVP蛋白累積加速擬南芥開花進程，此外早花基因ELF4、PRR73、APRR5等基因也是調控生物鐘和開花過程的重要因子[31～32]。因此，光周期途徑很可能是誘導杜仲成花的主要途徑。

轉錄因子家族在植物花器官發育、果實發育以及逆境脅迫等方面發揮重要作用，其中MADS-box基因是重要的花器官屬性決定基因。MADS-box家族成員FLC、SOC1、AGL24和SVP是植物成花誘導途徑春化途徑、自主途徑、赤霉素途徑和光周期途徑中的關鍵節點作用基因[33]。本文研究發現了4個TypeⅡMADS-box型MIKC亞家族成員在杜仲雄花芽2個發育時期差異表達，其中AGL3和AGL8在苞葉原基發育期表達下調而在雄蕊原基發育初期表達上調，而SOC1、FLC與AGL基因表達正好相反。FLC是一個開花抑制因子，春化途徑與自主途徑中發揮作用，在擬南芥中過表達FLC則可推遲開花[34]。開花整合子SOC1受FLC的負調控抑制植物成花，同時也受CO信號在長日照條件下誘導植物成花，而且也能接受赤霉素的調控誘導植物開花[35～36]。在ABC(DE)模型中，SEP類基因被命名為E類，而AGL3基因為第四個SEP基因，影響著四輪花器官和胚珠發育[37]。此外，在擬南芥中花原基中發現了MADS-box基因AGL8，在營養生長階段不表達，而在花序頂端分生組織、花梗和莖生葉中高表達[38]。所以，在杜仲雄花芽發育過程中SOC1和FLC在苞葉原基發育時期高表達調控花分生組織發育，然后SEP類基因被激活進而調控杜仲雄蕊原基發育，但這4個基因在杜仲雄蕊發育過程中的具體功能還需進一步實驗驗證。

4 結論

杜仲雄花芽2個發育時期獲得583個差異表達基因，雄蕊原基發育初期上調基因315個，下調基因267個。差異基因的GO和KEGG富集結果顯示，差異基因富集在營養生長向生殖生長的轉變、乙烯激素合成與代謝，核酸結合轉錄因子活性，光周期現象等生物過程以及碳水化合物代謝、植物激素信號轉導和植物生物節律等在花芽分化過程中起重要調控作用的代謝通路。根據差異基因分析結果，發現在杜仲雄花芽形態分化過程中除需要以碳水化合物作為基礎外，還受生長素、乙烯和脫落酸等激素調控，同時光周期途徑可能是杜仲成花誘導的重要途徑。此外，MADS-box家族成員FLC、SOC1、AGL3和AGL8與杜仲雄蕊器官發育密切相關。本研究結果從轉錄水平上分析了杜仲雄花芽發育的相關代謝途徑和基因，可在一定程度上解析杜仲雄花芽形態分化的分子調控模式與機制。