滇白前种子转录组测序及功能注释

王建秋， 王晓丽， 曹子林

(西南林业大学， 云南 昆明 650224)

滇白前(SileneviscidulaFranch)又名黏萼蝇子草，为石竹科(Cayophyllaceae)蝇子草属的一种多年生草本植物，主要生长在海拔1 200～3 200 m的四川、云南、贵州、西藏(东南)等地[1]。据研究，滇白前对铅、锌、镉具有很强的耐受性，是一种新的Pb/Zn/Cd共超富集植物，为土壤重金属复合污染地修复提供了新的种质资源[2]。目前，有关滇白前的研究主要集中于药用成分及抗癌作用上[3-4]。种子是植物生命周期下一代的载体[5]，对于这一具有重要药用价值的重金属超富集植物，其种子基因组结构及其功能是否有其特殊性未见报道。转录组是特定的细胞或组织在某一发育阶段或功能状态下转录出的所有RNA总和[6]，功能基因研究的主要方法是基于高通量测序的转录组分析[7]。转录组测序分析不仅能迅速获得基因转录组信息、功能注释和样品间的差异表达，而且能揭示特定生物学过程下，器官、细胞的分子机制[8]。本研究拟进行滇白前种子的转录组测序，分析测序质量并进行基因功能注释及转录因子预测，为深入挖掘滇白前的功能基因及分子机制等方面的研究奠定一定的基础。

1 材料与方法

1.1 样品采集

样品于2020年1月1日采自云南省怒江州兰坪县金顶镇凤凰山金顶铅锌矿区成熟的滇白前种子，地处北纬26°21′～27°02′，东经98°23′～99°28′，此矿区具有储量大、品位高、储存集中、埋藏浅等特点，是亚洲最大的铅锌矿床[2]。

1.2 RNA提取、文库的构建及测序

从种子中提取总RNA后，利用带有Oligo(dT)的磁珠与ployA进行A-T碱基配对，从总RNA中分离出mRNA；富集得到完整的RNA序列，通过加入fragmentation buffer,将mRNA随机断裂、磁珠筛选分离出300 bp左右的小片段；在逆转录酶作用下，加入六碱基随机引物，以mRNA为模板反转录合成一链cDNA，随后进行二链合成，形成稳定的双链结构；连接adaptor，用2%的琼脂糖胶回收目的条带大小合适的片段并通过PCR技术进行扩增；最后，用Illumina Novaseq 6000对短序列片段进行测序。

1.3 测序原始数据过滤及组装

测序获得原始数据后，使用SeqPrep软件将原始数据中含有大量的接头序列、低质量序列末端以及adapter和N的比例大于10%的碱基去除，从而得到高质量的质控数据(clean data)。使用Trinity软件[9]对所有clean data进行从头组装，使用TransRate[10]对从头拼接结果常见的错误(包括嵌合体、结构错误、组装不完整、碱基错误等)进行评估，给出每条contig的质量得分，并可将这些contigs评分整合获得整个组装结果的综合分数，对转录组获得的序列进行过滤和优化。使用CD-HIT软件[11]通过序列比对聚类(Cluster)的方法去除冗余、相似的序列，最后获得非冗余(non-redundant)的序列。BUSCO[12]利用单拷贝直系同源基因，评估基因组或转录组的组装完整性。基因组组装是与BUSCO的一致序列进行tBLASTn比对，然后再用Augustus对基因结构进行预测，而转录组组装则先寻找转录本的ORF编码框，然后再进行HMMER 3比对。

1.4 基因功能注释

用NCBI将Unigene序列比对到蛋白数据库Nr、SwissProt和COG(evalue<0.000 01)，用KOBAST软件将Unigene序列进行KEGG注释，用BLAST 2 GO软件将Unigene序列进行GO注释[13]。得到与给定Unigene具有最高序列相似性的蛋白，从而获得该Unigene的蛋白功能注释信息，并将注释结果进行分类统计。将预测的蛋白序列同相应的TF数据库(plant TFdb/animal TFdb)进行hmmscan比对，得到相应的转录因子家族。

2 结果与分析

2.1 滇白前种子转录组组装质量分析

对滇白前种子进行转录组测序后，通过严格质量控制和数据过滤，共获得50 423 900个高质量干净序列。利用Trinity软件，对这些高质量数据进行denovo组装，共获得有效序列片段32 916 084条，N 50为1 192 bp，平均长度为753.87 bp，GC核苷酸含量(%)为43.45%的滇白前Unigene 43 663条(不含N的组装片段)，长度分布于201～15 324 bp之间(表1)。通过用长度分布柱状图对组装后的基因做统计(图1)，长度为200～500 bp的有23 651条Unigene(54%)，长度为501～1 000 bp的有9 625条Unigene(22%)，长度为1 001～2 000 bp的有4 606条Unigene(11%)，长度为2 001～3 000 bp的有1 407条Unigene(3%)，长度超过3 000 bp的有797条Unigene(1%)。由此可得，随着基因长度的增加，基因数量呈下降趋势。

表1 Unigene基本信息Table 1 The basic information Unigenes

2.2 滇白前基因的总体注释情况

将比对到植物中的基因与公共数据库进行功能注释，分别有26 688、20 818、25 657、19 942、14 192条Unigenes比对到Nr、Swissprot、COG、GO、KEGG数据库，至少有29 177条Unigenes比对到一个数据库，并注释到功能，有14 486条Unigenes没有比对到数据库，说明有许多功能不清的基因。其中，Nr数据库注释到的基因最多，共有26 688条，占比为61.12%；KEGG数据库注释到的基因最少，共有11 527条，占比为26.4%(表2)。

表2 Unigene注释统计Table 2 Annotation statistics of Unigenes

2.3 滇白前基因的NR功能注释

共有26 688条Unigenes比对到NR数据库，与其他物种比对后的结果如图2。注释序列物种中注释到基因数量最多的是甜菜(Betavulgaris)，共有6 828条，占比为25.58%；其余依次为藜麦(Chenopodiumquinoa)、栓皮栎(Quercussuber)、菠菜(Spinaciaoleracea)，分别有5 997、4 738、4 272条Unigenes，占比分别为22.47%、17.75%、16.01%；其他289个物种共有4 853条Unigenes，占18.18%。

2.4 COG功能注释

将滇白前基因比对到COG数据库发现，有26 500条注释到了COG数据库中，占总Unigenes的61.12%，这些Unigenes被分为23类(图3)。有13 917个功能注释信息未知，未确定其准确的生物学功能，占所有功能注释信息的52.51%；翻译后修饰与转运注释到1 807条Unigenes，所占比例为6.82%；翻译、核糖体结构与生物合成注释到1 320条Unigenes，所占比例为4.98%；转录注释到1 238条Unigenes，所占比例为4.67%；胞内运输、分泌和囊泡运输注释到1 224条Unigenes，所占比例为4.62%；信号传导机制注释到1 073条Unigenes，占4.05%；最少的是细胞能动性和核酸结构的Unigenes，分别仅有7条(0.03%)和3条(0.01%)。这些结果表明，滇白前在蛋白质翻译后修饰与转运、翻译、核糖体结构与生物合成、转录、胞内运输、分泌和囊泡运输、信号转导机制等基因表达丰度较高。

2.5 GO功能注释

一共有19 942条Unigenes被注释到了GO数据库中，得到79 481个功能注释，全面性地描述了不同生物基因的生物学特征，共分为三大功能，分别是生物过程、细胞组分和分子功能，依次得到24 770、30 335、24 376个功能注释，占比分别为31.16%、38.17%、30.67%。表3中展示的是注释上的GO二级分类条目情况。3个本体细分为52个功能亚类，生物学过程类占其中23个功能亚类，细胞进程(8 260个注释)和代谢过程(7 051个注释)占比较大，其次是生物调控(2 466个注释)；分子功能类包括16个功能亚类，结合(10 327个注释)所占比例最多，其次是催化活性(10 240个注释)；细胞组分类包括13个亚类，“细胞”所占比例最多(9 017个注释)，其次是“膜”(6 205个注释)和“细胞器”(5 066个注释)。只有少数转录本被注释参与细胞杀死、行为、蛋白质标签、翻译调节器活动、营养物储存活动及细胞外基质。结果表明，以细胞过程、代谢过程、结合和催化活性相关的基因较多。对GO功能注释进一步分析，梳理得到许多抗逆性相关注释，分别是环境应激响应、蛋白激酶、激素代谢、渗透调节、氧化还原反应、受体和离子转运相关的注释(表4)。

表4 滇白前转录组中与逆境胁迫相关的基因Table 4 Adversity stress-related genes in the transcriptome of S. viscidula

表3 滇白前GO注释结果Table 3 GO Annotation of S. viscidula

2.6 KEGG功能注释

有11 527个Unigenes比对到KEGG数据库的参考代谢通路中，可将滇白前Unigenes归为6个类别，代谢相关的通路Unigenes有5 155条，占比最多为44.72%；其次是遗传信息处理相关的通路Unigenes有3 582条，占比为31.07%；环境信息处理相关的通路Unigenes有491条，占比为4.26%；细胞过程相关的通路Unigenes有968条，占比为8.40%；生物体系统相关的通路Unigenes有794条，占比为6.89%；而人类疾病相关的通路Unigenes有537条，占比为4.66%。代谢相关的通路细致分成11个亚类，其中排名前三的是碳水化合物代谢、氨基酸代谢、能量代谢，占比分别为11.46%、7.37%、6.03%，遗传信息加工和生物系统各分为4个亚类，环境信息处理和细胞过程各分为两类，人类疾病相关的通路分为3个亚类，其中“翻译”、“折叠、分类和降解”、“运输和代谢”以及“环境适应性”分别占15.47%、9.33%、6.11%、5.66%，其余11个亚类所占比例均小于5%(表5)。

表5 滇白前Unigenes的KEGG注释结果及分类Table 5 KEGG annotation and classification of Unigenes in S. viscidula

2.7 转录组Unigene的转录因子预测

滇白前转录组共预测到527个转录因子，可分为32个转录因子大家族。其中，MYB(65个，12.33%)属于最大家族，其次是ERF(54个，10.25%)，接着为bHLH(43个，8.16%)、C 2 C 2(38个，7.21%)、bZIP(35个，6.64%)、B 3(35个，6.64%)、C 3 H(29个，5.5%)、GRAS(24个，4.55%)、WRKY(24个，4.55%)和LBD(24个，4.55%)，其他13个家族共有基因132个，占25.05%(图4)。

2.8 SSR分析结果统计

采用MISA基于拼接所得转录本序列信息进行SSR分析(表6)，筛选得到43 663个Unigenes，识别的SSRs总数为4 730个，包含SSR的序列数目为3 995个，共有6种SSR重复类型，其中单碱基重复SSR 2 554条，占54%；二碱基重复SSR 459条，占9.7%；三碱基重复SSR 1 538条，占32.53%；四碱基重复(67条)、五碱基重复(22条)和六碱基重复(90条)共仅占3.78%。单碱基重复SSR中，发生频率最高的是A；二碱基重复SSR中，发生频率最高的是AG，最低的是CG；三碱基重复SSR中，发生频率最高的是ATC和AAG，频率最低的是ACG。

表6 滇白前SSR分析Table 6 Statistics of SSR analysis of S. viscidula

3 结论与讨论

高通量测序技术深受研究者欢迎，目前该技术广泛应用到植物特殊功能基因的挖掘与鉴定[14]。本研究利用Illumina Novaseq 6 000对滇白前种子进行转录组测序，产生clean reads的Q 20值为97.79%，Q 30值为93.25%以上。经de nove拼接组装及去冗余处理后，共得到43 663条Unigene，最大长度为15 324 bp，最小长度为201 bp，N 50为1 192 bp，平均长度为753.87 bp。组装结果明显高于白刺花(Sophoraviciifolia)种子(N 50为537 bp，平均长度为462 bp)[15]，而与油菜(Brassicanapus)(N 50为898 bp，平均长度659 bp)较为接近[16]。

COG、GO和KEGG注释是基因注释的重要途径，GO的功能分类对了解基因功能有重要作用，而KEGG数据库中的参考通路不仅可以推测基因功能，而且可以研究基因在不同代谢通路中的位置及作用[17]。通过COG功能分类可知，参与蛋白质翻译后修饰与转运、翻译、核糖体结构与生物合成的Unigene最多，分别是1 807条、1 320条，发现有13 917条未知功能基因，这有可能与序列片段过短、注释信息缺乏和获得新基因等有关系[18]。与COG数据库的比对，为了解基因功能及进化研究奠定了一定的基础[19]。GO功能共分为3个大类和52个亚类，其中参与细胞过程和代谢过程的Unigene数目最多，分别是8 260条和7 051条，并深入挖掘到从属于环境应激响应、蛋白激酶、植物内源激素、渗透调节、氧化还原、受体和离子转运的与植物抗逆相关的Unigenes。通过KEGG数据库和通路分析发现，涉及代谢途径的Unigene(5 155条)最多，通路数量最多的也是涉及代谢途径的，共101条。其中，嘌呤代谢(298)和嘧啶代谢(231)是核苷酸代谢的核心，被视为植物界内的管家功能，可以通过激活ABA代谢途径从而发挥对非生物逆境胁迫条件的应激保护作用[20]。植物激素信号转导代谢途径在KEGG中注释到196个Unigenes，植物在非生物逆境胁迫条件下，积累ABA，通过激素激活应激相关的基因参与基因表达和信号转导，从而获得承受恶劣环境的耐受能力[21]。

转录因子是读取并解释DNA中遗传“蓝图”的蛋白质组之一，其与DNA结合，可以帮助启动一个增加或减少基因转录的程序，因此转录因子对于许多细胞过程是至关重要的。转录因子在胁迫刺激下不断合成，将信号传递和放大，并通过与启动子特定区域相结合调控下游响应逆境相关基因的表达，从而引起植物生理生化的改变，在植物多种非生物逆境胁迫中起着关键作用[22]。根据转录因子预测的结果，滇白前中MYB_superfamily转录因子最多，其次是AP 2/ERF，接着是bHLH、C 2 C 2、bZIP、B 3、WRKY、NAC、C 2 H 2等转录因子。据研究，植物中的MYB、bZIP、WRKY、NAC、AP 2/ERF、bHLH、C 2 H 2等各类转录因子在响应高盐、干旱、寒冷等逆境胁迫中均有重要作用[23-28]。对这些转录因子进行统计,大约有261个(49.5%)可能与其抗逆适应相关。这也许是滇白前能在重金属浓度很高的土壤中正常生长，也是其能分布在3 200 m高海拔地方的原因[1-2]。简单序列重复SSR标记是目前比较理想的分子标记，滇白前转录组测定共筛选了3 995个SSRs序列，分为6种SSR重复类型，占比较多的是单碱基重复和三碱基重复。这为后续开展滇白前生长和抗逆分子机制研究提供了有力依据。

云南省素有“有色金属王国”之美誉，土壤重金属污染已成为有色金属开采行业面临的严重问题，而且污染土壤通常表现为多种重金属复合的特点[2]。滇白前是铅、锌、镉共超富集植物，是一种修复土壤重金属复合污染的宝贵的物种资源[2]。本研究以滇白前种子为研究材料，利用高通量测序技术构建了转录组数据库，得到了许多的转录本序列信息，并进行了功能注释和分类、转录因子和SSR等分析，揭示了滇白前种子转录组的整体特征。这为研究滇白前的繁殖机制、环境适应性机制等提供了依据。