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

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

關(guān)鍵詞 miRNA；西瓜；多倍體；抗逆性

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

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

1 材料與方法

1.1 材料

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

1.2 方法

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

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

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

1.3.1 sRNA數(shù)據(jù)分析 對測序獲得的原始序列進(jìn)行一定的處理，如去除低質(zhì)量的序列，去掉接頭或者去掉污染等，然后獲得干凈序列。之后將獲得的sRNA 干凈序列比對到GenBank和Rfam數(shù)據(jù)庫進(jìn)行注釋，然后通過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]計(jì)算出表達(dá)量和分布情況。接著，所有的clean reads比對到Rfam1 2.1 and miRBase 21.0數(shù)據(jù)庫，使用BGI軟件tag 2進(jìn)行進(jìn)一步注釋和分類[11]。然而，一些sRNA reads可能會(huì)有多個(gè)注釋分布，為了保證每一個(gè)特異性sRNA序列只能有一個(gè)注釋，遵循以下規(guī)律：rRNA etc（Genbank > Rfam）>保守miRNA>重復(fù)片段>外顯子>內(nèi)含子[8]。

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

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

2 結(jié)果與分析

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

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

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

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

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

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

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

在深度測序小RNA數(shù)據(jù)集中，有大量不同類型的無編碼RNAs （ncRNA）。經(jīng)過不同數(shù)據(jù)庫的比對后獲得核糖體RNA （rRNA）、轉(zhuǎn)移RNA （tRNA）、microRNA （miRNA）、小核RNA （snRNA）及核仁小分子RNA （snoRNA）（圖4）。在兩個(gè)數(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序列在兩個(gè)數(shù)據(jù)庫中的比例分別占11.07%和31.42%。在兩個(gè)數(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且比對上基因組的外顯子反義鏈、內(nèi)含子、基因間區(qū)的sRNAs，通過軟件Mireap鑒定新的miRNA。在2n和4n兩個(gè)數(shù)據(jù)庫中分別鑒定出已知miRNA為110個(gè)和172個(gè)，新miRNA為246和829個(gè)。對2n和4n數(shù)據(jù)庫中的已知miRNA和新miRNA差異分析獲得205個(gè)差異表達(dá)的已知miRNA和815個(gè)新miRNA。表達(dá)量分析顯示差異表達(dá)的已知miRNA表達(dá)量在2n數(shù)據(jù)庫中高于4n，而新miRNA表達(dá)量在2個(gè)庫中相差不大，說明miRNA的調(diào)控具有物種倍性特異性。

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

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

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

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

2.8 抗逆差異表達(dá)miRNA分析

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

3 討論

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

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

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

參考文獻(xiàn)

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

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

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

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

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

[6] 朱紅菊. NaCl脅迫對不同倍性西瓜幼苗影響及DNA甲基化分析[D]. 北京：中國農(nóng)業(yè)科學(xué)院，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調(diào)控基因表達(dá)機(jī)制研究進(jìn)展[J]. 作物學(xué)報(bào)，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高通量測序及生物信息學(xué)分析[J]. 植物學(xué)報(bào)，2014，49（1）：8-18.

[24] 南文智，吳 嫚，于霽雯，等. 利用高通量測序技術(shù)鑒定棉纖維發(fā)育相關(guān)miRNAs及其靶基因[J]. 棉花學(xué)報(bào)，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.