扁秆荆三棱不同组织混合样品的转录组分析

宁 宇, 郑彦超, 武高洁, 王义飞,4, 张曼胤

(1.中国林业科学研究院 湿地研究所，北京 100091；2.四川若尔盖高寒湿地生态系统国家定位观测研究站，四川 若尔盖 624500； 3.国家林业局华东林业调查规划设计院，浙江 杭州 310019；4.湿地生态功能与恢复北京市重点实验室，北京 100091； 5.河北衡水湖湿地生态系统国家定位观测研究站，河北 衡水 053000)

扁秆荆三棱(Bolboschoenusplaniculmis)是河漫滩、湖岸、稻田等湿地类环境中的常见杂草。该物种在自然状态下能通过地下根茎和球茎等结构不断进行营养繁殖而产生大量无性系分株，这种过程被称为克隆生长[1-2]。克隆生长赋予该物种生理整合、内部分工等高效利用资源和适应环境的能力，使其经常在自然群落中占据绝对优势地位[2-3]。已有的研究通常将其视为水稻田中的恶性杂草而研究其防除技术[4-6]。钱希[7]建立了预测扁秆荆三棱再生苗和球茎生长量的线性回归方程，指出机械和人工移除虽然暂时抑制生长，却促进了再生苗的蔓延。康学耕和唐恩全等[8]认为其无性繁殖平均10 d一代，并且当年生的种子当年不能萌发，至次年6月中旬萌发出土形成实生苗，形成“即时型”和“延后型”相结合的爆发，对稻田生产极具威胁。该植物虽然在稻田中通常被视为恶性杂草，但其在不同类型的系统中却可能具有多种应用潜力。在自然湿地中，其球茎由于淀粉含量较高，成为国家一级保护动物白鹤(Grusleucogeranus)在我国东北地区湿地停歇时的主要植物性食物[9]，并且有研究表明它可能是一种历史悠久的传统食物[10]，具有食品加工方面的应用潜力。该物种还具有较好的去污染及净化水体能力，可以作为人工湿地的构建材料[11-12]。但是，目前的工作多属于经典生物学和实验生态学范畴，在分子水平的信息比较匮乏，不利于生物学机理的深入研究；因此，亟待补充相关数据，为农林相关工作提供科学依据。

近年来，基于高通量方法的转录组测序正成为植物生理生态研究领域的有力工具。转录组测序可以获得某一物种特定组织或器官在某一状态下的几乎所有转录本序列信息，进而得到的基因功能注释、蛋白质编码区序列、基因的表达量、代谢途径等大量信息，为进一步研究提供基础数据和重要参考[13-14]。该方法不仅能够用于有参考基因组序列的物种研究，如水稻[15-16]、葡萄[17]等，也能用于无参考基因组序列的物种[18-20]，并且在探讨基因的差异性表达(Differential expression)和共表达(Co-expression)情况等方面具有较强功效。已有的研究显示，扁杆荆三棱的地下球茎网络系统在不同环境下具有很强的可塑性，根茎的长度、球茎的数量和营养成分都可能存在差异[3, 21-22]，从而造成不同的种群扩张策略。对该物种进行转录组分析，有助于获取其生长调节相关的基因信息，并为进一步深入研究不同环境条件下的表达模式提供基础。

针对现有的研究中扁秆荆三棱转录组信息缺失的情况，本研究采用高通量测序技术对扁秆荆三棱的茎、叶、根茎和球茎的混合样品进行转录组测序，结合生物信息分析对获得的Unigene开展功能注释、代谢通路分析等方面的研究，并探讨其潜在的SSR位点情况，以期为扁秆荆三棱遗传多样性分析、功能基因挖掘和利用、分子辅助育种等方面的研究提供数据和理论支持。

1 材料和方法

1.1 试验材料

1.1.1 采集地点 样品采集于野外自然种群，具体位置为河北衡水湖湿地生态系统国家定位观测研究站湖岸生态监测区域内(115°38′39.84″，37°39′10.08″)。采样时间为2018年7月份，该时间内扁杆荆三棱多处于克隆生长的起始阶段，生长迅速[22]，较为适合采集不同组织部位的样品。采样种群所在区域属暖温带大陆季风气候区，整个湖区的地形为浅碟状洼地，湖岸为自然平地，土壤以潮湿的黏土为主要类型。保护区内全年平均降水量为506.3 mm，主要集中于7-8月份，全年日照平均总时数为2 571.7 h[23]。

1.1.2 取样方法 选择长势均匀的扁秆荆三棱种群，在其范围内随机设置5个1 m×1 m的小样方，每个小样方内随机选择3个单株，分别采集其茎、叶、球茎和根茎部分，分开盛装并做好标记。总共采集15个扁秆荆三棱的单株，每个单株的取样量为5 g左右。各部分采样后立即使用去离子水冲洗去除杂质，然后封装放入足量的干冰环境中暂时保存，并在24 h内存于-80 ℃低温条件下，以避免RNA降解。

1.2 试验方法

1.2.1 RNA的提取、建库与测序 分别取茎、叶、根茎和球茎1 g的材料，采用RNeasy Plant Mini Kits(Qiagen, Inc., Valencia, CA, USA)试剂盒提取总RNA，对总RNA样品使用DNase(TaKaRa, Japan)纯化后，从每类组织的产物中取10 μg等量混合组成RNA池，然后使用Oligo(dT)磁珠分离和富集poly(A)mRNA，向富集后的mRNA中加入Fragmentation Buffer 并进行打断，再以片段化的mRNA为模板，用随机引物(Random hexamers)合成第一链cDNA，加入缓冲液、dNTPs、RNaseH 和DNA polymerase 合成第二链cDNA。经过试剂盒纯化回收、黏性末端修复后，在cDNA 的3′ 末端加上碱基“A”并连接接头，用琼脂糖凝胶电泳进行片段大小选择，然后进行PCR 扩增得到测序文库。构建好的文库用Agilent 2100 Bioanalyzer质检合格后，使用Illumina Solexa HiSeq 2000进行双端测序(Paired-end)。测序所获原始数据开放存储于figshare(doi：10.6084/m9.figshare.8150792)

1.2.2 测序数据的处理与拼装 测序得到的原始reads并不都是有效reads，在后续组装和分析前，首先对Raw reads进行预处理，移除含poly-N、接头和低质量的片段，获得Clean reads。然后使用Trinity，采用Paired-end的拼接方法，将具有一定长度overlap的reads 连成更长的片段，得到不含未知碱基的组装片段，之后利用TGICL软件聚类去冗余延伸得到一套Unigene，并最终选取每个Loci下最长的转录本作为确定的Unigene。

1.2.3 Unigene序列的同源比对及功能分类 把拼接获得的所有Unigene通过Blast比对(Evalue<1e-5)到主要蛋白数据库中，以获得其通路和功能等方面的注释。采用的数据库包括NR(非冗余蛋白数据库)、Swiss-Prot(蛋白质序列数据库)、GO(基因本体数据库)、KEGG(京都基因与基因组百科全书)和KOG(蛋白质真核直系同源数据库)等，得到跟给定Unigene具有最高序列相似性的蛋白(如有并列，取第1条)，从而得到试验中所获得的Unigene序列的功能注释信息。

1.2.4 转录组中的SSR分析 使用MISA软件(http://pgrc.ipk-gatersleben.de/misa/)搜索所获得的Unigene上的潜在SSR位点。搜索条件为2个碱基的motif(基序)，需要至少6个重复，3个碱基的motif需要至少5个重复，4～6个碱基的motif至少需要4次重复。如果2个SSR序列的距离短于100 bp，就会被认为是同一个潜在SSR标记。此外，为排除polyA结构的干扰，剔除了单核苷酸重复的情况。

2 结果与分析

2.1 测序与Denovo组装结果

基于Paired-end测序获得了大量的样本数据。为排除测序错误率的影响，使用滑动窗口法去除低质量片段，质量阈值20(错误率=1%)，窗口大小5 bp，长度阈值35 bp，并切除reads中含N的部分序列。最终得到了有效RNA-seq数据共0.32亿条，GC平均含量为44.76%，共32 417 032条有效reads，平均长度95.77 bp。对比质控前数据，有效率为87.11%。测序数据质量控制前后对比情况见表1。

表1 测序数据质控前后对比Tab.1 Comparison of sequence data after quality control

注：表中reads数量为双端测序的结果。

Note: reads amount in the table is from pair-end sequencing.

将样本的有效reads合并进行de novo拼接，对拼接序列去重复，最终得到了79 098个长度大于200 bp的转录本，大小73 Mb。取每个Loci下最长的转录本作为Unigene，得到了59 788个Unigene，大小50 Mb，其最大长度、平均长度和N50分别为10 923，842,1 402 bp，长度大于500 bp和1 000 bp的Unigene数目分别为30 171和18 542。图1展示了所有Unigene的长度分布情况。

2.2 Unigene的注释

将所获得的Unigene序列比对到公共数据库(NR、 Swiss-Prot、 KEGG、Go、KOG)，进行序列相似性分析(Evalue

图1 所获得的Unigenes长度分布情况Fig.1 The length distribution of acquired Unigenes

2.3 扁秆荆三棱Unigene的KOG分类

将获得的Unigene与KOG功能分类数据库进行比对，结果表明，有13 334个Unigene 被注释在25种功能分类中(图2)。其中，归类到信号传导机制(Signal transduction mechanisms)的Unigene数量最多(2 836个，4.74%)，其次是蛋白质翻译后修饰与转运、分子伴侣(Posttranslational modification, protein turnover, chaperones)(2 465个，4.12%)及一般功能(General function prediction only)(2 142个，3.58%)。最小的类别是胞外结构(Extracellular structures)(33，<1‰)和细胞运动(Cell motility)(6，<1‰)。

A.RNA 加工与修饰；B.染色质结构与变化；C.能量产生与转化；D.细胞周期调控与分裂，染色体重排；E.氨基酸运输与代谢；F.核苷酸运输与代谢；G.碳水化合物运输与代谢；H.辅酶运输与代谢；I.脂类运输与代谢；J.翻译，核糖体结构与生物合成；K.转录；L.复制、重组与修复；M.胞壁/膜生物发生；O.蛋白质翻译后修饰与转运、分子伴侣；P.无机离子运输与代谢；Q.次生产物合成，运输及代谢；R.一般功能；S.功能未知；T.信号传导机制；U.胞内分泌与膜泡运输；Z.细胞构架。

A.Nucleotide transport and metabolism; 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 biogenesis; O.Posttranslational modification, protein turnover, chaperones; P.Inorganic ion transport and metabolism; Q.Secondary metabolites biosynthesis, transport and catabolism; R.General function prediction only; S.Function unknown; T.Signal transduction mechanism; U.Intracellular trafficking, secretion, and vesicular transport; Z. Cytoskeleton.

图2 Uingenes的KOG分类情况
Fig.2 KOG classification of Unigenes

2.4 扁秆荆三棱Unigene的GO分类

对Unigene进行GO功能分类预测，并有24 670个Unigene被注释上GO分类，并按照基因的分子功能，参与的生物过程和所处的细胞位置进行分类。如图3所示，在生物过程(Biological process)分类中，主要聚集于代谢过程(Metabiolic process)(28.66%)和细胞过程(Cellular process)(23.97%)；在细胞组分(Cellular component)分类中，前三位依次是细胞(Cell)(16.00%)、细胞部分(Cell part)(16.00%)和细胞器(Organelle)(10.32%)；在分子功能(Molecular function)分类中，主要聚集于结合(Binding)(25.85%)和催化活性(Catalytic activity)(23.12%)。

图3 Uingenes的GO分类Fig.3 GO classification of Unigenes

2.5 扁秆荆三棱Unigene的KEGG通路分析

扁秆荆三棱共有13 141个Unigene映射到KEGG数据库中并得到注释，获得的参考代谢通路(Pathway)数目为291。其中包含Unigene最多的是代谢通路是核糖体(ko03010)，共有405条Unigene，其次是剪切体(ko03040)，共有360条Unigene。Pathway富集性分析的前10个Pathway数据列于表2。参与淀粉与糖代谢通路(ko00500)的Unigene有326条，共统计筛选出60条参与淀粉合成的Unigene，编码9个关键酶(表3)，其中5个Unigene编码α-糖苷酶；11个Unigene编码β-呋喃果糖苷酶；13个Unigene编码己糖激酶；6个Unigene编码果糖激酶；4个Unigene编码葡萄糖-6-磷酸异构酶；5个Unigene编码葡萄糖磷酸变位酶；6个Unigene编码葡萄糖-1-磷酸腺苷酰基转移酶；6个Unigene编码淀粉合成酶及4个Unigene编码1，4-α-葡聚糖分支酶，另外参与植物激素信号转导的Unigene为279条，参与根茎伸长生长的赤霉素受体GID1的Unigene有4条，与赤霉素(GA)生物合成过程中的关键酶赤霉素3-β-双氧化酶(GA3ox)同源的Unigene有5条，编码细胞分裂素合酶的Unigene有5条，共计14条Unigenes。

表2 扁秆荆三棱中包含Unigenes数目最多的10个代谢通路Tab.2 Top ten metabolic pathways involving B.planiculmis Unigenes

表3 扁秆荆三棱中地下茎伸长及球茎淀粉合成相关基因Tab.3 The genes and transcription factors related to sprouting and starch synthesis in B.planiculmis

2.6 扁秆荆三棱转录组中潜在SSR位点的分布

利用MISA软件对扁秆荆三棱的59 788条Unigene 序列进行扫描，扫描的碱基容量为50 341 770。结果表明，12 823个Unigene 中含有潜在SSR位点，包含多于1个SSR位点的Unigene有2 531个，发生频率(含有SSR的Unigenes 数量与总Unigenes 数量之比)为21.45%。在剔除单核苷酸重复类型后，总共有9 698个SSR位点，平均0.193个SSR/kbp，SSR的分布频率(SSR的个数与总Unigenes 的数量比)为16.22%，其中二核苷酸motif(6 971个)和三核苷酸motif(2 332个)是主要的串联序列类型，分别占SSR总数的71.88%和24.05%,两者合计占所有SSR 位点的95.93%。所有潜在SSR位点的重复单元详细分布情况见表4。

表4 扁秆荆三棱转录组中潜在SSR位点重复单元的分布Tab.4 Occurrence frequency of different microsatellites motifs in B.planiculmis transcriptome

3 结论与讨论

扁秆荆三棱为稻田和湿地环境中的常见杂草，隶属于莎草科(Cyperaceae)三棱草属(Bolboshoenus)，属于无参考基因组序列的物种，目前尚未发现扁杆荆三棱或其同属物种的转录组研究先例。具体到本研究，在研究策略上，采用茎叶、球茎、根茎混合样品进行转录组测序，有助于在节约试验成本的基础上发掘到更多的转录本[24-25]。在数据质量上，通过高通量测序的技术手段，首次对扁秆荆三棱的转录组信息进行获取和分析。获得了59 788个Unigene，其最大长度为10 923 bp、平均长度为842 bp，N50为1 402 bp(即50%的数据分布在长度不小于1 402 bp的Unigene上)，注释比率达到58.91%。以KEGG的代谢途径数据库为依据，可鉴别出291个代谢通路，并筛选到与地下茎伸长生长相关的Unigene 9条，与地下球茎中淀粉合成相关的Unigene 60条。此外，还获得了潜在的SSR位点9 698个，分布在12 823个Unigene中，发生频率为21.45%，其中二碱基和三碱基的motif为主要类型。类比已有的克隆植物转录组研究，如毛竹(Phyllostachysedulis) (均长612 bp，注释比率69.75%)[26]、巨龙竹(Dendrocalamussinicus)(均长723 bp，注释比率69.22%)[27]，本研究的Unigene长度较为进步，有助于获得较完整的基因表达信息。对于未获得注释的Unigene，可能是由于长度较短无法匹配，也有可能是非编码序列或者是新的基因[28]，需要进一步的试验予以验证。在分析思路上，GO和KOG的功能分类对初步了解基因的功能起着重要作用，而KEGG数据库中的参考pathway不仅可以推测基因的功能，还可以研究基因在不同代谢通路中所在位置及作用，三者相辅相成，成为新物种中发掘功能基因的重要手段[29-30]。本研究通过KEGG 数据库中的pathway分析筛选与淀粉合成相关的基因，同时结合所筛选到的Unigene在相关数据库中予以功能注释，进一步确保了所获得基因的可靠性。

本研究结果有助于理解扁杆荆三棱的地下根(球)茎系统的发展过程。对克隆植物矢竹(Pseudosasajaponica)的研究表明，地下根茎的生长发育受到激素的调节，其中赤霉素在启动地下茎节间快速伸长上具有重要的作用，而细胞壁生长及细胞骨架组织等相关功能基因则在维持节间快速伸长生长上具有重要作用[31]，扁秆荆三棱中也检测到了调控地下茎伸长生长的赤霉素代谢途径相关基因，但尚未发现细胞生长的相关基因。这可能是由于该物种在7月份正处于克隆生长的起始阶段，地下茎系统尚在形成过程中有关。该发现为防治该植物种群的过度扩张提供了新思路。在通路分析中，核糖体、剪接体和内质网上的蛋白加工居于前三位，表明可能在进行较大量的蛋白质合成工作，这与其地下根茎系统的快速扩张需求是相符合的[9, 22]。嘌呤代谢通路也位于前列，作为核苷酸合成、能量物质及一些初级产物如蔗糖、多糖及磷脂类物质的前体物质，嘌呤及其代谢途径在许多生物学过程中具有重要的作用[32]。这也从另一方面印证了该植物地下茎系统的能量与物质需求。此外，本研究的数据也检测到较多淀粉合成相关的基因及通路，与木薯(Manihotesculenta)[33]、马铃薯(Solanumtuberosum)[34]等主要以地下块茎储存淀粉的植物相一致。已有的研究认为扁秆荆三棱球茎中的淀粉是一种“储备”行为，可支持其度过环境胁迫时期，并为球茎萌苗的快速生长提供能量[3, 22]，结果进一步表明这种储备从地下系统扩张的早期就已经开始，为进一步挖掘该植物球茎的饲养价值或野生动物保护价值提供了支持。在研究展望方面，由于本研究采用的是不同组织材料等量混合后测序和组装，虽然有助于得到更多的转录组数据，成本控制也较优秀，但是由于可变剪切等特异性机制[35]，本研究结果在解释组织差异化表达方面具有局限性。在未来的研究中，应注重对该植物不同生长发育阶段，或不同组织的转录组情况进行考察，具体探讨其对水淹、盐渍和种间竞争等生境中主要限制条件的响应机制。