二倍體西瓜及其同源四倍體葉片sRNA表達譜分析

龍婭麗 江雪飛 周鵬 徐子健 孫夢利

摘 要 植物多倍化后的基因表達平衡重建受到sRNA調節(jié)，miRNA是植物多倍化中基因表達進行反式調控的主要小分子成員。本研究利用新一代高通量測序對西瓜二倍體（2n）及其同源四倍體（4n）葉片的sRNA進行表達譜測序分析。結果表明：（1）2n中sRNA以21 nt為主，4n中則以24 nt為主。（2）在2n和4n鑒定出已知miRNA分別為110個和172個，新miRNA分別為246和829個，其中差異表達的已知miRNA 205個，新miRNA 815個。（3）GO功能顯著性富集分析表明差異表達的miRNA在細胞組成、代謝過程及催化活性上顯著富集。（4）KEGG代謝分析表明，差異表達的已知miRNA主要富集在初級代謝和次級代謝途徑上，新miRNA則主要富集在植物激素信號轉導途徑、RNA轉運途徑和ABC轉座子途徑上。（5）差異表達的miRNA中與植物抗逆脅迫相關的已知miRNA 25個，新miRNA 62個，且80.46%抗逆miRNA在西瓜二倍體中上調表達，是四倍體中的4倍以上。根據(jù)miRNA的負向調控原理推測這些上調表達的miRNA在西瓜二倍體中可能會抑制部分抗逆基因的表達，影響二倍體的抗逆性。本研究在轉錄后調控調控水平上為西瓜多倍體較二倍體擁有更強抗逆性提供了理論依據(jù)，同時對西瓜多倍體及二倍體之間的基因表達譜分析做了補充。

關鍵詞 miRNA；西瓜；多倍體；抗逆性

中圖分類號 S651 文獻標識碼 A

Abstract After plant polyploidy， the regulation of gene expression balance is regulated by sRNA， and miRNA is the major small molecule in trans-regulation of gene expression during plant polyploidy. In this study， the expression profiles of sRNA in watermelon diploid （2n） and its tetraploid （4N） leaves were sequenced by using a new generation of high-throughput sequencing. Sequencing results showed that： （1） In 2n， sRNA was dominated by 21 nt， while sRNA was dominated by 24 nt in 4n. （2） In 2n and 4n， the number of identified miRNA was 110 and 172， and the new miRNA was 246 and 829 respectively， while the number of the known miRNA in differential expression was 205， and the number of new miRNA was 815. （3） GO functional enrichment analysis showed that the differential expression miRNA was significantly enriched in cell composition， metabolic process and catalytic activity. （4） KEGG metabolic analysis showed that differential expression known miRNA was mainly enriched in primary metabolism and secondary metabolic pathways. New miRNA was mainly concentrated on plant hormone signal transduction pathway， RNA transport pathway and ABC transposon pathway. （5） Among the differential expression miRNA， the number of known miRNA and new miRNA related to plant stress was respectively 25 and 62， and 80.46% resistant miRNA was down regulated in tetraploid watermelon. Only 19.54% of the up-regulated expression in diploid watermelon was expressed in miRNA. According to the principle of miRNA's negative regulation， these miRNA could inhibit the expression of some anti stress genes in diploid watermelon. Therefore， this study would provide a basis for the study of watermelon polyploid with stronger resistance than diploid in post transcriptional regulation， and the gene expression profiles between polyploid and diploid watermelon will be supplemented.

Key words miRNA； watermelon； tetraploid； stress resistance

doi 10.3969/j.issn.1000-2561.2018.04.008

西瓜（Citrullus lahatus），尤其是利用二倍體和四倍體雜交獲得的三倍體無籽西瓜，因食用方便、果實大和品質好而備受消費者歡迎[1-2]。多倍體西瓜不僅品質好，且多數(shù)情況下對環(huán)境的適應性優(yōu)于其二倍體。2000年郭啟高等[3]進行西瓜試管苗快繁研究時發(fā)現(xiàn)，四倍體西瓜對營養(yǎng)逆境有更強的適應性。2003年劉文革等[4]研究發(fā)現(xiàn)多倍體西瓜的耐鹽性顯著強于二倍體，同時還發(fā)現(xiàn)通過預冷處理后西瓜多倍體的耐低溫能力強于二倍體；對不同倍性西瓜進行淹水脅迫處理發(fā)現(xiàn)：西瓜幼苗的耐淹能力同樣表現(xiàn)為多倍體強于二倍體。2009年劉文革等[5]對18個同基因型不同倍性西瓜品種的苗期枯萎病抗性進行接種鑒定，結果表明，同基因型的西瓜二倍體植株最先進入發(fā)病期，同源多倍體西瓜對枯萎病的抗性優(yōu)于同源二倍體。2012年朱紅菊等[6]再次證明西瓜多倍體耐鹽性強于二倍體。因此清楚了解多倍體和二倍體之間抗逆調節(jié)分子機制的異同，有利于西瓜多倍體育種獲得優(yōu)質、高抗、高產的西瓜。

植物多倍化后的基因表達平衡重建受到sRNA調節(jié)[7-8]。sRNA作為一種非編碼RNA的產物以多種形式參與多倍體基因表達調控的多個方面，成為多倍化中基因表達的重要調節(jié)因子。miRNA是植物多倍化中基因表達進行反式調控的主要小分子成員，由其介導的基因沉默是轉錄后水平基因調控的重要途徑之一[7-9]。由于miRNA往往位于基因調控上游，其表達量變化將影響到相應的靶基因，可以認為是導致植物多倍化后表型變化主要原因之一，植物多倍化中這種負調控抑制是多倍體后代生長活力和適應性增強的重要分子基礎[7]。西瓜也不例外，大量的研究表明西瓜多倍體在多數(shù)情況下的抗逆性要優(yōu)于二倍體，但是，西瓜多倍體和二倍體之間轉錄后調控水平上的差異如何，目前還未見到相關報道。因此，以抗逆性差異較大的西瓜二倍體及其同源四倍體葉片為材料，采用二代深度測序的方法，對西瓜二、四倍體sRNA表達譜進行分析，明確他們之間miRNA表達量的差異，這有利于清楚了解西瓜多倍體和二倍體之間抗逆調節(jié)分子機制的異同，為西瓜多倍體抗逆育種提供轉錄后調控水平上的分子依據(jù)。

1 材料與方法

1.1 材料

由中國熱帶農業(yè)科學院熱帶作物品種資源研究所提供的二倍體（FR-32-1B-2n）和同源四倍體（FR-32-1B-4n）西瓜為材料，其中同源四倍體由二倍體通過秋水仙素誘導獲得。該西瓜材料種植在中國熱帶農業(yè)科學院熱帶作物品種資源研究所的日光溫室內，正常管理。

1.2 方法

1.2.1 取樣方法 待西瓜長到坐果期時，分別采集2～3個植株的正常功能葉片3～5片并混合作為1個生物學樣品，經液氮速凍后置于超低溫（-80 ℃）冰箱里保存?zhèn)溆谩?/p>

1.2.2 樣品RNA提取、文庫構建及測序 按照美國AXYGEN公司RNA小量制備試劑盒內植物葉片組織RNA提取方法提取樣品總RNA，不同長度大小的片段用PAGE膠分離，將長度大小在18～30 nt之間的片段進行回收。設計引物進行PCR擴增，完成sRNA文庫構建，該文庫的質量和產量經檢測合格后在Illumina Hiseq 2000平臺上進行測序。

1.3 數(shù)據(jù)處理

1.3.1 sRNA數(shù)據(jù)分析 對測序獲得的原始序列進行一定的處理，如去除低質量的序列，去掉接頭或者去掉污染等，然后獲得干凈序列。之后將獲得的sRNA 干凈序列比對到GenBank和Rfam數(shù)據(jù)庫進行注釋，然后通過Basic Local Alignment Search Tool （BLAST）移除rRNA、tRNA、snRNA和snoRNA。隨后，所有的干凈序列被比對到西瓜基因組（ftp：//www.icugi.org/pub/genome/watermelon/97103/v1/），利用Short Oligo Alignment Program （SOAP） （http：//soap.genomics.org.cn）[10]計算出表達量和分布情況。接著，所有的clean reads比對到Rfam1 2.1 and miRBase 21.0數(shù)據(jù)庫，使用BGI軟件tag 2進行進一步注釋和分類[11]。然而，一些sRNA reads可能會有多個注釋分布，為了保證每一個特異性sRNA序列只能有一個注釋，遵循以下規(guī)律：rRNA etc（Genbank > Rfam）>保守miRNA>重復片段>外顯子>內含子[8]。

1.3.2 miRNAs鑒定及靶基因預測 為了鑒定保守的miRNAs，所有的特異性序列被比對到miRBase 21.0 [12]。另外，使用Mireap 程序（https：//sourceforge.net/projects/mireap）對帶有莖環(huán)結構的未知miRNAs進行鑒定。具體鑒定方法參考相關文獻利用TargetFinder[13-15] （https：//github.com/carringtonlab/TargetFinder）和psRobot program [16]（http：//omicslab. genetics.ac.cn/psRobot）工具對milRNAs的靶基因進行預測，具體步驟參見Allen等[17]和Schwab等[18]方法。靶基因功能注釋參見Gong等[19]方法。

1.3.3 差異表達miRNAs分析 利用DEGseq R包對miRNAs差異表達進行分析[20]， miRNAs表達水平使用RPM（reads per million reads）值進行計算，用FDR （False Discovery Rate）值進行校準[21]。根據(jù)Kang等[22]的方法，F(xiàn)DR ≤ 0.001、P value ≤ 0.001和|log2 RPM Ratio| ≥1定義為顯著差異表達的miRNA。

2 結果與分析

2.1 西瓜二倍體及四倍體中sRNA數(shù)據(jù)質量及長度分布

西瓜二倍體和同源四倍體葉片的小RNA測序文庫，分別命名為2n和4n。2個文庫經IlluminaHiSeqTM2000測序分別獲得的Raw reads 12 254 705條和11 846 305條。去除低質量序列后，獲得Clean reads分別均為11 984 042條和11 600 832條，分別占Raw reads的97.98%和98.19% （表1）。表明西瓜二倍體和同源四倍體sRNA文庫構建及測序數(shù)據(jù)質量較高，且所得干凈序列的數(shù)量滿足試驗需要。每個樣本所比對上的片段長度為20～24 nt，2n中尤以21 nt的片段最多，4n中則以24 nt的為主（圖1），說明西瓜倍性發(fā)生變化后，其體內的sRNA長度分布也受到了影響、發(fā)生了改變。

2.2 西瓜二倍體及四倍體中公共及特有序列統(tǒng)計分析

比較2n和4n兩個數(shù)據(jù)庫發(fā)現(xiàn)，兩個數(shù)據(jù)庫中共有的特異性sRNAs僅有7.53% （226 004），但是2個數(shù)據(jù)庫中共有的總sRNA卻有很大的比例（70.04%，16 519 602），此外兩個數(shù)據(jù)庫中特有的特異性sRNAs超過了40%，但是特有的總sRNA卻少于20% （圖2）。說明兩個數(shù)據(jù)庫中出現(xiàn)了大量的特有序列，這可能是倍性變化造成的。

2.3 sRNA干凈序列與西瓜基因組比對

所有的總sRNAs和特異性sRNAs序列被比對到西瓜基因組，僅有序列完全匹配才可進行下一步分析。由于基因組之間存在的差異，結果發(fā)現(xiàn)僅有一部分序列能完全比對上，2個庫中的總sRNAs比對上的序列都在80%以上，而特異性sRNAs比對的情況兩個庫相差比較大，2n的僅有44.33%完全比對上，4n比對上的序列則高達74.48%（表2）。sRNAs片段在各條染色體上的分布情況顯示：無論是在2n數(shù)據(jù)庫還是4n數(shù)據(jù)庫中，正鏈上片段數(shù)分布最多的都在1號染色體和4號染色體上，負鏈上片段數(shù)分布最多的則是11號染色體（圖3）。

2.4 西瓜二倍體及四倍體中 sRNA分類注釋

在深度測序小RNA數(shù)據(jù)集中，有大量不同類型的無編碼RNAs （ncRNA）。經過不同數(shù)據(jù)庫的比對后獲得核糖體RNA （rRNA）、轉移RNA （tRNA）、microRNA （miRNA）、小核RNA （snRNA）及核仁小分子RNA （snoRNA）（圖4）。在兩個數(shù)據(jù)庫中的總sRNA序列比例最大的是未注釋的sRNAs。在2n和4n數(shù)據(jù)庫中分別約有44.15% （5 291 139）和30.82% （3 575 214）的序列被注釋到RNAs （ncRNAs），包括miRNAs、rRNA、tRNAs、snoRNAs和snRNAs。在所有的sRNA中，miRNA序列在兩個數(shù)據(jù)庫中的比例分別占11.07%和31.42%。在兩個數(shù)據(jù)庫中的特異性sRNA序列中情況類似，但未注釋到的sRNA比例更高，注釋到的rRNA數(shù)量減幅最大，最明顯。在4n數(shù)據(jù)庫中，注釋到的sRNA序列數(shù)量大幅度下降，下降比例高于2n數(shù)據(jù)庫（圖4）。

2.5 已知miRNA與新miRNA鑒定及差異分析

通過blast或bowtie將sRNA和miRBase數(shù)據(jù)庫比對，鑒定出已知miRNA。另外，對未注釋上任何RNA且比對上基因組的外顯子反義鏈、內含子、基因間區(qū)的sRNAs，通過軟件Mireap鑒定新的miRNA。在2n和4n兩個數(shù)據(jù)庫中分別鑒定出已知miRNA為110個和172個，新miRNA為246和829個。對2n和4n數(shù)據(jù)庫中的已知miRNA和新miRNA差異分析獲得205個差異表達的已知miRNA和815個新miRNA。表達量分析顯示差異表達的已知miRNA表達量在2n數(shù)據(jù)庫中高于4n，而新miRNA表達量在2個庫中相差不大，說明miRNA的調控具有物種倍性特異性。

2.6 差異表達的已知miRNA及新miRNA的GO功能分析

根據(jù)miRNA與其靶基因間的對應關系，對miRNA的靶基因的集合分別進行GO （Gene Ontology）富集分析。結果顯示無論是已知miRNA還是新miRNA，GO功能都分為生物過程、細胞組分和分子功能3類，在生物過程中，miRNA靶基因主要集中在“細胞過程”、“代謝過程”和“信號過程”上；在細胞組分中，miRNA靶基因主要集中在“細胞”、“細胞器”和“代謝”上；在分子功能中則主要集中在“結合”和“催化活性”上（圖6～7）。說明miRNA介導調控的生物過程較廣泛。

2.7 差異表達已知miRNA及新miRNA的KEGG代謝通路分析

在生物體內，不同基因相互協(xié)調行使其生物學功能，通過Pathway顯著性富集分析，能確定候選靶基因參與的最主要生化代謝途徑和信號轉導途徑，助于更進一步了解基因的生物學功能。KEGG富集結果顯示，已知miRNA主要富集在初級代謝和次級代謝途徑上，新miRNA則主要富集在植物激素信號轉導途徑、RNA轉運途徑和ABC轉座子途徑上。這說明已知和新的miRNA在代謝途徑上的差異較大，他們分別在不同的代謝途徑上發(fā)揮著重要的調節(jié)作用。

2.8 抗逆差異表達miRNA分析

經過仔細查閱文獻核對，在西瓜二倍體和四倍體sRNA數(shù)據(jù)庫中篩選獲得25個和植物抗逆脅迫相關的已知miRNA，62個和植物抗逆脅迫相關新miRNA。對于已知的miRNA，有9個在2n中上調表達，16個下調表達；而對于新miRNA則只有1個在2n中下調表達，其余61個都在2n中上調表達（表3）。根據(jù)miRNA負調控靶基因的原理推測，大部分在西瓜中上調表達的miRNA會抑制靶基因在二倍體中的表達。相反，這些miRNA會促進靶基因在西瓜四倍體中的表達。這從轉錄調控水平上證明了西瓜四倍體有較強的抗逆性。

3 討論

對西瓜二倍體和同源四倍體葉片的高通量測序分別得到了11 984 042和11 600 832條clean reads，將這些clean reads與西瓜參考基因組進行對比分類注釋，結果顯示在西瓜二倍體（2n）和同源四倍體（4n）兩個數(shù)據(jù)庫中，種類和數(shù)量最多的為rRNA與miRNA，這與冬小麥基因組測序結果相同[23]，說明sRNA在植物進化中具有一定的保守性。2n和4n兩個數(shù)據(jù)庫中分別鑒定出已知miRNA為110個和172個，新miRNA為246和829個。明顯可以看出新miRNA的數(shù)量高于已知miRNA的數(shù)量。已有的研究顯示已知miRNA在進化上相對比較保守，一般只負責調控植物的生長發(fā)育等基礎代謝，在嚴酷的進化條件下保持著序列的保守性，而新miRNA多是隨著環(huán)境變化，特異性進化出的序列，參與植物的抗逆等特異性調控[24]。因此，本研究所鑒定出來的新miRNA對于西瓜屬的特異性調節(jié)具有重要的參考意義。不同植物間sRNA的分布差異較大，本研究的西瓜2n中以21 nt的片段最多，在楊樹、大豆和番茄中，21 nt的sRNA分布也為最多[25-27]；擬南芥、小麥和棉花中，24 nt的sRNA分布最多[28-30]，本研究的西瓜4n中以24 nt的sRNA分布為最多，說明在不同的植物中，sRNA特異性分布表達，即便是同一個品種不同倍性的植物sRNA也表現(xiàn)出了特異性分布表達的特點，表明植物的加倍對sRNA有很大的影響。

miRNA往往會通過表達量差異進而影響基因功能和代謝途徑上的差異，這對植物的多倍體進化起著重要的調控作用。對2個庫中差異表達的miRNA分析顯示，已知miRNA表達量在數(shù)據(jù)庫2n中高于4n，根據(jù)miRNA的反向調控作用原理分析，推測在西瓜二倍體中有更多的差異表達基因被抑制表達或是沉默了。miRNA的表達差異會對細胞內KEGG通路和GO生物過程產生重要影響[31]。本實驗結果顯示，差異表達的miRNA靶基因在生物過程中主要集中在“細胞過程”、“代謝過程”和“信號過程”上。差異表達的miRNA靶基因的KEGG富集結果顯示已知miRNA主要富集在初級代謝和次級代謝途徑上，新miRNA則主要富集在植物激素信號轉導途徑、RNA轉運途徑和ABC轉座子途徑上，這進一步說明了已知miRNA在進化上的保守性，新miRNA在植物進化上的特異性。

植物多倍體的抗逆性往往優(yōu)于二倍體，這離不開抗逆差異表達miRNA的逆調控作用，因為miRNA控制的基因往往參與重要的抗逆和發(fā)育過程[32]。本實驗結果分析顯示，有80%的miRNA在西瓜二倍體中上調表達，根據(jù)miRNA對靶基因的調節(jié)特點可知，這些miRNA調節(jié)的抗逆靶基因在二倍體中的表達量會低于四倍體，甚至有可能被完全沉默掉，這會進一步導致西瓜四倍體中高度表達的抗逆基因會比二倍體中多。這在其他物種的研究中已經被證實，如二倍體泡桐中miR156b、miR28a/b等miRNAs表達量較高，相反其調控的靶基因表達量在四倍體泡桐中較高，這就導致了四倍體泡桐的抗鹽性比二倍體泡桐要強[33]。因此，我們的結果為四倍體優(yōu)于二倍體的抗逆性提供了轉錄調控水平上的證據(jù)。

參考文獻

[1] 譚素英，黃秀強，劉文革. 三倍體無籽西瓜的優(yōu)越性及系列無籽西瓜新品種[J]. 中國西瓜甜瓜，1994（4）：22-23.

[2] 劉文革，何 楠，趙勝杰，等. 我國三倍體無籽西瓜的周年生產與栽培[J]. 長江蔬菜，2014（14）：1-6.

[3] 郭啟高，宋 明，楊天秀，等. 西瓜子葉組織培養(yǎng)中四倍體的產生及鑒定[J]. 西南大學學報（自然科學版），2000，22（4）：298-300.

[4] 劉文革. 不同染色體倍性西瓜（Citrullus lanatus）的遺傳變異和抗逆機理研究[D]. 西安：西北農林科技大學，2003.

[5] 劉文革，閻志紅，趙勝杰，等. 不同染色體倍性西瓜對枯萎病的抗性研究[J]. 長江蔬菜，2009，（18）：19-20.

[6] 朱紅菊. NaCl脅迫對不同倍性西瓜幼苗影響及DNA甲基化分析[D]. 北京：中國農業(yè)科學院，2012.

[7] Ha M，Lu J，Tian L，et al. Small RNAs serve as a genetic buffer against genomic shock in Arabidopsis interspecific hybrids and allopolyploids[J]. Proceedings of the National Academy of Sciences，2009，106（42）：17 835-17 840.

[8] Ng D W，Zhang C，Miller M，et al. cis- and trans-Regulation of miR163 and target genes confers natural variation of secondary metabolites in two Arabidopsis species and their allopolyploids[J]. Plant Cell，2011，23（5）：1 729-1 740.

[9] 趙旭博，李愛麗，毛 龍. 植物多倍化過程中小分子RNA調控基因表達機制研究進展[J]. 作物學報，2013，39（8）：1 331-1 338.

[10] Li R，Li Y，Kristiansen K，et al. SOAP： short oligonucleotide alignment program[J]. Bioinformatics，2008，24（5）：713-714.

[11] Yuan C，Wang X，Geng R，et al. Discovery of cashmere goat （Capra hircus） microRNAs in skin and hair follicles by Solexa sequencing[J]. BMC Genomics，2013，14（1）：511.

[12] Kozomara A，Griffiths-Jones S. miRBase：annotating high confidence microRNAs using deep sequencing data[J]. Nucleic Acids Res，2014，42（Database issue）：D68-D73.

[13] Zuker M. Mfold web server for nucleic acid folding and hybridization prediction[J]. Nucleic Acids Research，2003，31（13）：3406-3415.

[14] Fahlgren N，Howell M D，Kasschau K D，et al. High-throughput sequencing of Arabidopsis microRNAs：evidence for freq uent birth and death of MIRNA genes[J]. PLoS ONE，2007，2（2）：e219.

[15] Fahlgren N，Carrington J C. miRNA Target Prediction in Plants[J]. Methods in Molecular Biology，2009，592：51-57.

[16] Wu HJ，Ma Y K，Chen T，et al. PsRobot： a web-based plant small RNA meta-analysis toolbox[J]. Nucleic Acids Research，2012，40（Web Server issue）：W22-W28.

[17] Allen E，Xie Z，Gustafson A M，et al. microRNA-directed phasing during trans-acting siRNA biogenesis in plants[J]. Cell，2005，121（2）：207-221.

[18] Schwab R，Palatnik J F，Riester M，et al. Specific effects of microRNAs on the plant transcriptome[J]. Developmental Cell，2005，8（4）：517-527.

[19] Gong S，Hao Z，Meng J，et al. Digital gene expression analysis to screen disease resistance-relevant genes from leaves of herbaceous peony （Paeonia lactiflora Pall.） infected by Botrytis cinerea [J]. PLoS ONE，2015，10（7）：e0133305.

[20] Wang L，F(xiàn)eng Z，Wang X，et al. DEGseq： an R package for identifying differentially expressed genes from RNA-seq data[J]. Bioinformatics，2010，26（1）：136-138.

[21] Benjamini Y，Hochberg Y. Controlling the false discovery rate： A practical and powerful approach to multiple testing[J]. Journal of the Royal Statistical Society，1995，57（1）：289-300.

[22] Kang K，Zhong J，Jiang L，et al. Identification of microRNA-Like RNAs in the filamentous fungus Trichoderma reesei by solexa sequencing[J]. PLoS ONE，2013，8（10）：e76288.

[23] 黃 儒，蒼 晶，于 晶，等. 冬小麥小RNA高通量測序及生物信息學分析[J]. 植物學報，2014，49（1）：8-18.

[24] 南文智，吳 嫚，于霽雯，等. 利用高通量測序技術鑒定棉纖維發(fā)育相關miRNAs及其靶基因[J]. 棉花學報，2013，25（4）：300-308.

[25] Moxon S，Jing R，Szittya G，et al. Deep sequencing of tomato short RNAs identifies microRNAs targeting genes involved in fruit ripening[J]. Genome Research，2008，18（10）：1 602-1 609.

[26] Li B S，Qin Y R，Duan H，et al. Genome-wide characterization of new and drought stress responsive microRNAs in Populus euphratica[J]. Journal of Experimental Botany，2011，62（11）：3 765-3 779.

[27] Turner M，Yu O，Subramanian S. Genome organization and characteristics of soybean microRNAs[J]. BMC Genomics，2012，13：169.

[28] Kasschau K D，F(xiàn)ahlgren N，Chapman E J，et al. Genome-wide profiling and analysis of Arabidopsis siRNAs[J]. PLoS Biology，2007，5（3）：e57.

[29] Wei B，Cai T，Zhang R，et al. Novel microRNAs uncovered by deep sequencing of small RNA transcriptomes in bread wheat （Triticum aestivum L.） and Brachypodium distachyon （L.） Beauv[J]. Functional and Integrative Genomics，2009，9（4）：499-511.

[30] Wang Z M，Xue W，Dong C J，et al. A comparative miRNAome analysis reveals seven fiber initiation-related and 36 novel miRNAs in developing cotton ovules[J]. Molecular Plant，2012，5（4）：889-900.

[31] Tang M，Mao D，Xu L，et al. Integrated analysis of miRNA and mRNA expression profiles in response to Cd exposure in rice seedlings[J]. BMC Genomics，2014，15（1）：835.

[32] Khraiwesh B，Zhu J K，Zhu J. Role of miRNAs and siRNAs in biotic and abiotic stress responses of plants[J]. Biochim Biophys Acta，2012，1 819（2）：137-148.

[33] Fan G，Li X，Deng M，et al. Comparative analysis and Identification of miRNAs and their target genes responsive to salt stress in diploid and tetraploid paulownia fortunei seedlings[J]. PLoS ONE，2016，11（2）：e0149617.