不同刈割強(qiáng)度下羊草轉(zhuǎn)錄組研究

趙勁博,侯向陽(yáng)*,武自念,任衛(wèi)波,胡寧寧,郭豐輝,馬文靜(1.中國(guó)農(nóng)業(yè)科學(xué)院草原研究所,內(nèi)蒙古 呼和浩特 010010;2. 農(nóng)業(yè)部草地生態(tài)與修復(fù)治理重點(diǎn)實(shí)驗(yàn)室，內(nèi)蒙古 呼和浩特 010010)

羊草(Leymuschinensis)屬于禾本科賴草屬多年生根莖型草本植物，是歐亞大陸溫帶草原東部的重要建群種之一[1]，具有抗旱、抗寒，耐鹽堿、耐貧瘠等特點(diǎn)，在維持中國(guó)北方草原生態(tài)系統(tǒng)穩(wěn)定、生物多樣性、人工草地改良等方面有著重要意義。羊草持綠期長(zhǎng)、葉量多，適口性好，營(yíng)養(yǎng)價(jià)值高，其粗蛋白含量明顯高于其他禾本科牧草[2]，為多種牲畜喜食的高質(zhì)量牧草。此外，羊草草地還是很好的割草地，具有很強(qiáng)的再生性，可為牲畜越冬提供優(yōu)良飼草。

刈割是草地利用的重要方式之一，也是人為影響草地生態(tài)環(huán)境的重要因素。刈割減少了植物光合作用的面積，降低了碳和營(yíng)養(yǎng)物質(zhì)的積累，直接影響草原生產(chǎn)力[3]。刈割制度對(duì)飼草的再生生理和形態(tài)都會(huì)產(chǎn)生強(qiáng)烈的影響[4]，而不同刈割強(qiáng)度會(huì)對(duì)刈割效果和牧草再生狀態(tài)產(chǎn)生不一樣的作用，刈割留茬高度與牧草的干草產(chǎn)量與損耗率有著直接的關(guān)系[5]。有研究表明，隨著刈割強(qiáng)度的增大，羊草生物量、高度和密度呈下降趨勢(shì)[6]。比起一些較為矮小的草原植物，羊草對(duì)刈割具有更強(qiáng)的敏感性，程度較大的機(jī)械損傷對(duì)其生長(zhǎng)可造成更嚴(yán)重的影響[7]。刈割高度對(duì)羊草群落的影響在生理生態(tài)上研究較為廣泛，但在個(gè)體的分子水平上還涉及甚少。

自從人們發(fā)現(xiàn)RNA在基因組和蛋白組研究中的關(guān)鍵作用，轉(zhuǎn)錄本識(shí)別和基因表達(dá)量計(jì)算便逐漸成為了分子生物學(xué)中重要研究方法[8]。隨著高通量測(cè)序技術(shù)的不斷完善，轉(zhuǎn)錄組數(shù)據(jù)挖掘與分析方法日趨成熟，一些沒有基因組數(shù)據(jù)支持而又有著很高經(jīng)濟(jì)與生態(tài)價(jià)值的重要物種可以在測(cè)序成本和時(shí)間相對(duì)較少的條件下得到大批差異基因的表達(dá)情況，極大地幫助了相關(guān)物種生物信息數(shù)據(jù)庫(kù)的完善，為基因克隆驗(yàn)證與表型研究提供了充足原始資料。在羊草低溫、干旱、鹽堿、刈割、放牧等脅迫作用下的轉(zhuǎn)錄組研究中，一些抗性基因與關(guān)鍵通路的挖掘與分析，為改良草地培育優(yōu)質(zhì)牧草奠定了基礎(chǔ)[9]。羊草刈割與創(chuàng)傷的相關(guān)基因表達(dá)譜分析結(jié)果表明，刈割激活了茉莉酸生物合成途徑基因的表達(dá)，誘導(dǎo)創(chuàng)傷信號(hào)，但核糖體蛋白基因、細(xì)胞分裂和膨脹相關(guān)基因及木質(zhì)素生物合成基因表達(dá)卻受到了抑制，影響羊草早期的恢復(fù)生長(zhǎng)，而創(chuàng)傷處理下羊草脯氨酸合成的基因上調(diào)，引起脯氨酸的積累[10]。在對(duì)羊草刈割傷口對(duì)牛血清蛋白(BSA)響應(yīng)的轉(zhuǎn)錄組研究中，有關(guān)細(xì)胞氧化、細(xì)胞程序性死亡、葉片卷曲動(dòng)力和刈割創(chuàng)傷修復(fù)的基因有特異表達(dá)[11]。不同刈割處理對(duì)羊草的生長(zhǎng)和利用產(chǎn)生重要影響，不同刈割強(qiáng)度轉(zhuǎn)錄組研究對(duì)于了解羊草響應(yīng)刈割的關(guān)鍵分子機(jī)制及關(guān)鍵基因挖掘與應(yīng)用有指導(dǎo)性意義。

本研究通過(guò)羊草在留茬高度分別為0、2、4、8、12 cm這5個(gè)不同強(qiáng)度刈割后再生葉片的轉(zhuǎn)錄組測(cè)序，RNA-seq數(shù)據(jù)組裝，注釋，以及對(duì)差異基因查找、篩選和功能與代謝途徑富集分析，挖掘刈割相關(guān)的基因，并且根據(jù)基因表達(dá)量與刈割強(qiáng)度的相關(guān)性，提取出表達(dá)量隨刈割強(qiáng)度的增強(qiáng)而持續(xù)升高或降低的基因。通過(guò)分析與刈割有關(guān)的關(guān)鍵基因，從分子層面探索了羊草響應(yīng)不同程度機(jī)械損傷的機(jī)理，豐富了植物抗逆基因的數(shù)據(jù)庫(kù)，為刈割對(duì)羊草在生理生化方面的影響提供研究方向，而且為培養(yǎng)牧草及其他經(jīng)濟(jì)作物的耐刈優(yōu)良品種奠定理論基礎(chǔ)，從而為改善生態(tài)環(huán)境、發(fā)展人工草地、提高農(nóng)牧業(yè)生產(chǎn)力提供幫助。

1 材料與方法

1.1 試驗(yàn)材料

試驗(yàn)材料于2014年5月采集于中國(guó)農(nóng)業(yè)科學(xué)院草原研究所錫林浩特放牧試驗(yàn)站，(44°15′ N，116°42′ E，海拔1100 m)，小心將單株羊草地上部分連同地下根部挖出，移至培養(yǎng)室25 ℃，光照16 h，黑暗8 h條件下土培(營(yíng)養(yǎng)土∶蛭石=3∶1)擴(kuò)繁，挑選生長(zhǎng)大小一致、狀況良好的分蘗株移入多個(gè)盆中，每盆6～10株，培養(yǎng)備用。

1.2 試驗(yàn)設(shè)計(jì)

2016年1月開始，對(duì)分蘗株沿土壤表面刈割，進(jìn)行室內(nèi)土培并控制培養(yǎng)條件，保證各單株處于同一發(fā)育時(shí)期，同一生長(zhǎng)環(huán)境。當(dāng)植株生長(zhǎng)14 d后(分蘗期)，挑選生長(zhǎng)狀況一致的30盆植株使用直尺與剪刀進(jìn)行梯度刈割處理并取樣：刈割植株地上部分，分別留茬0，2，4，8，12 cm，依次記為A，B，C，D，E五組，以及不刈割對(duì)照組O，刈割3 d后取刈割組地上再生葉片以及對(duì)照組的全部葉片，每組6～10株植株(單盆)混合取樣保證測(cè)序上樣量，取樣后立即液氮速凍，-80 ℃保存?zhèn)溆谩?/p>

1.3 RNA提取及轉(zhuǎn)錄組測(cè)序

采用TRIzol法提取羊草葉片總RNA，并用瓊脂糖凝膠電泳、Nanodrop 2000分析RNA的純度和完整性，樣品檢測(cè)合格后，送至北京諾禾致源生物信息科技有限公司采用Illumina HiSeq-PE150進(jìn)行測(cè)序。

1.4 數(shù)據(jù)組裝與蛋白編碼框預(yù)測(cè)

獲得原始數(shù)據(jù)后，去除帶接頭(adapter)的reads；去除N(N表示無(wú)法確定堿基信息)的比例大于0.1%的reads；去除低質(zhì)量reads(質(zhì)量值Qphred≤20的堿基數(shù)占整個(gè)reads的50%以上的reads)。將過(guò)濾后的clean data使用軟件Trinity[12]進(jìn)行無(wú)參組裝，設(shè)置計(jì)算最小kmer的參數(shù)(--min_kmer_cov)為2，最短contig長(zhǎng)度(--min_contig_length)為200 bp，組裝得到轉(zhuǎn)錄本，挑出最長(zhǎng)轉(zhuǎn)錄本作為Unigene，隨后使用軟件TransDecoder[13]預(yù)測(cè)組裝結(jié)果的蛋白編碼框并將得到的蛋白序列用于注釋。

1.5 組裝結(jié)果注釋與功能通路富集

通過(guò)使用比對(duì)軟件BLAST[14]將得到的羊草轉(zhuǎn)錄本核酸序列和蛋白序列與蛋白核酸數(shù)據(jù)庫(kù)進(jìn)行比對(duì)(E≤10-5)，每條序列挑取最好的比對(duì)結(jié)果使用Trinotate軟件合并，得到了Unigene的注釋結(jié)果。用到的數(shù)據(jù)庫(kù)有：Nr(non-redundant protein database，非冗余蛋白數(shù)據(jù)庫(kù))、SwissProt(SwissProt protein database，蛋白質(zhì)序列數(shù)據(jù)庫(kù))、Pfam[15](protein families database，蛋白質(zhì)家族域數(shù)據(jù)庫(kù))、KEGG[16](Kyoto Encyclopedia of Genes and Genomes，東京基因與基因組百科全書)、Uniref90[17](全球蛋白資源數(shù)據(jù)庫(kù)參考資料庫(kù))、KOG(Clusters of Orthologous Groups from 7 eukaryotic complete genomes，7個(gè)真核生物蛋白質(zhì)直系同源數(shù)據(jù)庫(kù))。通過(guò)在線軟件WEGO[18]將得到的GO進(jìn)行統(tǒng)計(jì)分類，使用GOseq[19]以所有Unigene的GO注釋結(jié)果為背景對(duì)上下調(diào)差異基因做GO功能富集分析。利用本地軟件KOBAS 2.0[20]以所有Unigene的KEGG注釋結(jié)果為背景對(duì)所有顯著差異基因的KEGG代謝通路做富集分析。

1.6 差異基因與梯度基因的篩選

首先將測(cè)序得到的各處理clean reads利用比對(duì)軟件Bowtie[21]分別回貼到所有Unigene上，并通過(guò)RSEM[22]計(jì)算每條Unigene在測(cè)序reads上的表達(dá)量。不同處理兩兩之間的顯著性差異基因利用DEGseq[23]軟件包根據(jù)MA-plot方法與隨機(jī)抽樣模型MARS(Random Sampling model)[23]篩挑得到，設(shè)定的差異表達(dá)倍數(shù)的對(duì)數(shù)(log2fc)大于1，P值小于10-3。

梯度基因篩挑方法：先提取所有處理兩兩之間的顯著差異基因，再將其標(biāo)準(zhǔn)化后的表達(dá)量矩陣輸入到用于研究處理梯度樣本表達(dá)量變化的聚類軟件stem[24]中，設(shè)定最大聚類數(shù)(maximum number of model profiles)為50，不同梯度單元之間最大變化(maximum number of units a model profile)為2，從而得到與刈割高度有關(guān)的梯度差異基因。

2 結(jié)果與分析

2.1 羊草轉(zhuǎn)錄組數(shù)據(jù)組裝

Illumina高通量測(cè)序平臺(tái)(HiSeq-PE150)測(cè)序共得到139767803條原始reads，總長(zhǎng)度為41.94 Gb，去除接頭和低質(zhì)量片段后reads數(shù)為137626643條，大小為41.30 Gb，原始測(cè)序數(shù)據(jù)已上傳至NCBI(SRX2792325)。使用Trinity對(duì)所有6個(gè)樣品的轉(zhuǎn)錄組clean reads進(jìn)行無(wú)參組裝，共得到轉(zhuǎn)錄本270207條，總長(zhǎng)度為191632770 bp(191.6 Mb)，平均長(zhǎng)度為709.21 bp，N50為1075 bp(1 kb)。提取最長(zhǎng)轉(zhuǎn)錄本得到Unigene共180770條，總長(zhǎng)度為106133700 bp(106.1 Mb)，平均長(zhǎng)度為587.12 bp，N50為806 bp。轉(zhuǎn)錄本與Unigene長(zhǎng)度分布如圖1所示。

圖1 轉(zhuǎn)錄本與Unigene長(zhǎng)度分布Fig.1 Length distribution of transcripts and Unigenes橫坐標(biāo)為轉(zhuǎn)錄本或Unigene的長(zhǎng)度，縱坐標(biāo)為其數(shù)目。 Abscissa for the length of transcripts or Unigenes, ordinate for their number.

2.2 Unigene注釋結(jié)果

通過(guò)對(duì)Unigene的蛋白預(yù)測(cè)，得到具有蛋白編碼框的Unigene共54555條，然后使用Blastx與Blastp將核酸及蛋白序列與數(shù)據(jù)庫(kù)進(jìn)行比對(duì)，總共注釋到Unigene蛋白48097條，其中注釋到SwissPort的有32913條(60.33%)，Nr有46765條(85.72%)，KOG有38470條(70.52%)，Pfam有27219條(49.89%)，GO有28409條(52.07%)，KEGG有46461條(85.16%)，未注釋出的蛋白有6458條(11.84%)。KOG分類結(jié)果主要與信號(hào)轉(zhuǎn)導(dǎo)機(jī)制(signal transduction mechanisms)、翻譯后修飾和蛋白反轉(zhuǎn)、分子伴侶(posttranslational modification, protein turnover, chaperones)等相關(guān)(圖2)。有關(guān)細(xì)胞、細(xì)胞組分、結(jié)合、催化、細(xì)胞過(guò)程、代謝過(guò)程等結(jié)果在GO分類中較為顯著(圖3)。

圖2 KOG功能分類Fig.2 KOG function classification橫坐標(biāo)為KOG種類，縱坐標(biāo)為基因數(shù)目所占百分比。 Abscissa for KOG classification, ordinate for percent of genes number.

圖3 GO分類Fig.3 GO classification橫坐標(biāo)為GO種類，縱坐標(biāo)為基因數(shù)目所占百分比。 Abscissa for GO classification, ordinate for percent of genes number.

2.3 差異基因數(shù)目分析

根據(jù)羊草刈割不同高度的表達(dá)量矩陣，找到每個(gè)刈割高度處理的實(shí)驗(yàn)組(A、B、C、D、E)與不刈割對(duì)照組(O)表達(dá)量共有顯著差異的基因24829個(gè)，其中A(留茬0 cm)有11618個(gè)，上調(diào)基因有7867個(gè)，下調(diào)基因有3751個(gè)；B(留茬2 cm)有9479個(gè)，上調(diào)基因6190個(gè)，下調(diào)基因3289個(gè)；C(留茬4 cm)有9159個(gè)，上調(diào)基因6031個(gè)，下調(diào)基因3128個(gè)；D(留茬8 cm)有13375個(gè)，上調(diào)基因10415個(gè)，下調(diào)基因2960個(gè)；E(留茬12 cm)有7143個(gè)，上調(diào)基因2777個(gè)，下調(diào)基因4366個(gè)。5個(gè)實(shí)驗(yàn)組與對(duì)照組共有的顯著差異基因2579個(gè)，其中上調(diào)表達(dá)994個(gè)，下調(diào)表達(dá)1585個(gè)，這些基因極有可能與羊草響應(yīng)刈割有關(guān)。

通過(guò)使用DEGseq對(duì)不同留茬高度處理樣本進(jìn)行兩兩比較分析，共找到30231個(gè)顯著差異表達(dá)基因。然后利用stem軟件對(duì)這些基因進(jìn)行聚類，得到隨著刈割強(qiáng)度增強(qiáng)(留茬高度降低)，表達(dá)量隨之規(guī)律變化的差異基因。其中，表達(dá)量呈上升趨勢(shì)模式的差異基因有3499個(gè)(P=0)(圖4中39號(hào)，圖5左)；表達(dá)量呈下降趨勢(shì)模式的差異基因有1245個(gè)(P=0)(圖4中8號(hào)，圖5右)；表達(dá)量呈先上升后下降趨勢(shì)模式的差異基因有338個(gè)(P=5×10-8)(圖4中48號(hào))；表達(dá)量呈先下降后上升的模式結(jié)果不顯著(圖4中9號(hào))。

圖4 所有差異基因stem聚類Fig.4 Stem clustering of all DEGs 曲線代表差異基因隨脅迫程度增加表達(dá)量變化趨勢(shì)，彩色框與無(wú)色框分布代表顯著性和非顯著性聚類，左上角與左下角數(shù)字分別為趨勢(shì)類型編號(hào)和聚類顯著性。 Curve for DEGs expression changed trend with the increase of stress degree. Color and colorless box represent significant and nonsignificant clustering respectively. The number of top left and lower left corner represent trend type and clustering significant respectively.

圖5 隨刈割強(qiáng)度增強(qiáng)表達(dá)量呈上升與下降趨勢(shì)差異基因熱圖Fig.5 Heatmap of DEGs having up and down trend expression with the increase of mowing intensity左圖呈上升趨勢(shì)，右圖呈下降趨勢(shì)，表達(dá)量均經(jīng)過(guò)z-score標(biāo)準(zhǔn)化處理。 Left represents up trend expression, right represents down trend expression, expression are normalized by z-score.

2.4 羊草刈割差異基因功能與代謝通路富集分析

將所有GO注釋結(jié)果作為參考背景，將刈割實(shí)驗(yàn)組與不刈割對(duì)照組上調(diào)和下調(diào)顯著性差異基因的GO分別進(jìn)行富集處理，設(shè)定FDR小于0.05，得到刈割上調(diào)基因GO富集結(jié)果19條(表1)。在細(xì)胞組分方面，上調(diào)基因與胞外區(qū)、質(zhì)外體、葉綠體有關(guān)；在生物途徑上與碳水化合物代謝、損傷應(yīng)答、過(guò)氧化氫分解有關(guān)；在分子功能上與水解酶活性、脂肪酸結(jié)合、蛋白二硫氧化還原酶活性有關(guān)。刈割下調(diào)基因GO富集結(jié)果24條(表2)，在細(xì)胞組分上主要與膜有關(guān)；在生物途徑上與跨膜運(yùn)輸、蛋白磷酸化、水楊酸反應(yīng)有關(guān)；在分子功能上與蛋白激酶活性、多糖結(jié)合、跨膜轉(zhuǎn)運(yùn)體活性有關(guān)。

表1 刈割上調(diào)差異基因GO富集結(jié)果(FDR<0.05)Table 1 GO enrichment of mowing up-regulated DEGs (FDR<0.05)

FDR:錯(cuò)誤發(fā)現(xiàn)率 False discovery rate.

利用KOBAS對(duì)差異基因KEGG代謝通路進(jìn)行富集，得到BH法(Benjamini and Hochberg)校正P-value小于0.05的4條顯著代謝通路(表3)，分別為淀粉和蔗糖代謝、半乳糖代謝、苯丙烷生物合成、植物激素信號(hào)轉(zhuǎn)導(dǎo)。其中，在植物激素信號(hào)轉(zhuǎn)導(dǎo)方面，主要涉及脫落酸代謝中的脫落酸受體PYR/PYL家族(abscisic acid receptor PYR/PYL family)、蛋白磷酸酶2C(protein phosphatase 2C)；茉莉酸代謝中茉莉酮酸酯ZIM結(jié)構(gòu)域蛋白(jasmonate ZIM domain-containing protein)、茉莉素信號(hào)受體蛋白COI1(coronatine-insensitive protein 1)；水楊酸代謝中病程相關(guān)蛋白1(pathogenesis-related protein 1)等。

此外，對(duì)刈割強(qiáng)度相關(guān)的梯度基因進(jìn)行GO與KEGG富集，挑取GO富集結(jié)果(FDR<0.05)以及代謝通路富集結(jié)果(校正P-value小于0.05)，發(fā)現(xiàn)表達(dá)量隨刈割強(qiáng)度持續(xù)上升的基因主要與碳水化合物代謝、氨酰-tRNA合成、卟啉與葉綠素代謝、過(guò)氧化物酶體等有關(guān)(表4)；而表達(dá)量隨刈割強(qiáng)度持續(xù)下降的基因與蛋白激酶活性、葉綠體類囊體膜、光合作用、半胱氨酸與蛋氨酸代謝、硫代謝等有關(guān)(表5)；表達(dá)量隨刈割強(qiáng)度先上升后下降的基因與光合作用的天線蛋白(antenna proteins)、光捕獲、光系統(tǒng)、膜組成、葉綠素等有關(guān)。

表2 刈割下調(diào)差異基因GO富集結(jié)果(FDR<0.05)Table 2 GO enrichment of mowing down-regulated DEGs (FDR<0.05)

表3 刈割差異基因KEGG富集(校正P-value<0.05)Table 3 KEGG enrichment of mowing DEGs (Corrected P-value<0.05)

表4 表達(dá)量隨刈割強(qiáng)度增強(qiáng)呈上升趨勢(shì)差異基因KEGG富集(校正P-value<0.05)Table 4 KEGG enrichment of DEGs having up trend expression with the increase of mowing intensity (Corrected P-value<0.05)

表5 表達(dá)量隨刈割強(qiáng)度增強(qiáng)呈下降趨勢(shì)差異基因KEGG富集(校正P-value<0.05)Table 5 KEGG enrichment of DEGs having down trend expression with the increase of mowing intensity (Corrected P-value<0.05)

圖6 刈割差異基因脫落酸與茉莉酸信號(hào)轉(zhuǎn)導(dǎo)代謝通路富集結(jié)果Fig.6 ABA and JA signal transduction metabolic pathway enrichment of mowing DEGs紅色代表上調(diào)，綠色代表下調(diào)。Red represents up-regulated, green represents down-regulated.

3 討論

隨著高通量測(cè)序技術(shù)的成熟與完善，越來(lái)越多無(wú)參考基因組的物種實(shí)現(xiàn)了轉(zhuǎn)錄組測(cè)序，并成功完成了在各種試驗(yàn)處理下物種轉(zhuǎn)錄表達(dá)水平變化的研究。羊草作為基因組數(shù)據(jù)龐大的異源四倍體草本植物，并且是一種具有抗寒抗旱耐鹽耐刈等多種抗逆性的優(yōu)良牧草，大量挖掘其關(guān)鍵功能基因是一項(xiàng)重要而艱巨的任務(wù)。目前與羊草抗逆相關(guān)的轉(zhuǎn)錄組研究已經(jīng)取得了很多的突破性進(jìn)展[10,25-26]。本研究通過(guò)不同刈割強(qiáng)度下羊草RNA轉(zhuǎn)錄表達(dá)情況，找到了與刈割相關(guān)的基因，并對(duì)其中部分基因的功能與代謝進(jìn)行了進(jìn)一步的分析研究。

刈割對(duì)植物在表型上造成可見的機(jī)械損傷，在其內(nèi)部的分子水平上也會(huì)有不同程度的影響。在面對(duì)環(huán)境壓力或生物脅迫時(shí)，植物會(huì)通過(guò)產(chǎn)生一些信號(hào)分子，在細(xì)胞內(nèi)部發(fā)生關(guān)鍵的轉(zhuǎn)錄變化進(jìn)而改變其生理狀態(tài)[27]，以應(yīng)對(duì)外來(lái)的刺激與干擾。植物激素如生長(zhǎng)素、脫落酸、茉莉酸、水楊酸等在植物抵御生物或非生物脅迫中起著重要作用。參與茉莉素形成的各種復(fù)合物廣泛分布于植物中，影響著諸如花粉、果實(shí)形成，根系生長(zhǎng)，莖葉卷曲等生長(zhǎng)發(fā)育過(guò)程[28]，茉莉素介導(dǎo)的抗性信號(hào)途徑參與了有關(guān)機(jī)械損傷、病原菌侵染、昆蟲噬咬等應(yīng)對(duì)脅迫的抗逆反應(yīng)[29]。COI1是一種茉莉素信號(hào)受體蛋白[30]，可與其他多種蛋白組合形成SCFCOI1泛素連接酶復(fù)合體[31]，從而促進(jìn)茉莉素信號(hào)響應(yīng)；JAZ蛋白家族是茉莉素信號(hào)途徑中的一類抑制蛋白,是SCFCOI1泛素連接酶復(fù)合體的一類直接底物[32]。而本研究結(jié)果顯示(圖6)刈割處理導(dǎo)致COI1基因表達(dá)量下調(diào)，編碼JAZ蛋白家族基因表達(dá)量上調(diào)，預(yù)測(cè)刈割對(duì)處理3 d后的羊草莖葉的茉莉素信號(hào)傳遞可能有抑制作用，從而對(duì)植物生長(zhǎng)發(fā)育造成影響。Chen等[10]在羊草轉(zhuǎn)錄組試驗(yàn)的結(jié)果顯示刈割處理后24 h內(nèi)茉莉酸合成途徑相關(guān)基因表達(dá)被激活，由此推測(cè)茉莉酸調(diào)控羊草在一定時(shí)間刈割脅迫下的抗逆反應(yīng)可能主要依靠促進(jìn)茉莉酸合成來(lái)實(shí)現(xiàn)，而茉莉素信號(hào)轉(zhuǎn)導(dǎo)途徑卻不一定會(huì)被激活甚至是被抑制。脫落酸在植物種子萌發(fā)、芽的休眠、氣孔開閉等生理過(guò)程中起著調(diào)控的作用，也介導(dǎo)了許多植物的環(huán)境逆境反應(yīng)[33]。PYR/PYL/RCAR是ABA的一類受體蛋白[34-35]，在脫落酸信號(hào)轉(zhuǎn)導(dǎo)途徑中起到正向調(diào)節(jié)的作用[36]，其與脫落酸結(jié)合后，會(huì)使有負(fù)調(diào)控作用的PP2C蛋白活性被抑制[37]，進(jìn)而促進(jìn)脫落酸信號(hào)轉(zhuǎn)導(dǎo)。經(jīng)過(guò)對(duì)數(shù)據(jù)的分析顯示(圖6)，與不刈割相比，刈割處理使編碼羊草PYR/PYL/RCAR蛋白家族的基因表達(dá)量呈現(xiàn)上調(diào)趨勢(shì)，而編碼PP2C蛋白的基因表達(dá)量下調(diào)，說(shuō)明刈割可能會(huì)導(dǎo)致脫落酸更易與受體結(jié)合，進(jìn)而促使脫落酸反應(yīng)。從刈割上調(diào)差異基因的GO富集結(jié)果中得到了一些與損傷響應(yīng)(GO:0009611)有關(guān)的基因，這些基因大多與擬南芥(Arabidopsisthaliana)應(yīng)對(duì)損傷的基因有相似性，如參與乙烯生物合成的ACS6(At4g11280)基因和與衰老有關(guān)的SEN1(At4g35770)基因，前者在擬南芥機(jī)械損傷試驗(yàn)[38]下鑒定為非依賴于COI基因上調(diào)表達(dá)，后者為依賴COI基因而抑制表達(dá)，這與本試驗(yàn)結(jié)果COI1基因表達(dá)量下調(diào)相符，說(shuō)明這兩個(gè)基因與COI1基因的表達(dá)關(guān)系在羊草與擬南芥中可能是相同的。羊草刈割處理下的上調(diào)基因還比對(duì)到了擬南芥的FAR1(At5g22500)基因、CYP94C1(At2g27690)基因、TIFY11B(At1g72450)基因、WRKY40(At1g80840)基因、BSMT1(At3g11480)基因，這些基因在損傷處理下的擬南芥中也為上調(diào)表達(dá)[39-43]。

由于刈割強(qiáng)度的差別，可能導(dǎo)致植物體內(nèi)相關(guān)基因的轉(zhuǎn)錄表達(dá)出現(xiàn)差異。隨著刈割強(qiáng)度的增強(qiáng)，部分基因表達(dá)量呈現(xiàn)上升趨勢(shì)，如超氧化物歧化酶(SOD)、過(guò)氧化氫酶(CAT)、過(guò)氧化物酶(POD)等與抗氧化酶活性有關(guān)的基因，這說(shuō)明羊草在遭受刈割這種逆境脅迫后，會(huì)通過(guò)加強(qiáng)自身表達(dá)合成過(guò)氧化物酶的過(guò)程來(lái)清除有害物質(zhì)，且刈割強(qiáng)度越高這種抗性越強(qiáng)，類似的結(jié)果在不同損傷高度處理冷蒿(Artemisiafrigida)的研究[38]中也有體現(xiàn)，該研究表示隨著機(jī)械損傷強(qiáng)度增加，冷蒿葉片的SOD、CAT和POD活性也隨之增加，可見植物會(huì)通過(guò)加強(qiáng)自身的抗氧化能力來(lái)適應(yīng)更高強(qiáng)度的損傷脅迫。而部分基因則隨著刈割強(qiáng)度增強(qiáng)表達(dá)量逐漸下降，如與光合作用有關(guān)的基因，這可能是由于刈割導(dǎo)致了羊草葉片減少，且留茬高度越低，其葉面積就越小，進(jìn)而造成光合作用逐漸下降，相關(guān)基因的表達(dá)也受到了進(jìn)一步限制，如與葉綠體類囊體膜、光系統(tǒng)Ⅱ復(fù)合體、光合電子傳遞鏈等結(jié)構(gòu)與功能有關(guān)的基因表達(dá)量呈下降趨勢(shì)。

本研究針對(duì)羊草刈割強(qiáng)度設(shè)置了不同留茬高度作為處理實(shí)驗(yàn)組，通過(guò)轉(zhuǎn)錄組測(cè)序得到大量的轉(zhuǎn)錄本數(shù)據(jù)。在常規(guī)分析的基礎(chǔ)下，重點(diǎn)對(duì)注釋到其他物種已知序列的關(guān)鍵基因進(jìn)行了研究，而對(duì)于一些功能不明確或沒有注釋結(jié)果但表達(dá)量有規(guī)律性的基因，還有待在結(jié)構(gòu)與功能等方面進(jìn)行進(jìn)一步的分析探索，從而使數(shù)據(jù)挖掘與利用更深入全面，以助于探求羊草在機(jī)械損傷逆境下最本質(zhì)的機(jī)理。

4 結(jié)論

通過(guò)轉(zhuǎn)錄組測(cè)序與生物信息學(xué)分析，獲得總長(zhǎng)度達(dá)192.6 Mb，N50為1075 bp的轉(zhuǎn)錄本270207條，經(jīng)過(guò)蛋白編碼框預(yù)測(cè)與序列注釋，共得到有注釋結(jié)果的Unigene 48097條。得到刈割實(shí)驗(yàn)組與不刈割對(duì)照組上調(diào)基因994個(gè)，主要富集到了碳水化合物代謝過(guò)程、損傷應(yīng)答、過(guò)氧化氫分解過(guò)程等功能；下調(diào)基因1585個(gè)，主要富集到了跨膜運(yùn)輸、蛋白磷酸化、水楊酸反應(yīng)等功能。隨刈割強(qiáng)度增強(qiáng)，表達(dá)量呈上升趨勢(shì)的基因有3499個(gè)，富集分析結(jié)果顯示刈割強(qiáng)度可能與氨酰-tRNA生物合成、脂肪酸生物合成、卟啉與葉綠素代謝、氨基糖和核苷酸糖代謝、過(guò)氧物酶體有關(guān)；隨刈割強(qiáng)度增強(qiáng)，表達(dá)量呈下降趨勢(shì)的基因有1245個(gè)，主要與光合作用、半胱氨酸和蛋氨酸代謝、硫代謝有關(guān)。

References：

[1] Li Y H. The divergence and convergence ofAneurolepidiumchinesesteppe andStipgrandissteppe under the grazing influence in Xilin river valley, Inner Mongolia. Chinese Journal of Plant Ecology, 1988, 12(3): 189-196.

李永宏. 內(nèi)蒙古錫林河流域羊草草原和大針茅草原在放牧影響下的分異和趨同. 植物生態(tài)學(xué)報(bào), 1988, 12(3): 189-196.

[2] Zhao H. The Study of Nutrition Value of TheLeymuschinensisfor Dairy Cattle and Assiociative Effects of Feeds. Daqing: Heilongjiang Bayi Agricultural University, 2008.

趙鶴. 羊草對(duì)奶牛營(yíng)養(yǎng)價(jià)值及其日糧組合效應(yīng)的研究. 大慶： 黑龍江八一農(nóng)墾大學(xué), 2008.

[3] Ferraro D O, Oesterheld M. Effect of defoliation on grass growth. A quantitative review. Oikos, 2002, 98(1): 125-133.

[4] Davidson J L, Milthorpe F L. Leaf growth in dactylis glomerate following defoliation. Annals of Botany, 1966, 30(2): 173-184.

[5] Wang Z F, Wang D J, Yu H Z,etal. Effects of cutting time and stubble height on hay yield and quality ofLeymuschinensismeadow. Pratacultural Science, 2016, (2): 276-282.

王志鋒, 王多伽, 于洪柱, 等. 刈割時(shí)間與留茬高度對(duì)羊草草甸草產(chǎn)量和品質(zhì)的影響. 草業(yè)科學(xué), 2016, (2): 276-282.

[6] Zhong Y K, Bao Q H. The effects of different mowing intensity on natural grassland. Grassland of China, 1999, (5): 15-18.

仲延凱, 包青海. 不同刈割強(qiáng)度對(duì)天然割草地的影響. 中國(guó)草地, 1999, (5): 15-18.

[7] Han L, Guo Y J, Han J G,etal. A study on the diversity and aboveground biomass in aLeymuschinensismeadow steppe community under different cutting intensities. Acta Prataculturae Sinica, 2010, 19(3): 70-75.

韓龍, 郭彥軍, 韓建國(guó), 等. 不同刈割強(qiáng)度下羊草草甸草原生物量與植物群落多樣性研究. 草業(yè)學(xué)報(bào), 2010, 19(3): 70-75.

[8] Conesa A, Madrigal P, Tarazona S,etal. A survey of best practices for RNA-seq data analysis. Genome Biology, 2016, 17(1): 181.

[9] Hou X Y. Advances and prospects of grassland plant basic biology research. China Basic Science, 2016, 18(2): 67-76.

侯向陽(yáng). 草原植物基礎(chǔ)生物學(xué)研究進(jìn)展與展望. 中國(guó)基礎(chǔ)科學(xué), 2016, 18(2): 67-76.

[10] Chen S, Cai Y, Zhang L,etal. Transcriptome analysis reveals common and distinct mechanisms for sheepgrass (Leymuschinensis) responses to defoliation compared to mechanical wounding. Plos One, 2014, 9(2): e89495.

[11] Huang X, Peng X, Zhang L,etal. Bovine serum albumin in saliva mediates grazing response inLeymuschinensisrevealed by RNA sequencing. BMC Genomics, 2014, 15(1): 1126.

[12] Grabherr M G, Haas B J, Yassour M,etal. Full-length transcriptome assembly from RNA-Seq data without a reference genome. Nature Biotechnology, 2011, 29(7): 644-652.

[13] Haas B J, Papanicolaou A, Yassour M,etal. De novo transcript sequence reconstruction from RNA-seq using the Trinity platform for reference generation and analysis. Nature Protocols, 2013, 8(8): 1494-1512.

[14] Altschul S F, Gish W, Miller W,etal. Basic local alignment search tool. Journal of Molecular Biology, 1990, 215(3): 403-410.

[15] Punta M, Coggill P C, Eberhardt R Y,etal. The Pfam protein families database. Nucleic Acids Research, 2012, 36: 263-266.

[16] Kanehisa M, Goto S, Sato Y,etal. KEGG for integration and interpretation of large-scale molecular data sets. Nucleic Acids Research, 2012, 40: D109-D114.

[17] Wu C H, Apweiler R, Bairoch A,etal. The universal protein resource (UniProt): an expanding universe of protein information. Nucleic Acids Research, 2006, 34: D187-D191.

[18] Ye J. WEGO: a web tool for plotting GO annotations. Nucleic Acids Research, 2006, 34(Web Server issue): W293-W297.

[19] Young M D, Wakefield M J, Smyth G K,etal. Gene ontology analysis for RNA-seq: accounting for selection bias. Genome Biology, 2010, 11(2): R14.

[20] Xie C, Mao X, Huang J,etal. KOBAS 2.0: a web server for annotation and identification of enriched pathways and diseases. Nucleic Acids Research, 2011, 39(Web Server issue): W316-W322.

[21] Langmead B, Trapnell C, Pop M,etal. Ultrafast and memory-efficient alignment of short DNA sequences to the human genome. Genome Biology, 2009, 10(3): R25.

[22] Li B, Dewey C N, Li B,etal. RSEM: accurate transcript quantification from RNA-Seq data with or without a reference genome. BMC Bioinformatics, 2011, 12(1): 323.

[23] Wang L, Feng Z, Wang X,etal. DEGseq: an R package for identifying differentially expressed genes from RNA-seq data. Bioinformatics (Oxford, England), 2010, 26(1): 136-138.

[24] Ernst J, Barjoseph Z. STEM: a tool for the analysis of short time series gene expression data. BMC Bioinformatics, 2006, 7(7): 1-11.

[25] Sun Y, Wang F, Wang N,etal. Transcriptome exploration inLeymuschinensisunder saline-alkaline treatment using 454 pyrosequencing. Plos One, 2013, 8(1): e53632.

[26] Zhao P, Liu P, Yuan G,etal. New insights on drought stress response by global investigation of gene expression changes in sheepgrass (Leymuschinensis). Frontiers in Plant Science, 2016, 7(147): 954.

[27] Devoto A, Ellis C. Expression profiling revealsCOI1 to be a key regulator of genes involved in wound- and methyl jasmonate-induced secondary metabolism, defence, and hormone interactions. Plant Molecular Biology, 2005, 58(4): 497-513.

[28] Turner J G, Ellis C, Devoto A. The jasmonate signal pathway. Plant Cell, 2002, 14(Suppl 1): S153-S164.

[29] Zhou W. Molecular Mechanism of Jasmonate-regulated Plant Defense Controled by JAV1. Beijing: Tsinghua University, 2013.

周武. JAV1調(diào)控茉莉素介導(dǎo)的植物抗性反應(yīng)的分子機(jī)理研究. 北京： 清華大學(xué), 2013.

[30] Yan J, Zhang C, Gu M,etal. TheArabidopsisCORONATINE INSENSITIVE1 protein is a jasmonate receptor. Plant Cell, 2009, 21(8): 2220.

[31] Adams E, Turner J. Illuminating COI1: a component of theArabidopsisjasomonate receptor complex also interacts with ethylene signaling. Plant Signaling & Behavior, 2010, 5(12): 1682-1684.

[32] Thines B, Katsir L, Melotto M,etal. JAZ repressor proteins are targets of the SCF(COI1) complex during jasmonate signalling. Nature, 2007, 448(7154): 661-665.

[33] Chen X Y, Tang Z C. Plant Physiology and Moleculer Biology. Beijing: Higher Education Press, 2012: 553-554.

陳曉亞, 湯章城. 植物生理與分子生物學(xué). 北京： 高等教育出版社, 2012: 553-554.

[34] Ma Y, Szostkiewicz I, Korte A,etal. Regulators of PP2C phosphatase activity function as abscisic acid sensors. Science, 2009, 324(5930): 1064.

[35] Park S Y, Fung P, Nishimura N,etal. Abscisic acid inhibits type 2C protein phosphatases via the PYR/PYL family of START proteins. Science, 2009, 324(5930): 1068-1071.

[36] Gonzalez-Guzman M, Pizzio G A, Antoni R,etal.ArabidopsisPYR/PYL/RCAR receptors play a major role in quantitative regulation of stomatal aperture and transcriptional response to abscisic acid. The Plant Cell, 2012, 24(6): 2483-2496.

[37] Schweighofer, Alois, Hirt,etal. Plant PP2C phosphatases: emerging functions in stress signaling. Trends in Plant Science, 2004, 9(5): 236.

[38] Jia L, Liu M M, Zhang H Q,etal. Antioxidant defense system responses ofArtemisiafrigidato mechanical damage. Journal of Zhejiang A & F University, 2016, 33(3): 462-470.

賈麗, 劉盟盟, 張洪芹, 等. 冷蒿抗氧化防御系統(tǒng)對(duì)機(jī)械損傷的響應(yīng). 浙江農(nóng)林大學(xué)學(xué)報(bào), 2016, 33(3): 462-470.

[39] Domergue F, Vishwanath S J, Joubès J,etal. Three Arabidopsis fatty acyl-coenzyme A reductases, FAR1, FAR4, and FAR5, generate primary fatty alcohols associated with suberin deposition. Plant Physiology, 2010, 153(4): 1539-1554.

[40] Koo A J K, Cooke T F, Howe G A. Cytochrome P450 CYP94B3 mediates catabolism and inactivation of the plant hormone jasmonoyl-L-isoleucine. Proceedings of the National Academy of Sciences of the United States of America, 2011, 108(22): 9298.

[41] Yan Y, Stolz S, Chételat A,etal. A downstream mediator in the growth repression limb of the jasmonate pathway. Plant Cell, 2007, 19(8): 2470-2483.

[42] Walley J W, Coughlan S, Hudson M E,etal. Mechanical stress induces biotic and abiotic stress responses via a novel cis-element. Plos Genetics, 2007, 3(10): 1800-1812.

[43] Chen F, D’Auria J C, Tholl D,etal. AnArabidopsisthalianagene for methylsalicylate biosynthesis, identified by a biochemical genomics approach, has a role in defense. Plant Journal, 2003, 36(5): 577-588.