不同刈割強度下羊草轉錄組研究

趙勁博,侯向陽*,武自念,任衛波,胡寧寧,郭豐輝,馬文靜(1.中國農業科學院草原研究所,內蒙古 呼和浩特 010010;2. 農業部草地生態與修復治理重點實驗室，內蒙古 呼和浩特 010010)

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

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

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

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

1 材料與方法

1.1 試驗材料

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

1.2 試驗設計

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

1.3 RNA提取及轉錄組測序

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

1.4 數據組裝與蛋白編碼框預測

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

1.5 組裝結果注釋與功能通路富集

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

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

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

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

2 結果與分析

2.1 羊草轉錄組數據組裝

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

圖1 轉錄本與Unigene長度分布Fig.1 Length distribution of transcripts and Unigenes橫坐標為轉錄本或Unigene的長度，縱坐標為其數目。 Abscissa for the length of transcripts or Unigenes, ordinate for their number.

2.2 Unigene注釋結果

通過對Unigene的蛋白預測，得到具有蛋白編碼框的Unigene共54555條，然后使用Blastx與Blastp將核酸及蛋白序列與數據庫進行比對，總共注釋到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分類結果主要與信號轉導機制(signal transduction mechanisms)、翻譯后修飾和蛋白反轉、分子伴侶(posttranslational modification, protein turnover, chaperones)等相關(圖2)。有關細胞、細胞組分、結合、催化、細胞過程、代謝過程等結果在GO分類中較為顯著(圖3)。

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

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

2.3 差異基因數目分析

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

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

圖4 所有差異基因stem聚類Fig.4 Stem clustering of all DEGs 曲線代表差異基因隨脅迫程度增加表達量變化趨勢，彩色框與無色框分布代表顯著性和非顯著性聚類，左上角與左下角數字分別為趨勢類型編號和聚類顯著性。 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 隨刈割強度增強表達量呈上升與下降趨勢差異基因熱圖Fig.5 Heatmap of DEGs having up and down trend expression with the increase of mowing intensity左圖呈上升趨勢，右圖呈下降趨勢，表達量均經過z-score標準化處理。 Left represents up trend expression, right represents down trend expression, expression are normalized by z-score.

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

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

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

FDR:錯誤發現率 False discovery rate.

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

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

表2 刈割下調差異基因GO富集結果(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 表達量隨刈割強度增強呈上升趨勢差異基因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 表達量隨刈割強度增強呈下降趨勢差異基因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 刈割差異基因脫落酸與茉莉酸信號轉導代謝通路富集結果Fig.6 ABA and JA signal transduction metabolic pathway enrichment of mowing DEGs紅色代表上調，綠色代表下調。Red represents up-regulated, green represents down-regulated.

3 討論

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

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

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

本研究針對羊草刈割強度設置了不同留茬高度作為處理實驗組，通過轉錄組測序得到大量的轉錄本數據。在常規分析的基礎下，重點對注釋到其他物種已知序列的關鍵基因進行了研究，而對于一些功能不明確或沒有注釋結果但表達量有規律性的基因，還有待在結構與功能等方面進行進一步的分析探索，從而使數據挖掘與利用更深入全面，以助于探求羊草在機械損傷逆境下最本質的機理。

4 結論

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

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.

李永宏. 內蒙古錫林河流域羊草草原和大針茅草原在放牧影響下的分異和趨同. 植物生態學報, 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.

趙鶴. 羊草對奶牛營養價值及其日糧組合效應的研究. 大慶： 黑龍江八一農墾大學, 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.

王志鋒, 王多伽, 于洪柱, 等. 刈割時間與留茬高度對羊草草甸草產量和品質的影響. 草業科學, 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.

仲延凱, 包青海. 不同刈割強度對天然割草地的影響. 中國草地, 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.

韓龍, 郭彥軍, 韓建國, 等. 不同刈割強度下羊草草甸草原生物量與植物群落多樣性研究. 草業學報, 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.

侯向陽. 草原植物基礎生物學研究進展與展望. 中國基礎科學, 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調控茉莉素介導的植物抗性反應的分子機理研究. 北京： 清華大學, 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.

陳曉亞, 湯章城. 植物生理與分子生物學. 北京： 高等教育出版社, 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.

賈麗, 劉盟盟, 張洪芹, 等. 冷蒿抗氧化防御系統對機械損傷的響應. 浙江農林大學學報, 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.