基于轉錄組測序的高羊茅分蘗與株高相關差異表達基因分析

楊志民，邢瑞，丁鋆嘉，莊黎麗

（南京農業大學草業學院，江蘇 南京 210095）

分蘗和株高是禾本科牧草和草坪草的兩個重要表型特征。高植株多分蘗的牧草品種地上生物量大，可以獲得高產；矮植株多分蘗的草坪草品種能夠快速成坪，成坪后草坪質量高，耐踐踏，少修剪，可節約養護成本。此外，在草坪草遭受諸如干旱、鹽堿和熱激等非生物脅迫后，分蘗對于其恢復再生也十分重要。植物分蘗發育過程受遺傳和環境信號的調控，分蘗的發育一般分為兩步——葉腋分生組織（axillary meristems，AMs）和分蘗芽的起始以及分蘗芽的伸長［1—2］。葉腋分生組織位于植物的葉腋處，在植物的整個生命周期內，AMs會生長出分枝或是分蘗。在模式植物腋芽缺失突變體的試驗中，一些基因已經被證實是調控AMs發育的關鍵基因［3］，例如，番茄（Solanum lycopersicum）的Gibberellic acid insensitive（GAI），Repressor of GA 1-3 mutant，Scarecrow（GRAS）同源轉 錄 因 子LATERAL SUPPRESSOR（LS）［4］、水 稻（Oryza sativa）的MONOCULM 1（MOC1）［5］和 擬 南 芥（Arabidopsis thaliana）的LATERAL SUPPRESSOR（LAS）［6］。MOC1功能缺失會導致分蘗數和生殖枝數下降［5］。miR164負調控LAS，可直接促進CUP SHAPED COTYLEDON1（CUC1）和CUC2轉錄因子的降解［7］。REGULATOR OF AXILLARY MERISTEMS1（RAX1）是另一種重要的AMs起始因子，可直接正向調控CUC2在邊界區域的特殊表達。RAX1還可以和RAX2、RAX3協同調控AMs的起始［8］。TEOSINTE BRANCHED1，CYCLOIDEA，PCF（TCP）轉錄因 子家族的FINE CULM 1/PsBRANCHED1/BRANCHED 1（FC1/PsBRC1/BRC1）等基因在腋芽處特異表達，通過整合各種通路如激素響應和光形態建成等負調控植株分枝［9—11］。目前對調控腋芽起始和伸長的分子機制主要集中在模式植物，其他植物中側枝/分蘗發育分子機制的研究還有很多欠缺。

植株高度受遺傳和內源激素等多種因素調控。赤霉素（gibberellin，GA）可以打破植物休眠、促進種子萌發、莖節伸長和調節莖稈形態。在西瓜（Citrullus lanatus）中過表達ClTCP14a和ClTCP15可以提高內源GA含量，促使植株增高［12］。過表達GA 2oxs會使植物內源活性GA含量下降，使植株矮化，例如，在菠菜（Spinacia oleracea）中過表達SoGA 2ox3［13］、在擬南芥中過表達AtGA 2ox7和AtGA 2ox8［14］和在矮牽牛（Petunia hybrid a）中過表達PhGA 2ox1［15］，植株均出現矮化表型。在GA信號轉導通路中，GID1作為GA受體可與GA和DELLA蛋白結合形成復合體。DELLA蛋白是GA信號傳導途徑的阻遏因子，當GA濃度較高時，DELLA蛋白會被蛋白酶體降解失活，植物出現正常的GA響應反應［16］。吲哚乙酸（indole-3-acetic acid，IAA）是最先發現的植物激素，參與胚胎形成、器官分化、莖與根的伸長生長等發育過程，具有頂端優勢效應，可抑制側芽的生長。IAA可以和油菜素甾醇（brassinosteroids，BR）互作產生累加的生物學效應，可以促進細胞的伸長與分化［17］。并且IAA可以誘導DWARF4的表達，促進BR的生物學合成［18］。BR的一個重要功能是促進細胞的伸長與分裂進而促進植株地上部的生長發育。擬南芥BR缺失突變體［19］與水稻BR缺失突變體均呈現矮化表型，并且水稻在外源施用BR可以恢復正常表型［20］。

高羊茅（Festuca arundinacea）是應用最廣泛的冷季型禾草，既可用于建植草坪，又是優良牧草。“Kentucky-31”（K31）是典型的高羊茅牧草品種，植株高大，但分蘗能力較弱；“Regenerate”是草坪型培育品種，優點是植株相對較矮且分蘗能力強。兩個品種明顯不同的表型是由哪些基因的差異表達造成的目前尚不清楚。高羊茅的基因組較大（約6 GB），當前已有一些表達數據標簽（expressed sequence tag，EST）及基于非生物脅迫如抗熱、抗重金屬等研究所獲得的轉錄組序列［21—22］，但有關高羊茅分蘗和株高方面的轉錄組序列研究，目前尚未有報道。下一代測序技術和RNA測序可以分析轉錄組序列，挖掘差異表達基因（differentially expressed genes，DEGs），甚至可以檢測尚無基因組參考的物種的基因序列。本研究目的在于獲取更多的高羊茅轉錄組序列信息，并且比較“K31”和“Regenerate”在生長發育尤其是分蘗和株高發育上的分子機制差異，以便進一步挖掘調控高羊茅腋芽發育和植株高度的基因。

1 材料與方法

1.1 植物材料、培養條件及處理

本研究的試驗材料為高羊茅品種“K31”和“Regenerate”，由美國羅格斯大學植物生理課題組提供。試驗于2015年2月在南京農業大學草坪科學實驗室進行。將兩種高羊茅種子統一置于濕潤的濾紙上，并在黑暗條件下低溫春化。春化后，將種子置于光下萌發。取15株7 d齡的幼苗移栽到直徑20 cm的塑料盆中，培養基質采用陶粒土。所有植物材料都栽培于環境可控的人工氣候室內（MT 8070iE，河南旭邦），光周期設置為白天14 h/黑夜10 h，空氣濕度控制在70%，光密度為300μmol·m—2·s—1，晝/夜溫度為25oC/15oC，每周補充一次Hoagland營養液。

為了獲得高羊茅de novo拼接轉錄組序列，分別采集“K 31”和“Regenerate”各類植物樣本（各品種1個生物學重復），包括：植物全株（2周齡）、葉片、根、分蘗節、腋芽、葉耳和莖（3月齡）。將所有樣本混合并儲存在—80oC冰箱用于提取RNA。

為了獲得“K 31”和“Regenerate”分蘗節中差異表達基因，在3月齡植株分蘗節部位取樣（各品種3個生物學重復），將樣品液氮速凍后于—80oC冰箱儲存。

1.2 總RNA提取、cDNA文庫構建與轉錄組測序

材料總RNA的提取使用RNeasy Plant Mini Kit（Qiagen，美國）試劑盒。采用NanoDrop（Nano-Drop Technologies，美國）和Agilent 2100 Bioanalyzer（Agilent Technologies，美國）檢測樣本純度和濃度，采用瓊脂糖凝膠電泳檢測RNA完整性，以保證進行轉錄組測序的樣品符合標準。

采用Illumina TruSeq RNA Sample Prep Kit（Illumina Inc.，美國）和Ultra RNA Library Prep Kit for Illumina構建cDNA文庫，再用Beckman AMPure XP beads（Beckman Coulter，美國）純化PCR產物。使用Agilent 2100對文庫大小進行檢測，保證測序文庫質量。文庫構建完成后用Bio-RAD CFX 96進行熒光定量PCR（qPCR）準確定量。采用HiSeq 2500結合Miseq 300測序平臺（Illumina Inc.，美國）進行“Regenerate”和“K31”轉錄組深度測序用于拼接組裝，采用HiSeq 2000測序平臺對“Regenerate”和“K31”分蘗節進行差異表達基因分析。

1.3 高羊茅轉錄組序列組裝、注釋及分析

使用FASTX-Toolkit軟件，將測序得到的原始數據進行過濾，去除接頭和低質量reads（Q＜30，length＜50 bp）。然后通過Trinity軟件對clean reads進行組裝獲得轉錄本序列，利用CAP3軟件得到一套最終的單基因簇（unigenes）。利用BlastX軟件將unigenes的序列比對到NCBI非冗余蛋白數據庫（non-redundant protein database，NR）、蛋白質數據庫（universal protein，Uniprot）、東京基因與基因組數據庫（kyoto encyclopedia of genes and genomes，KEGG）和蛋白質直系同源數據庫（cluster of orthologous groups，COG）中，得到unigenes的蛋白功能注釋信息，然后利用WEGO軟件進行unigenes的基因本體數據庫（gene ontology，GO）功能分類統計。

1.4 DEGs篩選注釋及Real-time PCR驗證

取兩個高羊茅品種的分蘗節部位轉錄本進行差異表達分析，以“K31”為對照組，“Regenerate”為處理組，應用RPKM（reads per kilobase per million reads）方法分析基因表達水平，以差異倍數（fold change，FC）＞1或＜—1且P值≤0.05為篩選標準，之后對獲得的DEGs進行GO功能分類統計。在DEGs中挑選出與高羊茅生長發育相關的關鍵基因進行熒光定量PCR驗證。采用2—ΔΔCT法分析基因的相對表達量，以α-Tubblin（Accession No.GT 051159）（F：ATGCTTTCGTCTTATGCCC；R：CTCTTGGTTTTGATGGTTGC）作 為 內 參 基 因，用Primer Premier Software（Version 5.0）設計定量PCR引物序列（表1），每個樣本3次生物學重復，每個生物學重復做3次技術重復。

表1 PCR驗證引物Table 1 PCR validation primers

1.5 植物表型觀測

每個品種設計3個生物學重復，每個重復內有20株植物材料。從萌發后4周開始分別統計株高直到第7周結束，每周測量一次。將栽培基質到葉尖的最大長度定為株高。

統計種子萌發后4～7周的分蘗數和腋芽數，在體視顯微鏡（Olympus，日本）下觀察腋芽并拍照。

1.6 數據分析

采用一般線性統計模型（SAS 9.0，美國）對生長和生理數據進行方差分析（ANOVA）。采用Fisher最小顯著性檢驗（LSD，P＜0.05）進行均數分離。

2 結果與分析

2.1 兩個高羊茅品種的表型差異

本研究中，“K31”和“Regenerate”在生長發育方面有不同的表型特征（圖1A，B）。從種子萌發后4周開始統計分蘗數和腋芽數、株高、莖稈直徑、葉細胞大小等生長數據。萌發后第4周及第5周時，“K31”和“Regenerate”的分蘗數沒有顯著差異，但是在第6周和第7周時差異顯著（圖1C）。在株高方面，萌發后第4周，“K 31”植株高度為24.53 cm，顯著高于“Regenerate”的16.1 cm。萌發后第7周，“K 31”已經長至50.78 cm，而“Regenerate”為28.01cm（圖1D）。兩者莖葉（此處莖指的是層層包裹的葉鞘）生長速率也有明顯不同，“K31”的莖比“Regenerate”更為粗壯，萌發后第4周時“K31”莖稈直徑為1.33 mm，“Regenerate”直徑為1.22 mm。萌發后第6周時，“K31”的直徑為2.35 mm，“Regenerate”為1.91 mm（圖1E）。表皮細胞面積同樣有顯著差異，萌發后第4周時，“K 31”細胞面積為8887.29μm2，“Regenerate”細胞面積為6536.32μm2；第7周時，細胞面積縮小，“K31”細胞面積為6674.48μm2，“Regenerate”細胞面積為5629.03μm2（圖1F）。

圖1 “Regenerate”和“K 31”植株的表型差異比較Fig.1 Compar ison of phenotypic differ ence between“K 31”and“Regenerate”plants

2.2 組裝與測試結果分析

2.2.1高羊茅測序結果和組裝 “K31”和“Regenerate”兩個品種高羊茅組成的RNA-seq文庫通過Hiseq 2500結合Miseq 300測序平臺進行測序。“K31”獲得76673214條raw reads，“Regenerate”獲得88475512條raw reads（表2），原始測序數據已提交到NCBI的序列片段歸檔（sequence read archive，SRA）數據庫（Bioproject ID：PRJNA 669484）。將原始序列數據中的接頭污染數據、低質量數據以及重復冗余數據去除后，對“K31”的76533900條clean reads和“Regenerate”的88345910條clean reads進行后續分析。Clean reads分別占比97.60%和97.81%（表2），證明所得到的數據是優質數據。“K31”和“Regenerate”的clean reads共獲得77872條unigenes，unigenes長度范圍從200 bp跨越到13402 bp，平均長度為611.086 bp（N50＝920 bp，N90＝262 bp）（表3）。Unigenes的GC含量為53.2%（表3），接近水稻的GC含量（55%）。

表2 de novo RNA-seq詳細序列信息Table 2 Detail sequence information for de novo RNA-seq

表3 樣品RNA-Seq數據匯總Table 3 Summary of sample based on the RNA-Seq data

2.2.2高羊茅轉錄組unigenes的物種分布與COG功能注釋 利用BlastX將組裝得到的77872條unigenes分別與NR、Uniprot/Swiss-Prot、KEGG、COG和GO數據庫（E-value＜1e—5）進行比對的結果顯示：共有60552條unigenes獲得蛋白功能注釋，占總unigenes的77.76%。其中有8206條unigenes在所有數據庫都有獲得功能注釋（圖2，表4）。從注釋匹配的物種分布來看，來自二穗短柄草（Brachypodium distachyon）的注釋最多，15882條注釋，占總注釋的26.50%；其次為大麥（Hordeum vulgure），13594條注釋，占總注釋的22.68%。說明在已知基因組信息的物種中，高羊茅與二穗短柄草和大麥的序列相似度較大。其他與高羊茅序列相似度較大的物種依次為山羊草（Aegilops tauschii）（17.58%）、烏拉爾圖小麥（Triticum urartu）（11.87%）、粳稻（Oryza sativaJaponica Group）（4.53%）和玉米（Zea mays）（4.50%）等。除此之外，有388條unigenes（0.65%）沒有比對到任何物種（圖3）。高羊茅unigenes與COG數據庫對比結果表明：共有13767條unigenes被歸納到24個不同功能分類當中。其中，有3801條unigenes為一般功能預測蛋白（general function prediction only），占總功能序列的4.88%，在所有功能基因中占比例最高；然后依次是蛋白質翻譯后修飾—蛋白質周轉—分子伴侶（posttranslational modification，protein turnover，chaperones，1970 unigenes，2.53%）、轉錄（transcription，1766 unigenes，2.27%）和信號轉導機制（signal transduction mechanisms，1628 unigenes，2.09%）。只有少數序列歸納到細胞能動 性（cell motility，43 unigenes，0.06%）和 核 結 構（nuclear structure，6 unigenes，0.01%）（圖4）。

圖3 物種分布結果非冗余注釋Fig.3 Species distr ibution of the results of non-r edundant（NR）annotation

圖4 高羊茅unigenes COG注釋分類Fig.4 Cluster of orthologous group（COG）classif ication of unigenes in F.arundinacea

表4 數據庫BLAST分析結果Table 4 BLAST analysis results against important public databases

圖2 獲得注釋的unigenes在5個數據庫的分布Fig.2 Distribution of unigenes annotated in five databases

2.2.3高羊茅轉錄組unigenes的GO注釋及分類高羊茅轉錄組unigenes的GO功能注釋所分析的結果顯示：共有44447條unigenes（57.08%）獲得功能注釋，分布在48個亞類中。其中，有28316、37017和31216個unigenes分別注釋到細胞組分、分子功能和生物過程（圖5）。在細胞組分的12個功能亞類中，其中以細胞（cell）和細胞組分（cell part）較多，二者均占該類別的36.1%。其次為細胞器（organelle）和大分子復合物（macromolecular complex），分別占該類別的25.6%和7.0%。在分子功能所包括的13個亞類中，捆綁（binding）和催化活性（catalytic activity）占比最多，分別占比33.3%和26.8%。在生物過程的23個功能亞類中，代謝過程（metabolic process）和細胞過程（cellular process）占比較多，分別為31.0%和30.4%。總的來看，大量unigenes分布在細胞、細胞組分、細胞過程和代謝過程等亞類中。

圖5 高羊茅unigenes GO注釋分類Fig.5 Gene ontology（GO）classification of unigenes in F.arundinacea

2.2.4高羊茅轉錄組KEGG代謝通路分析 通過對比KEGG數據庫，高羊茅轉錄組測序結果顯示：共有15146條unigenes獲得了注釋，占總數的19.45%，分布在31個代謝通路上（圖6）。按照獲得注釋的unigenes數目進行排序，翻譯（translation，1648 unigenes，10.88%）所獲注釋最多，然后依次是信號轉導（signal transduction，1410 unigenes，9.31%）、碳水化合物代謝（carbohydrate metabolism，1404 unigenes，9.27%）、折疊，分類和降解（folding，sorting and degradation，1350 unigenes，8.91%）、氨基酸代謝（amino acid metabolism，1103 unigenes，7.28%）、運輸和分解代謝（transport and catabolism，917 unigenes，6.05%）和內分泌系統（endocrine system，861unigenes，5.68%）等代謝通路。感覺系統（sensory system，35 unigenes，0.23%）和信號分子與互作（signaling molecules and interaction，12 unigenes，0.08%）所獲注釋較少。獲得注釋的代謝通路幾乎代表了高羊茅發育的所有代謝過程，為進一步開展高羊茅的功能基因組學研究提供了豐富的數據資源。

圖6 高羊茅unigenes KEGG注釋分類Fig.6 KEGG classif ication of unigenes in F.arundinacea

2.3 高羊茅差異表達轉錄本分析

2.3.1高羊茅DEGs的數目分析 基于以上的高羊茅轉錄組測序獲得的從頭拼接組裝序列，我們分別對“K31”和“Regenerate”的分蘗節部位進行RNA-Seq分析，以期揭示調控兩個高羊茅品種生長發育的DEGs。以差異倍數＞1或＜—1且P值≤0.05為篩選標準，共獲得3014條DEGs。其中“Regenerate”相較于“K31”有1573條基因上調，1441條基因下調。

2.3.2差異基因GO分類與功能富集 對兩個高羊茅品種的DEGs進行GO富集分析結果顯示：3014條DEGs的功能注釋被分為生物過程、分子功能和細胞組分三大類。三大功能分類又可進一步劃分為68個功能亞類。其中生物過程包括34個功能亞類；細胞組分包括16個功能亞類；分子功能最少，只有18個功能亞類（圖7僅顯示三大類中富集最顯著的10個亞類）。在生物過程所涵蓋的34個亞類中，以細胞組分生物發生（cellular component biogenesis）、細胞大分子復合物亞基組織（cellular macromolecular complex subunit organization）、DNA包裝（DNA packaging）、染 色 體 組 織（chromosome organization）和 染 色 質 組 裝 或 拆 卸（chromatin assembly or disassembly）等10個亞類最為富集。細胞組分的16個功能亞類中，細胞內非膜細胞器（intracellular non-membranebounded organelle）、蛋白質-DNA復合物（protein-DNA complex）、核小體（nucleosome）、染色體組分（chromosomal part）和染色質（chromatin）等7個亞類最為富 集，其 次是細胞組 分（cell part）、細 胞（cell）和大 分子復合物（macromolecular complex）3個亞類。分子功能的18個功能亞類中，以肽鏈內切酶抑制劑活性（endopeptidase inhibitor activity）、酶抑制劑活性（enzyme inhibitor activity）、天冬酰胺酶活性（asparaginase activity）等10個亞類最為富集（圖7）。綜合分析，高羊茅差異基因的GO功能富集，以生物過程相關基因的差異性最為明顯。

圖7 高羊茅DEGs的GO功能分析Fig.7 Functional categorization of assembled unigenes based on gene ontology（GO）classification in crowns of F.arundinacea

2.4 關鍵通路基因轉錄表達分析

2.4.1轉錄因子 轉錄因子作為能調控基因表達功能的蛋白質分子，在植物的生長發育過程中以及應對環境脅迫時都起到了重要的作用。根據BLAST數據分析，許多蛋白可能是來自不同家族的轉錄因子。在所有DEGs中，發現有42個轉錄因子出現差異表達，這些轉錄因子來自不同的家族，包括MADS、TCP、ARF、AP2/ERF和bHLH等多種類型（表5）。MADS-box家族基因功能主要與細胞分化及其相關發育過程有關。有6個MADS轉錄因子存在差異表達，其中有5個轉錄因子（Unigene003445；Unigene008312；Unigene032888；Unigene067434；Unigene049455）上調，1個轉錄因子（Unigene016122）下調。TCP是植物特有的轉錄因子，研究中發現了5個TCP轉錄因子（Unigene032517；Unigene063740；Unigene010687；Unigene022674；Unigene010684）存在差異表達，且全部下調。3個ARF轉錄因子存在差異表達，其中有2個轉錄因子（Unigene027288；Unigene003404）下調，1個轉錄因子（Unigene003016）上調。AP2/ERF轉錄因子是植物特有的一類轉錄因子家族，1個AP2轉錄因子（Unigene038874）下調，6個ERF轉錄因子存在差異表達，其中5個轉錄因子（Unigene077012；Unigene066363；Unigene000226；Unigene059623；Unigene051413）上調，1個轉錄因子（Unigene016424）下調。植物bHLH轉錄因子參與調解各種信號轉導及合成代謝途徑，共有3個bHLH轉錄因子呈現差異表達，2個轉錄因子（Unigene067511；Unigene067617）上調，1個轉錄因子（Unigene017047）下調。

表5 差異表達轉錄因子Table 5 Differentially expressed transcription factors

2.4.2激素代謝及信號轉導 本研究中，有85條unigenes與脫落酸（abscisic acid，ABA）、生長素、乙烯（ethylene，ETH）、細胞分裂素（cytokinin，CTK）、赤霉素、茉莉酸（jasmonic acid，JA）、油菜素甾醇和水楊酸（salicylic acid，SA）等8種植物激素相關，并且在“Regenerate”和“K31”之間存在差異表達。其中與ABA、IAA和ETH相關的unigenes條數較多。還有8條unigenes參與植物激素互作，IAA-ETH（Unigene003404）、ABA-ETH（Unigene066363、Unigene000226）、ETH-JA（Unigene023000）、ABA-GA（Unigene012266）、ABA-JA（Unigene 046376）、IAA-BR-CTK（Unigene016424）和IAA-ETH-JA（Unigene058662）。在ABA相關的DEGs中有18個基因上調，8個基因下調。下調基因中有3個醛氧化酶基因（Unigene072636、Unigene030806、Unigene 074427）。上調基因中有2個富含半胱氨酸的受體樣蛋白激酶（Unigene060943、Unigene034794）、2個AP2/ERF轉錄因子（Unigene066363、Unigene000226）和1個轉錄抑制因子（Unigene000154）。與IAA相關的DEGs中有10個上調基因，9個下調基因。其中有1個生長素原初響應基因（Aux/IAA）（Unigene046724）和1個核糖體蛋白基因（Unigene046971）下調；上調基因中有1個糖基水解酶基因（Unigene074290）和1個dof蛋白（Unigene016356）。3個生長素應答因子中2個基因下調（Unigene003404、Unigene027288），1個基因上調（Unigene003016）。與ETH相關的DEGs有12個上調基因，9個下調基因。其中有3個脫氫酶基因（Unigene003249、Unigene071207、Unigene071208）和1個ACC氧化酶基因（Unigene064634）上調；2個剪接因子基因（Unigene004991、Unigene 004995）和1個熱激蛋白基因（Unigene037362）下調。與GA相關的6個DEGs全部上調表達，其中有3個擴張素前體（Unigene037534、Unigene031989、Unigene074120）和1個開花促進因子（Unigene026443）。之外還有7個DEGs與JA相關（3個基因上調，4個基因下調）；3個DEGs與CTK相關（2個上調，1個下調）；2個DEGs與SA相關（1個上調，1個下調）；1個DEGs與BR相關（1個下調）（表6）。

表6 植物激素相關差異表達基因Table 6 Plants hormone related DEGs

續表Continued Table

2.5 DEGs的熒光定量PCR驗證

為了驗證通過RPKM值統計比較得到的DEGs數據，挑選5種類型（轉錄因子、激素相關基因、擴張素基因、致病相關蛋白和未知功能基因），共24個DEGs進行實時定量PCR（real-time quantitative PCR，qRT—PCR）驗證。選擇FaACTIN4（ACTIN-related protein 4）作為熒光定量PCR驗證的內參基因。qRT—PCR檢測到的所有基因的表達水平與RNA-Seq數據具有較強的一致性（圖8），表明這些RNA-Seq數據準確有效，可用于高羊茅的基因表達圖譜分析。

圖8 DEGs的熒光定量PCR驗證Fig.8 Validation of RNA-seq by q RT-PCR

3 討論

Illumina Hiseq高通量轉錄組測序分析結果顯示：共獲得77872條unigenes，與NR、Uniprot、COG、GO和KEGG等數據庫進行比對，共60552條unigenes獲得功能注釋。此外從“Regenerate”和“K31”兩個高羊茅品種的分蘗節部位篩選出3014個DEGs，“Regenerate”相較于“K31”，上調基因數為1573個，下調基因數為1441個。這些上調和下調的基因對于揭示高羊茅參與株高和分蘗等形態發育的調控基因及其相關代謝途徑具有重要作用。通過GO富集分析，明確了DEGs的基因功能以及參與調控的信號通路。

本研究中發現42個存在差異表達的轉錄因子，包括MADS、TCP、MYB、WRKY、bHLH等19種類型。其中，MADS家族基因廣泛存在于植物、動物和真菌中，它們在生長發育和信號轉導中起著重要作用。在植物中，它們幾乎參與到植物發育的各個方面，最主要體現在調控開花時間和花器官形成方面，以及根和種子的發育。著名的ABCDE模型解釋了在花的不同部位，不同轉錄因子對花器官形成的作用，其中許多轉錄因子來自MADS家族［23］。此外，前人研究發現OsMADS57可以影響植物分蘗/分枝的形成，它可以調節D14（獨腳金內酯受體），并與OsTB1相互作用，共同調控水稻分蘗［24］。MADS轉錄因子還可以通過調控植物體內激素水平影響植物發育，OsMADS22負調控BR的合成［25］，而BR在植物生長發育過程中發揮重要的作用，在植物的不同器官均含有BR［26］，它可以調控細胞伸長和分裂、光形態建成和維管束發育等過程［27］。BR受體過表達擬南芥出現葉片表皮和葉肉細胞體積增大現象［28］，本研究中“Regenerate”的表皮細胞體積小于“K 31”，且BR相關差異基因（Unigene016424）呈下調模式，這可能與BR對葉片細胞體積的調控相關。關于MADS參與草類植物生長發育的報道相對較少，本研究所鑒定的6個MADS-box unigenes是否影響了兩種高羊茅在表型上的差異還有待進一步研究。

TCP是植物特有的轉錄因子，以玉米的TEOSINTE BRANCHED 1（TB1）、金魚草（Antirrhinum majus）的CYCLOIDEA（CYC）和水稻的PROLIFERATING CELL FACTORS1/2（PCF 1、PCF 2）的首字母命名。TCP在調控植物細胞周期及根、莖、葉、花、種子等器官發育過程中發揮著重要的作用［29］。TCP家族通常在快速生長的組織器官中表達較多，說明這種轉錄因子與植物細胞的分裂分化與組織器官的發育有密切關系。在擬南芥中，AtTCP20能夠與CYCB1的基因元件相結合，進而調控細胞的分裂與生長［30］。AtTCP18是玉米TB1的同源基因，能夠抑制腋芽的發生，AtTCP18缺失突變體呈現大量側枝的表型，并且參與生長素調控側枝發育的信號網絡［9］。AtTCP14和AtTCP15參與細胞分裂素的信號轉導途徑，并且AtTCP14和AtTCP15雙突變體植株出現矮化性狀，這與GA-DELLA的信號轉導途徑相關［31—32］。從我們的RNA-Seq數據可知，TCP轉錄因子的表達量全部下調，結合“Regenerate”相較于“K 31”分蘗更多，植株更矮的表型，TCP轉錄因子的差異表達可能是兩種高羊茅出現表型差異的原因之一。

WRKY是高等植物中第七大轉錄因子家族，普遍參與植物的生長發育以及脅迫響應機制。WRKY轉錄因子家族每一個成員的功能不盡相同，在植物生長發育方面的影響相對復雜，很多研究發現WRKY轉錄因子既可以調控開花時間、促進植株發育也可以誘導植株矮化。在前人的研究中，大豆（Glycine soja）GsWRKY 20能夠控制植物的開花時間，過表達GsWRKY 20植株比野生型開花更早［33］。OsWRKY 53缺失突變體出現葉傾角變小和植株變矮等生長受阻性狀［34］，但過表達OsWRKY 78可促進水稻莖稈伸長［35］。擬南芥中過表達OsWRKY 72植物出現側枝生長增強的現象［36］。在擬南芥中發現AtWRKY 71能夠促進開花并且調節植株分枝發育［37］，過表達AtWRKY 15可促使葉片擴增［38］。在棉花（Gossypium hirsutum）中過表達GhWRKY 15會出現莖稈伸長的性狀［39］。但是，過表達OsWRKY 11水稻出現了植株矮化的表型，這一過程主要受CTK、IAA、GA和BR的調控［40］。過表達OsWRKY 13［41］和OsWRKY 89［42］使水稻生長發育遲緩，植株變矮。在我們的研究中共發現2個差異表達WRKY轉錄因子，均呈現上調的表達模式，結合“K 31”和“Regenerate”的性狀表現，在株高和分蘗上的差異可能與WRKY轉錄因子相關，但是由于WRKY轉錄因子家族成員眾多，且在不同的物種中起到的作用也不盡相同，WRKY在高羊茅生長發育過程起到的作用還需進一步挖掘。

對多種物種的研究表明，芽中ABA的含量與芽活力呈負相關［43］。外源ABA抑制了豌豆（Pisum sativum）［44］、擬南芥［45］和番茄芽的生長［46］，ABA生物合成抑制劑促進了月季（Rosa hybrida）芽的生長［47］。有研究認為，ABA與玉米素核苷（trans-zeatin-riboside，ZR）可協同調控分蘗發育，ZR/ABA低時，玉米分蘗發育受到抑制，ZR/ABA高時，促進玉米分蘗的發育［48］。IAA和CTK共同調控植物分蘗芽的休眠與發育，兩者的調控作用相反。IAA的極性運輸被認為是抑制分枝發育的主要因素，可以抑制芽的起始［49］和伸長［50］。CTK可以解除頂端優勢，誘導花芽分化，促進分蘗芽的萌發。Aux/IAA是生長素原初響應基因，可以對生長素做出快速響應。在水稻中過表達OsIAA 1導致水稻對生長素敏感性減弱，出現植株矮化，側根增多的現象［51］。在擬南芥中過表達楊樹（Populus）Ptr IAA 14.1基因，植株呈現矮化和側枝增多的表型［52］。前人的多項研究表明ABA和IAA相關基因可以參與調控植株高度和分蘗發育，在研究中注釋到的ABA相關DEGs主要參與ABA代謝和信號轉導，IAA相關DEGs不參與代謝，主要參與IAA的信號轉導。因此可以在之后的研究中從激素代謝和信號轉導的角度進一步探討具體是哪些DEGs直接參與影響“K31”和“Regenerate”之間的表型差異。

GA可以促進種子萌發、莖稈伸長、葉片伸展和抑制分蘗芽的發育。前人研究發現，在水稻中過表達OsGA 2ox1［53］、OsGA 2ox5［54］和OsGA 2ox6［55］會使植物體內活性GA含量下降，出現植株變矮、分蘗增多的性狀。外源噴施GA還可以有效控制植物分蘗的發生［56］。本研究共得到6個與GA相關的DEGs，全部呈上調表達模式。其中有3個是擴張素前體，1個是糖基水解酶。擴張素蛋白能夠改變細胞壁結構，使細胞壁松弛，促進細胞伸展與分裂；糖基水解酶在植物中主要起到調控生長發育和響應環境脅迫的功能［57］。這些基因間的差異表達可能是“K 31”和“Regenerate”表型差異的部分原因。GA的信號轉導通路與GRAS家族的轉錄因子密切相關，GRAS家族轉錄因子參與了植物側枝形成、光信號轉導、莖稈發育和赤霉素信號轉導等多種生理生化過程。DELLA蛋白是GRAS蛋白家族中的亞家族，參與GA與其他多種激素的信號轉導，在調控植物生長發育以及應對環境脅迫中具有重要意義。DELLA蛋白可以負向調控GA的信號轉導抑制植物的生長，而GA通過降解DELLA蛋白起到促進植物生長的作用［16］。本研究注釋到1個GRAS家族轉錄因子，在水稻的研究中發現OsCIGR（chitin-inducible gibberellin-responsive）在高植株中的表達量高于矮植株［58］。該轉錄因子在“Regenerate”較“K31”的差異分析中呈下調表達模式，這與“Regenerate”和“K31”的株高性狀一致。所以可以推測CIGR轉錄因子可能是引起“Regenerate”和“K31”株高產生差異的原因之一。

4 結論

本研究通過Illumina Hiseq高通量測序技術構建了不同高羊茅品種的轉錄組文庫，初步獲得了“K 31”和“Regenerate”的差異基因，挖掘到MADS、TCP和WRKY等差異表達轉錄因子以及眾多與ABA、IAA和GA等激素代謝和信號轉導相關的DEGs。為探究高羊茅生長發育的分子機制提供了基礎數據和理論支持，后續可對還未明確功能和轉導通路的候選基因進行進一步的功能驗證。