长叶红砂耐盐特有基因序列Rt-st11787在植物抗盐基因工程中的应用

摘要

本发明涉及一种特有盐生植物长叶红砂(Reaumuria trigyna Maxim.)的耐盐性状基因Rt-st11787。本发明选用内蒙古东阿拉善-西鄂尔多斯地区特有耐盐植物长叶红砂为试材，运用Illumina/Solexa转录组深度测序获得长叶红砂耐盐相关基因，对候选长叶红砂耐盐基因进行生物信息学分析，最终筛选出NaCl处理前后在转录组水平表达量变化明显的一种耐盐性状基因Rt-st11787。这给利用基因工程手段改造植物从而提高其抗盐能力提供了方向和靶点，对于植物在高NaCl胁迫环境下的育种试验具有重要意义。

说明书

技术领域

本发明属于植物生物技术领域，具体地说，涉及一种特有盐生植物长叶红砂的耐盐性状基因。本发明选用内蒙古东阿拉善-西鄂尔多斯地区特有耐盐植物长叶红砂为试材，运用Illumina/Solexa转录组深度测序获得长叶红砂耐盐基因，对候选长叶红砂耐盐基因进行生物信息学分析，了解基因的功能及其履行的生物学功能。

背景技术

长叶红砂(Reaumuria trigyna Maxim.)，又名黄花红砂、黄花琵琶柴，隶属于柽柳科(Tamaricaceae)琵琶柴属(Reaumuria Linn.)，为古地中海残遗珍稀植物，是亚洲中部地区东阿拉善-西鄂尔多斯特有种，为该地区重要牧草之一，其生存环境极端恶劣，为典型的高原干旱半干旱季风气候区，该地区终年干旱、低温、高盐，年均降水量140.9～302.2mm，年均温度6.0～9.2℃，土壤含盐量可达0.4％。长叶红砂是典型的泌盐盐生植物，叶表面和幼茎表面分布着大量的8细胞盐腺，这一结构可有效的将植物体内过量的盐分排出体外，研究表明该植物可在500mM NaCl胁迫下正长生长。随着人类活动范围的不断扩大如基础设施建设用地、过度放牧等对生长在该地区的植物提出了巨大的生存挑战。现在该物种已被列为内蒙古自治区珍稀濒危植物。对长叶红砂研究的意义不仅仅在于对濒危物种的保护上，更重要的是对于这一生存在极端环境下的植物的基因水平上的开发和利用。目前对长叶红砂的研究主要集中在形态学、耐盐生理机理、激素及水分调节及生物多样性上，而关于该物种盐诱导及特有耐盐基因的克隆及利用方面仍均无报道。

诞生于20世纪70年代的Sanger法是最早广泛应用的DNA测序技术，也是完成人类基因组计划的基础。但随着科学技术的不断发展，2005年以来，以Roche公司的454技术、Illumina公司的Solexa技术和ABI公司的SOLiD技术为标志的新一代测序技术相继诞生。新一代测序技术又称作深度测序技术，主要特点是测序通量高、测序时间和成本显著下降。把这种高通量测序技术应用到由mRNA逆转录生成的cDNA上，从而获得来自不同基因的mRNA片段在特定样本中的含量，这就是mRNA测序或mRNA-seq。同样原理，各种类型的转录本都可以用深度测序技术进行高通量定量检测，统称作RNA-seq或RNA测序。

目前，在已经推出的几种新一代测序平台中，Illumina/Solexa测序平台上的RNA-seq应用最广。该技术已被应用到科学研究的各个领域，如2009年发表在Science上的有关家蚕驯 化相关基因的研究中得到354个驯化候选基因，其中包括产丝重要相关基因Sgf-1；水稻研究领中发现在家养水稻及野生型水稻间有517个基因存在SNPs，目前452个源于野生稻的基因已被转入栽培稻中进行表达，其中24个基因在苗期表现出明显的表型差异，16个基因已经申请专利并被受理；此外，还运用该技术从基因组水平对青藏高原上藏人对高海拔的适应性给出了遗传学的解析，结果发现，西藏自治区2个村庄的50个藏民的外显子组，与40个北京汉人基因组数据进行比对，鉴定出一些高原适应候选相关基因。这些基因上发生的优异变异可被认为是西藏地区人类生存和繁衍后代的长期选择的重要结果。

据文献报道，关于RNA-seq技术的缺点尚未出现，但就该技术存在的难点和所面对的挑战进行了相关的描述。其中很重要的一点是高通量测序技术数据处理中的生物信息学挑战，也就是在测序所得数据做后续处理中怎样克服系统误差及解决偏好性的问题。以差异表达基因为基础的分析中，由于基因表达水平都是通过读段(reads)计数来估计的，表达水平较高或转录本较长的基因拥有更多的读段，更容易被多数统计方法识别为差异表达基因。这种偏好可能对后续分析带来影响。针对这种偏好性，Young等人发展了一种方法，较好的避开了这种偏好性。另外，在植物耐盐基因克隆及利用方明，目前大多集中在模式植物拟南芥、盐芥或水稻，小麦等作物上。对野生天然强抗耐盐植物的研究还比较少见，尤其是结合高通量测序技术系统性发掘和筛选通过环境长期筛选出的盐胁迫相关基因。本发明获得的耐盐性状基因正是利用新一代测序技术对内蒙古特有耐盐植物长叶红砂转录组进行深度测序所得。

发明内容

本发明的目的是利用二代测序技术illumina/solexa系统，对内蒙古东阿拉善-鄂尔多斯特有盐生植物长叶红砂(室内培养幼苗)转录组进行De novo深度测序。对测序结果进行生物信息学分析后，得到该物种转录组中全部与非生物胁迫相关的基因序列(全长或部分)。最终筛选出NaCl处理前后在转录组水平表达量变化最明显的耐盐相关基因序列。

本发明的实施方案是选用内蒙古东阿拉善-西鄂尔多斯地区特有耐盐植物长叶红砂为试材，运用Illumina/Solexa转录组深度测序获得长叶耐盐基因，对候选长叶红砂耐盐基因进行生物信息学分析，了解基因的功能及其履行的生物学功能。该方法除通量大、速度快、成本低的优点外，还可应用于无全基因组背景的物种的研究上，并且，该方法产生的海量数据，可以为后续深入研究某一物种特定代谢机理提供可靠信息。

本发明的一个目的在于提供新的长叶红砂耐盐性状的基因，定名为Rt-st11787，其序列为SEQ ID No：1或SEQ IDNo：2所示的序列。

本发明的另一个目的在于提供上述基因编码的多肽。

本发明的再一个目的在于提供含有上述基因序列的表达载体。

根据本发明的一方面，一种在NaCl胁迫下表达量明显上升的基因Rt-st11787，来源于长叶红砂(Reaumuria trigyna Maxim.)，名称为Rt-st11787，其核苷酸序列如序列表中SEQ ID NO：1所示，其编码区如序列表中SEQ ID NO：2所示，其编码的325个氨基酸序列，如序列表中SEQ ID NO：3所示。该基因是参与渗透胁迫通路MAPK信号级联系统前体MAPKK家族成员。

具体地说，本发明提供了含有下述序列之一的分离的多核苷酸：

(1)序列表中的SEQ ID NO：1和SEQ ID NO：2；

(2)与序列表中SEQ ID NO：1或SEQ ID NO：2限定的DNA序列具有90％以上同源性，且编码相同功能蛋白质的DNA序列；

上述涉及的多核苷酸还包括取代、缺失、和插入变体以及等位变体、剪接变体、片段、衍生物等，其中可以通过取代、缺失、插入或衍生一个或多个核苷酸。优选的这些多核苷酸是那些具有同样耐盐性状的功能。

本领域技术人员可以理解的是，上述的分离的多核苷酸也包括那些与SEQ ID NO：1或SEQ ID NO：2所示序列具有较高同源性的序列，例如同源性大于90％、甚至95％、甚至98％的序列；还包括那些在严谨条件下可与SEQ ID NO.1或SEQ ID NO：2所示序列杂交的序列；或者可与SEQ ID NO.1或SEQ ID NO：2所示序列互补的序列。

根据本发明的另一方面，本发明亦提供包括一种或多种上述多核苷酸的重组载体。在优选的实施方案中，这种重组载体包括上述的多核苷酸，它编码含SEQ ID NO：3的多肽，其含SEQ ID NO：1或SEQ ID NO：2的序列所示的多核苷酸。

本发明同时提供了将本发明的新基因Rt-st11787应用于植物抗盐基因工程中的技术方案。

本发明从长叶红砂中成功地分离获得了在NaCl胁迫下表达量明显上升的基因Rt-st11787，这给利用基因工程手段改造植物从而提高其抗盐能力提供了方向和靶点，对于植物在高NaCl胁迫环境下的育种试验具有重要意义。具体来说，本发明的具体实施方案之一是将Rt-st11787基因应用于植物基因工程以提高植物在高NaCl胁迫环境下生存能力的途径。

以上概括的描述了本发明，可通过参照本文提供的某些具体实施例进一步理解本发明，这些实施例仅是为了说明而不是限制本发明。

附图说明

图1：对建好的测序文库用Illumina HiSeq^TM 2000进行测序的流程图；

图2：unigene生物信息学分析的流程图，clean reads代表测序所得原始reads(读段)去 掉只含adaptor(接头)序列的部分，Unigenes代表经组装所得序列；

图3：将测序组装所得的65340条unigene通过COG注释到25个分子家族，纵轴代表Unigenes的数量，横轴各字母代表25种基因功能分类：其中A代表RNA加工和修改(RNA processing and modification)；B代表染色质结构和动力学(Chromatin structure and dynamics)；C代表能源生产和转换(Energy production and conversion)；D代表细胞周期调控，细胞分裂，染色体分离(Cell cycle control，cell division，chromosome partitioning)；E代表氨基酸转运与代谢(Amino acid transport and metabolism)；F代表核苷酸转运与代谢(Nucleotide transport and metabolism)；G代表碳水化合物运输和代谢(Carbohydrate transport and metabolism)；H代表运输和代谢的辅酶(Coenzyme transport and metabolism)；I代表脂质转运和代谢(Lipid transport and metabolism)；J代表翻译、核糖体结构和合成(Translation，ribosomal structure and biogenesis)；K代表转录(Transcription)；L代表复制，重组和修复(Replication，recombination and repair)；M代表细胞壁/膜/胞膜合成(Cell wall/membrane/envelope biogenesis)；N代表细胞运动(Cell motility)；O代表翻译后修饰(Post-translational modification)；P代表无机离子转运与代谢(Inorganic ion transport and metabolism)；Q代表生物合成，运输和代谢的次生代谢产物(Secondary metabolites biosynthesis，transport and catabolism)；R代表一般功能预测(General function prediction only)；S代表功能未知(Function unknown)；T代表信号转导机制(Signal transduction mechanisms)；U代表细胞内，分泌和膜泡运输胞(Intracellular trafficking，secretion，and vesicular transport)；V代表防御机制(Defense mechanisms)；W代表胞外结构(Extracellular structures)；Y代表核结构(Nuclear structure)；Z代表细胞骨架(Cytoskeleton)；

图4：在5032条差异表达基因，其中上调(T43相对于C21表达量增加)基因2370条，下调(C21相对于T43表达量增加)基因2662条，图中由两条曲线划分出三个区域，偏上区域代表上调差异表达基因，偏下区域代表下调差异表达基因，中间区域为未检测到的差异表达基因(Not detected expression genes)。

具体实施方式

实施例1.长叶红砂幼苗培养

长叶红砂种子于2008年9月采自内蒙古自治区东阿拉善-西鄂尔多斯地区，挑选饱满种子用10％次氯酸钠灭菌15min，灭菌ddH₂O冲洗三遍，播种于装有40ml MS培养基的150ml三角瓶中。暗培养72h后在25℃、湿度70％、16h:8h光照/黑暗的条件下培养。幼苗生长至10cm(大约15d)左右，转移至25cm高，直径5cm的大试管中培养，每两天换一次灭菌1/2Hoagland营养液，每次约50ml。挑选三株长势相近的长叶红砂幼苗为实验试材，NaCl处理 前，适当剪取三株幼苗茎叶；对幼苗进行100mM NaCl处理后，于0.5h、1h、2h、4h及8h收取适量幼苗，NaCl浓度也随即增加到200mM，8h后收取材料并重复该步直至收获到400mM NaCl处理材料为止，所有材料在收获后应迅速保存在液氮中备用。

实施例2.RNA-seq cDNA文库的制备

a.total RNA的提取

采用plant plus RNA regent(DP437，Tiangen，Beijing)根据用法说明进行总RNA的提取，对所得总RNA进行45分钟DNase I(TaKaRa)消化，取1.5μl进行1％琼脂糖凝胶电泳，Nanovue plus测定260/280OD值均在1.9～2.1之间、260/230均大于2.0。将对照及盐胁迫处理样品分别等量混合，得到两样本RNA池：C21、T43。最后对两样品做Agilent 2100检测。

b.测序文库的制备

总RNA用带有Oligo(dT)的磁珠富集mRNA后将mRNA片段化，以mRNA为模板，用六碱基随机引物(random hexamers)合成cDNA第一链，然后加入缓冲液、dNTPs、RNaseH和DNA polymerase I合成cDNA第二链。反转录产物经QiaQuick PCR试剂盒纯化后做末端修复、加poly(A)并连接测序接头，琼脂糖凝胶电泳进行片段大小选择后进行PCR扩增，建好的测序文库用Illumina HiSeq^TM 2000进行测序(流程图可见附图1)。

实施例3.测序结果组装及生物信息学分析

a.测序结果组装

使用短reads组装软件SOAPdenovo(组装软件)做转录组从头组装。SOAPdenovo首先将具有一定长度overlap(重叠)的reads(读段)连成更长的片段，这些通过reads overlap关系得到的不含N的组装片段我们称之称为Contig。

将reads比对回Contig，通过paired-end reads能确定来自同一转录本的不同Contig以及这些Contig之间的距离，SOAPdenovo将这些Contig连在一起，中间未知序列用N表示，这样就得到Scaffold。进一步利用paired-end reads对Scaffold做补洞处理，最后得到含N最少，两端不能再延长的序列，我们称之为Unigene。如果同一物种做了多个样品测序，则不同样品组装得到的Unigene可通过序列聚类软件TGICL做进一步序列拼接和去冗余处理，得到尽可能长的非冗余Unigene。

b.unigene生物信息学分析(分析流程见附图2)

通过blastx将Unigene序列比对到蛋白数据库nr、Swiss-Prot、KEGG和COG(evalue＜0.00001)，得到跟给定Unigene具有最高序列相似性的蛋白，从而得到该Unigene 的蛋白功能注释信息。

GO注释将组装所得的65340条unigene注释到三个ontology，分别是分子功能(molecular function)、所处的细胞位置(cellular component)、参与的生物过程(biological process)。对于GO注释而言，Gene Ontology(简称GO)是一个国际标准化的基因功能分类体系，提供了一套动态更新的标准词汇表(controlled vocabulary)来全面描述生物体中基因和基因产物的属性。GO总共有三个ontology(本体)，分别描述基因的分子功能(molecular function)、所处的细胞位置(cellular component)、参与的生物过程(biological process)。GO的基本单位是term(词条、节点)，每个term都对应一个属性。GO功能分析一方面给出差异表达基因的GO功能分类注释；另一方面给出差异表达基因的GO功能显著性富集分析。

COG注释将所得unigene注释到25个分子家族，主要包括细胞结构、生物化学代谢、分子过程、信号转导、基因表达以及刺激响应等，如附图3所示。对于COG注释而言，COG是Cluster of Orthologous Groups of proteins(蛋白相邻类的聚簇)的缩写。构成每个COG的蛋白都是被假定为来自于一个祖先蛋白，并且因此或者是orthologs或者是paralogs。Orthologs是指来自于不同物种的由垂直家系(物种形成)进化而来的蛋白，并且典型的保留与原始蛋白有相同的功能。Paralogs是那些在一定物种中的来源于基因复制的蛋白，可能会进化出新的与原来有关的功能。

对于KEGG(京都基因与基因组百科全书)而言，它是基因组破译方面的数据库。在后基因时代一个重大挑战是如何使细胞和有机体在计算机上完整的表达和演绎，让计算机利用基因信息对更高层次和更复杂细胞活动和生物体行为作出计算推测。为达到此目的，人们建立了一个在相关知识基础上的网络推测计算工具。在给出染色体中一套完整的基因的情况下，它可以对蛋白质交互(互动)网络在各种细胞活动起的作用作出预测。KEGG的Pathway数据库整合当前在分子互动网络(比如通道，联合体)的知识。

实施例4.差异表达基因的获得和其功能注释

经测序共得非冗余序列65340条，比较C21及T43两样本测序结果得到差异表达基因信息，差异表达基因筛选参照Audic S.等人发表在Genome Research上的基于测序的差异基因检测方法，公式如下：

$p\left(x\right)=\frac{e^{-\lambda }{\lambda }^x}{x!}$(λ为基因A的真实转录数)

假设观测到基因A对应的reads(测序所产生的短片段)数为x，已知在一个大文库中，每个基因的表达量只占所有基因表达量的一小部分，在这种情况下，p(x)的分布服从泊松分布。已知，样本一能比对到所有Unigene的总reads数为N1，样本二能比对到所有Unigene 的总reads数为N2，基因A在样本一中对应的reads数为x，在样本二中对应的reads数为y，则基因A在两样本中表达量相等的概率可由以下公式计算：

$2\underset{i=0}{\overset{i-y}{\Sigma }}p\left(i\right|x)$

或$2\times (1-\underset{i=0}{\overset{i-y}{\Sigma }}p\left(i\right|x))$(如果$\underset{i=0}{\overset{i-y}{\Sigma }}p\left(i\right|x)>0.5)$

$p\left(y\right|x)={\left(\frac{N_2}{N_1}\right)}^y\frac{(x+y)!}{x{!y!(1+\frac{N_2}{N_1})}^{(x-y-1)}}$

然后，对差异检验的p value作多重假设检验校正，通过控制FDR(False Discovery Rate)来决定p value的域值。假设挑选了R个差异表达基因，其中S个是真正有差异表达的基因，另外V个是其实没有差异表达的基因，为假阳性结果。希望错误比例Q＝V/R平均而言不能超过某个可以容忍的值(比如1％)，则在统计时预先设定FDR不能超过0.01。在得到差异检验的FDR值同时，同时也根据基因的表达量(RPKM值：说明如下)计算该基因在不同样本间的差异表达倍数。FDR值越小，差异倍数越大，则表明表达差异越显著。在我们的分析中，差异表达基因定义为FDR0.001且倍数差异在2倍以上的基因。

RPKM法：其计算公式为：

$\mathrm{RPKM}=\frac{{10}^{6}C}{\mathrm{NL}/{10}^3}$

设RPKM(A)为Unigene A的表达量，则C为唯一比对到Unigene A的reads数，N为唯一比对到所有Unigene的总reads数，L为Unigene A的碱基数。RPKM法能消除基因长度和测序量差异对计算基因表达的影响，计算得到的基因表达量可直接用于比较不同样品间的基因表达差异。

经计算，在限定条件内产生5032条差异表达基因，其中上调(T43相对于C21表达量增加)基因2370条，下调(C21相对于T43表达量增加)基因2662条，如附图4所示，图中偏上区域代表上调差异表达基因，偏下区域代表下调差异表达基因，中间区域为未检测到的差异表达基因(Not detected expression genes)。

对5032条差异表达基因做GO功能注释、KEGG Pathway分析(结果与上述类似)。经GO显著性富集分析及差异表达模式聚类分析，以及KEGG pathway显著性富集分析可得知差异表达基因富集的代谢通路及执行的主要生物学功能，从该结果中可得到直接被注释到与非生物胁迫响应(response to stress)相关的全部基因，统计共得到200余条，其中上调基因数目略大于下调基因。

通过上述生物学注释，我们筛选出三条长叶红砂在NaCl胁迫下表达量明显上升的gene，分别是参与渗透胁迫通路MAPK信号级联系统前体MAPKK家族成员，与非生物胁迫密切相 关的WRKY转录因子家族成员以及氧化呼吸酶。其中一个gene-Rt-st11787，为参与渗透胁迫通路MAPK信号级联系统前体MAPKK家族成员，详细信息如表1所示，此外还利用ESTscan软件预测了该序列的编码区(coding region)序列及可能翻译的氨基酸序列。

表1：Rt-st11787基因相关信息。

注：Depth：测序深度(等于reads碱基数/unigene序列长度；Coverage：覆盖度(等于unigene中有reads覆盖的碱基数/unigene序列长度)

以上实施例进一步说明了本发明的内容，但不应理解为对本发明的限制。在不背离本发明精神和实质的情况下，对本发明方法、步骤或条件所作的修改或替换，均属于本发明的范围。