二倍体西瓜及其同源四倍体叶片sRNA表达谱分析

龙娅丽 江雪飞 周鹏 徐子健 孙梦利

摘 要 植物多倍化后的基因表達平衡重建受到sRNA调节，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倍以上。根据miRNA的负向调控原理推测这些上调表达的miRNA在西瓜二倍体中可能会抑制部分抗逆基因的表达，影响二倍体的抗逆性。本研究在转录后调控调控水平上为西瓜多倍体较二倍体拥有更强抗逆性提供了理论依据，同时对西瓜多倍体及二倍体之间的基因表达谱分析做了补充。

关键词 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]。多倍体西瓜不仅品质好，且多数情况下对环境的适应性优于其二倍体。2000年郭启高等[3]进行西瓜试管苗快繁研究时发现，四倍体西瓜对营养逆境有更强的适应性。2003年刘文革等[4]研究发现多倍体西瓜的耐盐性显著强于二倍体，同时还发现通过预冷处理后西瓜多倍体的耐低温能力强于二倍体；对不同倍性西瓜进行淹水胁迫处理发现：西瓜幼苗的耐淹能力同样表现为多倍体强于二倍体。2009年刘文革等[5]对18个同基因型不同倍性西瓜品种的苗期枯萎病抗性进行接种鉴定，结果表明，同基因型的西瓜二倍体植株最先进入发病期，同源多倍体西瓜对枯萎病的抗性优于同源二倍体。2012年朱红菊等[6]再次证明西瓜多倍体耐盐性强于二倍体。因此清楚了解多倍体和二倍体之间抗逆调节分子机制的异同，有利于西瓜多倍体育种获得优质、高抗、高产的西瓜。

植物多倍化后的基因表达平衡重建受到sRNA调节[7-8]。sRNA作为一种非编码RNA的产物以多种形式参与多倍体基因表达调控的多个方面，成为多倍化中基因表达的重要调节因子。miRNA是植物多倍化中基因表达进行反式调控的主要小分子成员，由其介导的基因沉默是转录后水平基因调控的重要途径之一[7-9]。由于miRNA往往位于基因调控上游，其表达量变化将影响到相应的靶基因，可以认为是导致植物多倍化后表型变化主要原因之一，植物多倍化中这种负调控抑制是多倍体后代生长活力和适应性增强的重要分子基础[7]。西瓜也不例外，大量的研究表明西瓜多倍体在多数情况下的抗逆性要优于二倍体，但是，西瓜多倍体和二倍体之间转录后调控水平上的差异如何，目前还未见到相关报道。因此，以抗逆性差异较大的西瓜二倍体及其同源四倍体叶片为材料，采用二代深度测序的方法，对西瓜二、四倍体sRNA表达谱进行分析，明确他们之间miRNA表达量的差异，这有利于清楚了解西瓜多倍体和二倍体之间抗逆调节分子机制的异同，为西瓜多倍体抗逆育种提供转录后调控水平上的分子依据。

1 材料与方法

1.1 材料

由中国热带农业科学院热带作物品种资源研究所提供的二倍体（FR-32-1B-2n）和同源四倍体（FR-32-1B-4n）西瓜为材料，其中同源四倍体由二倍体通过秋水仙素诱导获得。该西瓜材料种植在中国热带农业科学院热带作物品种资源研究所的日光温室内，正常管理。

1.2 方法

1.2.1 取样方法 待西瓜长到坐果期时，分别采集2～3个植株的正常功能叶片3～5片并混合作为1个生物学样品，经液氮速冻后置于超低温（-80 ℃）冰箱里保存备用。

1.2.2 样品RNA提取、文库构建及测序 按照美国AXYGEN公司RNA小量制备试剂盒内植物叶片组织RNA提取方法提取样品总RNA，不同长度大小的片段用PAGE胶分离，将长度大小在18～30 nt之间的片段进行回收。设计引物进行PCR扩增，完成sRNA文库构建，该文库的质量和产量经检测合格后在Illumina Hiseq 2000平台上进行测序。

1.3 数据处理

1.3.1 sRNA数据分析 对测序获得的原始序列进行一定的处理，如去除低质量的序列，去掉接头或者去掉污染等，然后获得干净序列。之后将获得的sRNA 干净序列比对到GenBank和Rfam数据库进行注释，然后通过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数据库，使用BGI软件tag 2进行进一步注释和分类[11]。然而，一些sRNA reads可能会有多个注释分布，为了保证每一个特异性sRNA序列只能有一个注释，遵循以下规律：rRNA etc（Genbank > Rfam）>保守miRNA>重复片段>外显子>内含子[8]。

1.3.2 miRNAs鉴定及靶基因预测 为了鉴定保守的miRNAs，所有的特异性序列被比对到miRBase 21.0 [12]。另外，使用Mireap 程序（https：//sourceforge.net/projects/mireap）对带有茎环结构的未知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]。根据Kang等[22]的方法，FDR ≤ 0.001、P value ≤ 0.001和|log2 RPM Ratio| ≥1定义为显著差异表达的miRNA。

2 结果与分析

2.1 西瓜二倍体及四倍体中sRNA数据质量及长度分布

西瓜二倍体和同源四倍体叶片的小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文库构建及测序数据质量较高，且所得干净序列的数量满足试验需要。每个样本所比对上的片段长度为20～24 nt，2n中尤以21 nt的片段最多，4n中则以24 nt的为主（图1），说明西瓜倍性发生变化后，其体内的sRNA长度分布也受到了影响、发生了改变。

2.2 西瓜二倍体及四倍体中公共及特有序列统计分析

比较2n和4n两个数据库发现，两个数据库中共有的特异性sRNAs仅有7.53% （226 004），但是2个数据库中共有的总sRNA却有很大的比例（70.04%，16 519 602），此外两个数据库中特有的特异性sRNAs超过了40%，但是特有的总sRNA却少于20% （图2）。说明两个数据库中出现了大量的特有序列，这可能是倍性变化造成的。

2.3 sRNA干净序列与西瓜基因组比对

所有的总sRNAs和特异性sRNAs序列被比对到西瓜基因组，仅有序列完全匹配才可进行下一步分析。由于基因组之间存在的差异，结果发现仅有一部分序列能完全比对上，2个库中的总sRNAs比对上的序列都在80%以上，而特异性sRNAs比对的情况两个库相差比较大，2n的仅有44.33%完全比对上，4n比对上的序列则高达74.48%（表2）。sRNAs片段在各条染色体上的分布情况显示：无论是在2n数据库还是4n数据库中，正链上片段数分布最多的都在1号染色体和4号染色体上，负链上片段数分布最多的则是11号染色体（图3）。

2.4 西瓜二倍体及四倍体中 sRNA分类注释

在深度测序小RNA数据集中，有大量不同类型的无编码RNAs （ncRNA）。经过不同数据库的比对后获得核糖体RNA （rRNA）、转移RNA （tRNA）、microRNA （miRNA）、小核RNA （snRNA）及核仁小分子RNA （snoRNA）（图4）。在两个数据库中的总sRNA序列比例最大的是未注释的sRNAs。在2n和4n数据库中分别约有44.15% （5 291 139）和30.82% （3 575 214）的序列被注释到RNAs （ncRNAs），包括miRNAs、rRNA、tRNAs、snoRNAs和snRNAs。在所有的sRNA中，miRNA序列在两个数据库中的比例分别占11.07%和31.42%。在两个数据库中的特异性sRNA序列中情况类似，但未注释到的sRNA比例更高，注释到的rRNA数量减幅最大，最明显。在4n数据库中，注释到的sRNA序列数量大幅度下降，下降比例高于2n数据库（图4）。

2.5 已知miRNA与新miRNA鉴定及差异分析

通过blast或bowtie将sRNA和miRBase数据库比对，鉴定出已知miRNA。另外，对未注释上任何RNA且比对上基因组的外显子反义链、内含子、基因间区的sRNAs，通过软件Mireap鉴定新的miRNA。在2n和4n两个数据库中分别鉴定出已知miRNA为110个和172个，新miRNA为246和829个。对2n和4n数据库中的已知miRNA和新miRNA差异分析获得205个差异表达的已知miRNA和815个新miRNA。表达量分析显示差异表达的已知miRNA表达量在2n数据库中高于4n，而新miRNA表达量在2个库中相差不大，说明miRNA的调控具有物种倍性特异性。

2.6 差异表达的已知miRNA及新miRNA的GO功能分析

根据miRNA与其靶基因间的对应关系，对miRNA的靶基因的集合分別进行GO （Gene Ontology）富集分析。结果显示无论是已知miRNA还是新miRNA，GO功能都分为生物过程、细胞组分和分子功能3类，在生物过程中，miRNA靶基因主要集中在“细胞过程”、“代谢过程”和“信号过程”上；在细胞组分中，miRNA靶基因主要集中在“细胞”、“细胞器”和“代谢”上；在分子功能中则主要集中在“结合”和“催化活性”上（图6～7）。说明miRNA介导调控的生物过程较广泛。

2.7 差异表达已知miRNA及新miRNA的KEGG代谢通路分析

在生物体内，不同基因相互协调行使其生物学功能，通过Pathway显著性富集分析，能确定候选靶基因参与的最主要生化代谢途径和信号转导途径，助于更进一步了解基因的生物学功能。KEGG富集结果显示，已知miRNA主要富集在初级代谢和次级代谢途径上，新miRNA则主要富集在植物激素信号转导途径、RNA转运途径和ABC转座子途径上。这说明已知和新的miRNA在代谢途径上的差异较大，他们分别在不同的代谢途径上发挥着重要的调节作用。

2.8 抗逆差异表达miRNA分析

经过仔细查阅文献核对，在西瓜二倍体和四倍体sRNA数据库中筛选获得25个和植物抗逆胁迫相关的已知miRNA，62个和植物抗逆胁迫相关新miRNA。对于已知的miRNA，有9个在2n中上调表达，16个下调表达；而对于新miRNA则只有1个在2n中下调表达，其余61个都在2n中上调表达（表3）。根据miRNA负调控靶基因的原理推测，大部分在西瓜中上调表达的miRNA会抑制靶基因在二倍体中的表达。相反，这些miRNA会促进靶基因在西瓜四倍体中的表达。这从转录调控水平上证明了西瓜四倍体有较强的抗逆性。

3 讨论

对西瓜二倍体和同源四倍体叶片的高通量测序分别得到了11 984 042和11 600 832条clean reads，将这些clean reads与西瓜参考基因组进行对比分类注释，结果显示在西瓜二倍体（2n）和同源四倍体（4n）两个数据库中，种类和数量最多的為rRNA与miRNA，这与冬小麦基因组测序结果相同[23]，说明sRNA在植物进化中具有一定的保守性。2n和4n两个数据库中分别鉴定出已知miRNA为110个和172个，新miRNA为246和829个。明显可以看出新miRNA的数量高于已知miRNA的数量。已有的研究显示已知miRNA在进化上相对比较保守，一般只负责调控植物的生长发育等基础代谢，在严酷的进化条件下保持着序列的保守性，而新miRNA多是随着环境变化，特异性进化出的序列，参与植物的抗逆等特异性调控[24]。因此，本研究所鉴定出来的新miRNA对于西瓜属的特异性调节具有重要的参考意义。不同植物间sRNA的分布差异较大，本研究的西瓜2n中以21 nt的片段最多，在杨树、大豆和番茄中，21 nt的sRNA分布也为最多[25-27]；拟南芥、小麦和棉花中，24 nt的sRNA分布最多[28-30]，本研究的西瓜4n中以24 nt的sRNA分布为最多，说明在不同的植物中，sRNA特异性分布表达，即便是同一个品种不同倍性的植物sRNA也表现出了特异性分布表达的特点，表明植物的加倍对sRNA有很大的影响。

miRNA往往会通过表达量差异进而影响基因功能和代谢途径上的差异，这对植物的多倍体进化起着重要的调控作用。对2个库中差异表达的miRNA分析显示，已知miRNA表达量在数据库2n中高于4n，根据miRNA的反向调控作用原理分析，推测在西瓜二倍体中有更多的差异表达基因被抑制表达或是沉默了。miRNA的表达差异会对细胞内KEGG通路和GO生物过程产生重要影响[31]。本实验结果显示，差异表达的miRNA靶基因在生物过程中主要集中在“细胞过程”、“代谢过程”和“信号过程”上。差异表达的miRNA靶基因的KEGG富集结果显示已知miRNA主要富集在初级代谢和次级代谢途径上，新miRNA则主要富集在植物激素信号转导途径、RNA转运途径和ABC转座子途径上，这进一步说明了已知miRNA在进化上的保守性，新miRNA在植物进化上的特异性。

植物多倍体的抗逆性往往优于二倍体，这离不开抗逆差异表达miRNA的逆调控作用，因为miRNA控制的基因往往参与重要的抗逆和发育过程[32]。本实验结果分析显示，有80%的miRNA在西瓜二倍体中上调表达，根据miRNA对靶基因的调节特点可知，这些miRNA调节的抗逆靶基因在二倍体中的表达量会低于四倍体，甚至有可能被完全沉默掉，这会进一步导致西瓜四倍体中高度表达的抗逆基因会比二倍体中多。这在其他物种的研究中已经被证实，如二倍体泡桐中miR156b、miR28a/b等miRNAs表达量较高，相反其调控的靶基因表达量在四倍体泡桐中较高，这就导致了四倍体泡桐的抗盐性比二倍体泡桐要强[33]。因此，我们的结果为四倍体优于二倍体的抗逆性提供了转录调控水平上的证据。

参考文献

[1] 谭素英，黄秀强，刘文革. 三倍体无籽西瓜的优越性及系列无籽西瓜新品种[J]. 中国西瓜甜瓜，1994（4）：22-23.

[2] 刘文革，何 楠，赵胜杰，等. 我国三倍体无籽西瓜的周年生产与栽培[J]. 长江蔬菜，2014（14）：1-6.

[3] 郭启高，宋 明，杨天秀，等. 西瓜子叶组织培养中四倍体的产生及鉴定[J]. 西南大学学报（自然科学版），2000，22（4）：298-300.

[4] 刘文革. 不同染色体倍性西瓜（Citrullus lanatus）的遗传变异和抗逆机理研究[D]. 西安：西北农林科技大学，2003.

[5] 刘文革，阎志红，赵胜杰，等. 不同染色体倍性西瓜对枯萎病的抗性研究[J]. 长江蔬菜，2009，（18）：19-20.

[6] 朱红菊. NaCl胁迫对不同倍性西瓜幼苗影响及DNA甲基化分析[D]. 北京：中国农业科学院，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，Feng 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] 南文智，吴 嫚，于霁雯，等. 利用高通量测序技术鉴定棉纤维发育相关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，Fahlgren 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.