海南山栏稻抗旱基因转录组分析

翟楠鑫 迟会 夏玥琳 刘彩月 裴新梧 袁潜华

（1. 海南大学热带作物学院，海口 570228；2. 中国农业科学院生物技术研究所，北京 100081）

山栏亦作“砍山”或“砍山栏”，山栏稻是海南特有的旱稻，主要分布于保亭、白沙和乐东等山坡地，盛产于五指山地区。不同于水稻的育秧插秧，山栏稻是点播直播，主要种植在山坡地，产量低，加之种植规模小，品种老化，导致山栏稻种资源流失严重［1］。山栏稻具有抗旱、耐热、抗病、抗虫等优良特性，是珍贵的旱稻种质资源［2］。

山栏稻的农艺性状一般表现为植株偏高，穗形披散，结实率较高，粒型多为大长粒，叶片宽大而披垂，冠层密封分蘖力弱，穗大而弯垂等特点［2］。刘维俊等［4］研究了干旱胁迫条件下山栏稻与栽培稻品种苗期表型性状及生理差异，结果表明山栏稻苗期具有较强的抗旱性。唐力琼等［5］采用PEG人工模拟干旱的方法，研究20个山栏稻资源在芽期和苗期的抗旱性，发现山栏稻品种总体抗旱性较强。徐建欣等［3］对22个海南山栏稻品种进行了全生育期抗旱性鉴定，发现山栏稻品种不同生育期间所表现出的抗旱能力具有极显著的正相关关系，这表明山栏稻可以在早期（芽期或苗期）快速筛选抗旱品种。在山栏稻品种选育方面，研究人员也成功培育出了“热大99 W”［6］和“山栏陆1号”［7］等抗旱丰产品系。刘跃林等［8］通过基因芯片技术，研究海南山栏稻在干旱逆境条件下的基因表达谱，将其中的OsMS17在水稻中超表达，转基因植株的苗期和成株期的抗旱性明显增强。

早在2013年，本研究团队就对山栏稻展开了研究，刘欣欣等［9］克隆了海南山栏稻的HKT2基因片段，研究了海南山栏稻运输Na+和K+的分子机制；刘维俊等［4］研究了干旱胁迫下2个山栏稻与栽培水稻品种苗期表型性状及生理差异；吴丹等［10］通过田间试验设计覆盖秸秆旱作、水作和传统旱作3种栽培方式，发现山栏稻在水田环境下具有减少干旱胁迫、优化农艺性状和提高产量等诸多优点并建立了山栏稻水田灌溉种植的高产栽培方法［11］，为山栏稻从传统山区旱地种植引向普通旱作和水田种植提供科学的依据。我们在对山栏稻前期研究的基础上，通过RNA-seq（RNA Sequencing）Illumina 二代高通量测序技术，对山栏稻抗旱相关基因转录组进行表达分析及鉴定，以进一步明确海南山栏稻在干旱逆境条件下的基因表达谱，为解析山栏稻的抗旱机理奠定基础。

1 材料与方法

1.1 材料

山栏稻品种白沙糯，由海南大学热带作物学院生物多样性与水稻种质创新实验室收集保存，该品种遗传稳定，根系发达，耐旱性强。

1.2 方法

1.2.1 植物材料培育 山栏稻品种白沙糯在光照培养箱中水培，培养条件如下：相对湿度50%-70%；光照16 h，30℃/黑暗8 h，26℃培养。水培营养液用Yoshida水稻营养液配制，2 d更换一次。

1.2.2 建库和转录组测序 山栏稻白沙糯培养至三叶一心期后，对照组继续营养液培养，处理组加20%的PEG-6000（聚乙二醇6000）培养液模拟干旱胁迫处理。取0 h、24 h、48 h、72 h的水稻样本分别作为不同胁迫时段材料并提取总RNA。RNA样品检测合格后，用3个生物学重复的对照组和处理组（胁迫72 h）构建6个cDNA文库，库检合格后用Illumina HiSeq 2500进行测序。

1.2.3 数据过滤、重新序列组装以及unigenes的功能注释 将测序得到的原始序列Raw reads进行过滤得到Clean reads。取大于200 bp且序列最长的转录本作为unigene，将6个库混合拼接成一个转录本。序列拼接后与NR（Non-redundant）、PFAM（Protein Families）、Swiss-Prot（http：//www.expasy.ch/sprot/）和KEGG（Kyoto Encyclopedia of Genes and Genomes）等数据库进行比对。GO（Gene Ontology）功能注释基于NR和PFAM两部分的蛋白注释结果，软件为Blast2GO v2.5和WEGO。并根据它们参与的KEGG代谢通路进行分析，KEGG相关注释通过KEGG Automatic Annotation Server完成。

1.2.4 CDS预测基因和表达水平分析 CDS（Coding sequence）预测是按照NR、Swissprot、KEGG GENES的优先级顺序，将unigene比对到蛋白数据库。采用RSEM（RNA-Seq by Expectation-Maximization）软件，对bowtie的比对结果进行统计，得到每个样品比对到每个基因上的readcount数目，并对其进行FPKM（Fragments Per Kilobase of exon model per Million mapped reads）转换，进而分析基因的表达水平。

1.2.5 差异表达unigenes的验证和表达分析 在基因表达水平计算的基础上，我们采用DESeq进行分析。为了筛选到更多的与干旱胁迫相关的候选基因，对干旱胁迫处理前后的阈值设置为P-value＜0.1，|log2FoldChange|＞1。

为保证转录组测序的差异表达基因表达趋势可靠性，我们从测序获得基因中随机选取差异表达基因进行qRT-PCR验证其是否和转录组库趋势一致；另外选取了已报道干旱胁迫相关基因进行验证，鉴定它们与干旱胁迫的相关性。qRT-PCR采用SYBR Green结合相对定量法在ABI 7500 Real Time System上进行。选用OsActin作为内参，所有样本作3个重复。序列见表1。

表1 qRT-PCR分析特异性引物

2 结果

2.1 转录组拼接质量分析及参考序列比对

将Illumina测序结果去除衔接子和质量差的读数后，把Clean reads映射到参考基因组数据中，对照组的3个重复基因序列个数分别为：38 504 957、35 184 648和26 794 172个；胁迫组的3个重复基因序列个数分别为：32 620 538、29 253 860和30 195 571个。这6个库的测序错误率均小于0.05%，Q20均大于95.00%，GC含量也在50%左右，表明这6个库的测序结果可靠，能够用于后续分析。将Clean reads与参考序列进行比对（表2），比对到参考序列上的Reads数量占总数的92%左右，其中比对到多个位置的Reads占总数的2.6%以下。

表2 Reads与参考序列比对情况一览表

2.2 SNP与InDel

SNP（Single Nucleotide Polymorphisms），指 基因组水平上由单个核苷酸变异引起的DNA序列多态性。SNP主要为转换和颠换两种类型，一般二者的比例为2∶1。因为其广泛存在于基因组上，现在成为了遗传研究的主要标记。InDel（Insertion and Deletion）指基因组上小片段的插入和缺失，编码区的InDel往往会引起基因功能的重大改变。我们首先通过samtools和picard-tools等工具对比对结果进行染色体坐标排序，然后通过变异检测软件GATK分别进行SNP Calling和InDel Calling，最后进行过滤，得到可靠的SNP和InDel检测结果（图1）。

2.3 基因表达对比

使用DESeq对得到的数据进行基因差异表达分析，结果（表3）表明：山栏稻白沙糯干旱胁迫前后有2 791个差异表达的基因，其中1 669个表达上调，1 122个表达下调。与目前已知的水稻基因组数据对比，其中已知基因共2 650个，未知基因共141个。

图1 SNP和InDels的统计数据

表3 干旱胁迫前后差异表达基因列表

2.4 差异基因GO富集分析和KEGG富集分析

为研究山栏稻基因调控与干旱胁迫的关系，在NR功能注释基础上，我们对胁迫前后差异表达基因进行了GO分析。干旱胁迫前后差异表达基因共2 791个，有GO注释的差异基因有1 599个，被分配到了1 524个GO terms功能组，所占比例最大的为生物过程和分子功能。我们统计出这1 524个组中有25个GO terms与胁迫响应相关，其中对应激的反应所占比例最大（基因个数与所占比例分别为：164，10.26%，下同），其次从大到小分别为：对压力的响应（84，5.25%）、细胞对应激的响应（70，4.38%）、对化学应激响应（57，3.56%）、对非生物应激的反应（30，1.88%）、氧化胁迫响应（23，1.44%）和对水的反应（20，1.25%）。

如图2-A所示，纵坐标为富集的GO term，横坐标为该term中差异表达基因个数。带“*”为显著富集的GO term。在1 524个GO terms中达到显著富集有24个，可大致分为以下几类：蛋白质磷酸化、蛋白激酶活性、磷酸化与磷代谢过程、修饰过程、催化活性、对非生物刺激的反应、转移酶活性、离子结合、阳离子结合、几丁质分解代谢过程、水解酶活性、氨基聚糖分解代谢过程和对水的反应。

KEGG是系统分析基因功能、基因组信息的数据库，它有助于研究者把基因及表达信息作为一个整体网络进行研究。在生物体内，不同基因相互协调行使其生物学功能，通过Pathway显著性富集分析，找出与整个转录组背景相比显著性富集的Pathway。如（图2-B）纵轴表示Pathway名称，横轴表示Rich factor，点的大小表示此Pathway中差异表达基因个数多少，而点的颜色对应于不同的q-value范围，q-value的值越接近于零，说明富集程度越显著。可以看出干旱胁迫前后差异表达基因KEGG显著富集在植物激素信号转导，次生代谢物的生物合成，角质、木栓质和蜡生物合成中。

2.5 差异基因转录因子对比

转录因子（Transcription factor，简称TF）可以通过结合靶基因启动子的顺式作用元件，从而调控基因的表达。转录因子可以控制一系列下游胁迫反应，从而启动信号传导中的级联反应，使细胞产生相应的抗逆性，分析这些转录因子可以为水稻干旱胁迫下复杂的调控网络提供重要的信息。对差异基因进行筛选，发现可能与抗旱有关的转录因子，如表4所示为差异表达转录因子基因信息（部分），将库中与抗旱有关的转录因子进行分类统计得到如下结果：Znf（Zinc finger）家族转录因子在干旱胁迫后上调的数量最多，达41个，占相关上调转录因子总数的35.65%（表5）。

2.6 差异表达基因验证

为了验证测序结果的可靠性，我们从干旱胁迫相关的已知基因中随机选取6个基因并利用qRTPCR进行验证，这6个基因的相关注释信息以及在库中的表达量如表6所示；另外选取已报道的与干旱胁迫相关的基因RAB21和OsLEA3-2进行试验，验证干旱胁迫处理是否适当。图3是8个基因的qRT-PCR结果：RAB21和OsLEA3-2在胁迫后表达量显著升高，表明模拟干旱胁迫处理适当。基因Os02g0693700、Os07g0658600和Os11g0115100为已注释的上调基因；基因Os09g0351700和Os07g025-1900为已注释下调基因；基因Os01g0910800为未注释表达上调基因，且上调倍数明显。这6个基因与测序的表达趋势基本一致，说明本次研究的RNAseq高通量测序结果可靠。

图2 差异基因GO富集柱状图（A）及差异基因KEGG富集散点图（B）

2.7 干旱胁迫相关差异基因不同胁迫时段的表达检测

验证测序可靠后，运用qRT-PCR技术，我们从干旱胁迫相关的候选基因中选取了4个基因来检测这些基因在山栏稻白沙糯不同胁迫时期（0 h、24 h、48 h和72 h）的表达情况，这4个 基因的相关注释信息以及在库中的表达量如表6所示。其中基因Os01g0910800为未注释上调基因；基因Os02g0693700和基因Os07g0658600为上调已注释基因；基因Os07g0251900为下调已注释基因。从图4的结果可以看出：基因Os07g0658600、Os02g0693700和Os01g0910800随着胁迫时间的延长，表达逐渐上升与干旱胁迫正相关；而基因Os07g0251900随着胁迫时间的延长，表达逐渐下降与干旱胁迫负相关。

表4 差异表达转录因子基因信息（部分）

表5 不同库中差异表达转录因子分布

3 讨论

本研究对山栏稻白沙糯进行了干旱胁迫应答的转录组分析发现，所有差异表达的基因都为单基因，有2 791个差异基因受干旱胁迫诱导，其中1 586个基因在胁迫后表达上调，1 065个基因在胁迫后表达下调。转录因子可以控制一系列下游胁迫反应，从而启动信号传导中的级联反应并使细胞产生相应的抗逆性。目前与水稻干旱胁迫相关的转录因子主要有：bZIP、NAC、DREB/CBFs、Zinc finger蛋白、HSP和MYB等转录因子［12-17］。在本研究中获得的山栏稻抗旱转录组数据中，有184个差异表达转录因子，在干旱胁迫后上调的转录因子主要有：Znf、AP2/ERF、MYB、NAC、bZIP等家族转录的因子，其中Znf家族转录因子在胁迫后上调的数量最多，占相关上调转录因子总数的35.65%。这些上调的转录因子与田新杰［18］对海南野生稻和黄立钰［19］对水稻H471干旱处理后转录组测序的差异表达转录因子的类型相似，不同的是山栏稻干旱胁迫处理后差异表达的转录因子数量比前两者的数量更多，这也从侧面说明了山栏稻抗旱性更强。

山栏稻白沙糯差异基因的GO分类分析显示，1 599个单基因在功能上用GO术语注释，分属1 524个GO terms功能组。山栏稻白沙糯显著富集在25个terms，包括蛋白质磷酸化、非生物应激反应、几丁质分解代谢、蛋白激酶活性、金属离子结合和几丁质酶活性等terms。蛋白激酶主要起着催化蛋白质磷酸化的作用，而蛋白质的可逆磷酸化是细胞信号识别与转导的重要环节，与植物干旱、高盐、低温、激素等反应的信号传递有关［20］；大多数几丁质酶也可受某些非生物因素的诱导，其诱导表达是植物对逆境的一般反应［21］；金属离子如K+、Ca2+和Mg2+与作物的耐旱性极为密切，它们是维持作物细胞渗透压的最主要离子［22］。这些研究结果表明上述通路的确与抗旱有关，与本研究结果一致。

表6 干旱胁迫前后差异表达基因信息

图4 干旱胁迫相关基因不同胁迫时段的表达检测

KEGG通路分析结果显示，山栏稻白沙糯差异基因的富集在109个通路中。在植物激素信号转导，次生代谢物的生物合成，角质、木栓质和蜡生物合成这3个途径中显著富集。植物激素是在植物体内合成，对生长发育产生显著作用的微量有机物质。前人研究表明，干旱条件下各种植物激素调节通过相互影响协同作用来抵御干旱［23］。角质层覆盖在植物最表层，是所有陆生植物最重要的保护层之一。角质层包括角质和蜡质，在植物生长发育过程中可以保护植物免受环境压力，使其地上组织适应水分胁迫的环境及其它非生物和生物逆境，我们研究的结果也很好的印证了这一点。在水分胁迫下，植物组织中次生代谢产物的浓度常常上升，包括氰甙、其他硫化物、萜类化合物、生物碱、单宁和有机酸等［24］。对比上述研究及本试验结果，我们认为山栏稻白沙糯主要通过渗透调节物质（金属离子）、角质蜡质保护和激素调节来抵御干旱胁迫。

本研究首次利用高通量测序技术对海南旱稻山栏稻进行转录组分析及表达鉴定，我们发现有1 599个单基因在功能上有注释，有1 192个差异基因是未注释基因。其中以基因Os01g0910800为代表的若干未注释基因在干旱胁迫诱导后上调，且倍数非常高，说明这些基因可能受干旱胁迫诱导表达，与水分胁迫相关，这些未注释基因的功能有待我们去研究和探索。我们获得了海南山栏稻在干旱胁迫条件下的基因表达谱，期望为旱稻的进一步深入研究及利用提供理论依据。

4 结论

山栏稻白沙糯干旱胁迫后有2 791个基因差异表达，未注释基因1 192个。其中1 669个基因表达上调，1 122个基因表达下调；对差异基因中的转录因子进行分析，发现山栏稻白沙糯中有184个差异表达的转录因子，其中Znf家族转录因子是在干旱胁迫后上调的数量最多的，共41个。本研究发现山栏稻白沙糯通过植物激素、角质阻力增大和渗透压调节来提高抗旱性。一系列未注释基因（如基因Os01g0910800）随着胁迫时间的延长，表达逐渐上升且增长倍数明显。