基于RNA-Seq分析茶树新品系‘春绿’生物信息

林郑和，孔祥瑞，陈常颂,2*，李鑫磊，张亚真，游小妹，陈志辉，单睿阳，钟秋生

（1. 福建省农业科学院茶叶研究所/国家茶树改良中心福建分中心，福建 福州 350013；2. 国家土壤质量福安观测实验站，福建 福安 355015）

0 引言

【研究意义】伴随着高通量测序新技术（RNA-Seq）的广泛应用，越来越多的研究者利用转录组学研究茶树生长发育、抗性、物质代谢的分子调控机制，探索茶树基因表达的转录调控机制，并试图通过基因工程技术培育出优良品种。【前人研究进展】昆明植物研究所的研究团队发布了茶树大叶种‘云抗10号’基因组[1]，开启了茶树基因组学的研究序幕，之后安徽农业大学的茶学团队基于三代PacBio测序方法，获得了茶树小叶种‘舒茶早’的基因组序列[2]；华南农业大学也公布了茶树小叶种‘碧云’的染色体级别基因组序列[3]；同时中国农业科学院茶叶研究所也公布了‘龙井43’的染色体级别基因组序列[4]；之后‘黄棪’与‘铁观音’染色体级别基因组序列的公布[5,6]，探索茶树基因表达的转录调控机制，并试图通过基因工程技术培育出优良品种，现已成为茶树研究新兴热点之一。李娜娜等[7]以‘福鼎大白茶’和‘小雪芽’叶片为研究材料，发现小雪芽新梢白化的可能遗传机理存在于蛋白质翻译后的修饰过程，而且这些变异主要发生在细胞内。朱兴正等[8]采用PacBio平台进行‘云茶1号’全长转录组测，发现云茶1号中216个代谢途径分支，包括茶叶品质、活性物质代谢以及抗逆等相关基因。然而，茶树庞大的3 Gb大小基因组、超过70%的重复序列含量以及自交不亲和导致的高度杂合特性[7]，遗传结构十分复杂，遗传分析很困难。【本研究切入点】‘春绿’茶树新品系由福建省农业科学院茶叶研究所于2012—2022年，以‘福云6号’为母本、‘早春毫’为父本杂交选育而成，小乔木型，树姿直立，生长势中等。一芽一叶期3月7日；一芽二叶期3月10日；芽叶绿色，芽叶茸毛中等或较多，一芽三叶长8.2 cm，发芽较密，一芽三叶百芽重88.6 g。叶片为长椭圆形，叶面微隆，叶片质地软，叶尖渐尖，叶基楔形，叶缘微波，叶齿锐，叶齿密度中。春茶一芽二叶蒸青样水浸出物含量49.9%、茶多酚含量22.4%、游离氨基酸总量3.4%、咖啡碱含量4.0%。适制绿茶和白茶，在福建、湖北等多地有种植示范，是一个适制性广、高香且特早生优良茶树品系。【拟解决的关键问题】本研究利用转录组测序技术，分析春绿与福云6号的转录组测序，通过序列分析、基因功能注释、基因功能分类、差异表达基因筛选等，旨在为进一步挖掘春绿茶树新品系次生代谢相关功能基因和抗逆基因奠定基础，为新品种的选育提供参考依据。

1 材料与方法

1.1 试验材料

供试茶树品种福云6号与新品系春绿，种植于福建省农业科学院茶叶研究所二号山同一地块，树龄均为5年，施肥、喷药等茶园管理均一致。春季统一采摘新梢第一轮一芽二叶，迅速用液氮固样，送上海欧易生物医学科技有限公司测序。

1.2 表观性状观测及绿茶制作与审评

表观性状按陈亮等[9]的标准观测。采摘2022年春茶第一批新梢（达到一芽一叶标准），绿茶制作参照传统烘青绿茶加工工艺：鲜叶萎凋→杀青→揉捻→烘干。详细工艺参数为：采用自然萎凋方式（水筛摊凉，1.5斤/个水筛），萎凋时间约6 h；采用滚筒杀青，温度为280℃，投叶量为3斤/锅（2个水筛的鲜叶），时间约90 s。出锅后摊凉10～15 min（待热气散发后），放入6CR-15的揉捻机揉捻10～12 min，迅速解块，毛火120℃约10 min，足火80℃烘箱中烘干至含水率低于5%。

样品感官审评参照GB/T 23776—2018《茶叶感官审评方法》，审评项目有外形、香气、汤色、滋味和叶底五项。审评结果采用加权评分法，汤色10%、外形占20%、滋味30%、香气30%、叶底10%。再结合记录的评审术语，综合评定品质。由4名专业审评人员进行密码审评。

1.3 RNA的提取与文库构建

茶树叶片总RNA的提取参考本实验室改良的Tri-2Reagent法，凝胶电泳检测RNA完整性；测定260 nm和280 nm波长下的OD值，合格后进行混匀。RNA混匀后使用DNase消化DNA，用带有Oligo（dT）的磁珠富集mRNA（去除rRNA）；加入打断试剂将mRNA打断成短片段，用六碱基随机引物合成一链cDNA，然后配制二链合成反应体系合成二链cDNA，并使用试剂盒纯化双链cDNA；纯化的双链cDNA再进行末端修复、加A尾并连接测序接头，然后进行片段大小选择，最后进行PCR扩增；构建好的文库用Agilent 2100 Bioanalyzer质检合格后，使用Illumina HiSeqTM 2500测序仪进行测序，产生125 bp或150 bp的双端数据。具体参照林郑和等[10]的方法。

1.4 Illumina测序、数据质检和参考基因组比对

通过Illumina HiSeq 测序平台，对文库片段进行双末端（Paired-end，PE）测序，得到原始下机数据fastq 文件。使用FastQC 软件的默认参数对fastq 文件进行测序质量统计，使用Cutadapt 软件进行reads的过滤和去除低质量序列（overlap≤10 bp，20%的碱基错误率），最终得到高质量的clean reads。再通过 Bowtie2 和 Tophat2 将 clean reads 比对至茶树参考基因组[2]，得到SAM/BAM 比对结果文件，具体方法参照文献[11]。

1.5 差异表达基因筛选及富集分析

分别以各测序样本的总表达量为内标，对2个测序样本的RPKM进行Fisher-test差异检验。以FDR（错误发现率）≤0.05及Fold-Change（表达差异倍数）≥2为条件进行差异基因筛选，并进行差异转录组的GO和KEGG富集分析，以判定差异转录组主要影响的生物学功能或者通路。

1.6 差异表达基因的实时荧光定量PCR（qPCR）检测

为了验证RNA-Seq数据的结果，随机选取14个ungene进行RT-qPCR分析。RT-qPCR引物对采用 primer Premier 5.0（Premier Biosoft，Palo Alto，CA，USA）设计，内参基因为GAPDH（CSA024857.1），见表1。用RNA纯化试剂盒（中国天根）纯化总RNA，按照说明书使用TaKaRa（中国大连）的 PrimeScript II 1st Strand cDNA synthesis kit进行第一链cDNA的合成。使用 SYBR Premix Ex Taq II（Tli RnaseH Plus）试剂盒（TaKaRa，日本）在Lightcycler 480 Real-Time PCR检测系统（Roche，德国）中进行反应。采用3 μL 模板 cDNA、1 μL 正向引物（5 pmol）、1 μL反向引物（5 pmol）和 5 μL SYBR Green混合物（Qiagen，Hilden，Germany） 进 行Real-time qPCR。PCR过程包括95℃预热5 min，95℃ 10 s，60℃ 20 s，45个循环，然后生成解离曲线（95℃15 s，60℃ 60 s，95℃ 15 s）。每个实验设 3个重复，用2-ΔΔCt法计算各基因的相对表达量。对不同样品进行归一化处理，以GAPDH基因（accession number：CSA024857.1）为内标，以福云6号（CK）的叶片为参比样本，设置其基因表达量为1。具体的方法参照文献[10]。

表1 引物序列及特征Table 1 Primer pair used for RT-qPCR analysis

1.7 数据处理

利用软件SIMCA-P 11.5对所测定的数据进行分析，筛选出差异基因，并对两组样本进行主成分分析（PCA）。

2 结果与分析

2.1 表观性状与品质差异分析

从表2可以发现，春绿叶片呈稍上斜，福云6号的呈水平或者稍下垂，二者有明显的区别；叶身内折也较福云6号大，叶齿锐深与福云6号的钝浅有明显区别，嫩梢茸毛较福云6号少，萌发期大约晚3 d。从表3可看出，新品系春绿香气与滋味得分明显高于福云6号。春绿品质明显优于福云6号。

表2 表观性状分析Table 2 Apparent trait analysis

表3 感官品质审评表Table 3 Sensory evaluation form

2.2 测序原始序列数量及质量

测序共得到春绿与福云6号品种（系）叶片转录组原始序列，通过计算每千个碱基的转录每百万映射RPKM（Reads Per Kilobase of exon model per Million mapped reads）的读取量来分析RNASeq数据的转录谱。构建了12个cDNA文库，包括福云6号叶片（FY6H1-6）与春绿（CL1-6）各6个生物重复，共获得81.51 G的Clean data，各样本的有效数据量分布在6.59～7.2 G，Q30 碱基分布在89.6%～90.86%，从每个库生成的原始读取数从47 385 674到51 271 456。Clean reads和Q20（测序错误率低于1%）的百分比分别超过97%和95%（表4），碱基G和C数占总碱基数的44%以上，平均GC含量为45.00%。通过将reads比对到参考基因组上，得到各个样本的基因组比对情况，比对率为90.91%～91.52%，不能测序的核苷酸N为0，说明转录组测序质量整体较高。

表4 测序数据统计结果Table 4 Summary of RNA-Seq data on Fuyun No. 6 (FY6H) and Chunlv (CL)

2.3 样本的PCA分析与差异基因分析

从图1-A可以看出，该组原始试验数据所得在两个主成分PC1、PC2中良好地呈现，第一主成分的贡献率为99.26%，第二主成分的贡献率为0.22%，两种贡献率和为99.46%。通过样品分布散点图发现，2个品种的信息可以得到很好的区分，说明拟合性良好，能够直观地反映样本间的差异性。通过对2个品种的基因测序，以福云6号为对照，在错误发现率（false discovery rate，FDR）＜0.05且|log2 fold change（FC）|＞1的筛选条件下，发现春绿中有7 042个差异基因（DEGs），其中3 273个表达上调，3 769个表达下调，见图1-B。

图1 春绿与福云6号PCA分析（A）与差异基因分析（B）Fig. 1 Principal component and differential gene analyses on CL and FY6H

2.4 GO功能富集分析

春绿与福云6号样本的DEGs主要分为三类，包括23个基于生物过程的GO组、20个基于细胞成分的GO组和21个基于分子功能的GO组（图2）。通过GO（Gene Ontology）富集分析表明，差异显著的生物过程功能基因中（前10）有68条基因显著上调、268条显著下调；参与细胞功能组成基因中有67条显著差异基因上调、1 444条显著差异下调；参与分子功能基因中发现220条显著差异基因上调、623显著差异基因下调。以上发现，参与细胞功能组成差异基因最多。GO项富集分析显示，20个DEGs在生物碱代谢过程（GO：0009820）富集，26个DEGs在营养活动（GO：0045735）中富集，393个DEGs在ATP结合中富集（GO：0005524），13个DEGs富集在次生代谢过程（GO：0019748），3个DEGs富集在精胺酸分解过程(GO：0046208)，28个DEGs富集在葡糖基转移酶活性（GO：0035251），3个DEGs富集在苯甲酸葡萄糖基转移酶活性(GO：0052641)，8个DEGs富集在水杨酸糖苷转移酶(GO：0052640)。28个DEGs富集在葡糖基转移酶活性。

图2 春绿与福云6号差异显著富集GO分类Fig. 2 Significantly enriched GO groups in CL and FY6H

2.5 KEGG代谢通路分析

对春绿与福云6号样本进行KEGG代谢通路分析（图3），有8 657个基因注释到191个代谢通路中，这些基因中有1 277个显著差异基因。其中Unigenes排名前20代谢通路有：40条DEGs参与卵母细胞减数分裂（ko04114），45条DEGs参与RNA降解（ko03018），29条DEGs参与脂肪酸代谢（ko01212），17条DEGs参与不饱和脂肪酸的生物合成（ko01040），40条DEGs参与药物代谢-细胞色素P450（ko00982），40条DEGs参与细胞色素P450外源生物代谢（ko00980），27条DEGs参与氨酰基-tRNA生物合成（ko00970），66条DEGs参与苯丙素的生物合成（ko00940），13条DEGs参与氮代谢（ko00910），13条DEGs参与类胡萝卜素生物合成（ko00906），11条DEGs参与油菜素甾醇生物合成（ko00905），10条DEGs参与生物素代谢（ko00780），37条DEGs参与光合作用生物的碳固定（ko00710），27条DEGs参与甲烷代谢（ko00680），8条DEGs参与苯乙烯退化（ko00643），5条DEGs参与氨基苯甲酸酯退化代谢（ko00627），53条DEGs参与谷胱甘肽代谢（ko00480），14条DEGs参与氨酸代谢（ko00380），23条DEGs参与精氨酸和脯氨酸代谢（ko00330），14条DEGs参与角质、亚伯碱和蜡的生物合成（ko00073）。

图3 春绿与福云6号KEGG代谢通路富集显著差异基因表达模式Fig. 3 Top enriched KEGG pathways of DEGs in CL and FY6H

2.6 差异表达基因的qPCR检测

为验证转录组测序结果的可靠性，从转录组数据鉴定到的差异基因中选取了14个基因进行实时荧光定量PCR分析。其将RT-qPCR得到的基因表达趋势与转录组测序结果进行比较，发现基因的表达趋势基本是一致的（图4），说明转录组数据是可靠的。

图4 差异基因的表达验证Fig. 4 Expression verification of differential genes

3 讨论与结论

近年来，随着基因组测序技术的进步，茶树研究已经步入基因组时代。首先是大叶种茶树云抗10号的基因组序列，拉开了茶树基因组学研究的帷幕[1]，之后基于三代测序补洞的方法测序获得了茶树小叶种品种舒茶早的基因组序列[2]，接着乌龙茶品种黄棪[5]和铁观音[6]的染色体级别基因组发布等。传统的茶树育种和改良方式（新品种选育需15～20年）已很难满足目前市场多元化的需求[12]，分子生物学技术介入茶树育种，能加速理想性状的茶树种质的分子育种。本研究通过对茶树（春绿与福云6号）进行转录组测序，获得了大量转录因子Unigene，可为后期茶树新品系（春绿）抗性研究、次生代谢研究、新品种选育等提供理论依据。

本研究通过测序分别获得春绿为47 367 063.33条Clean reads，福云6号为48 130 043.67条Clean reads，福云6号比春绿有着更长的序列，通过将reads 比对到参考基因组上，得到各个样本的基因组比对情况，比对率为90.91%～91.52%，说明转录组测序质量整体较高。通过对差异表达基因的GO富集分析表明，新品系春绿与福云6号在生物过程、细胞组分和分子功能三大类别上均分布有差异基因。其中393个DEGs在ATP结合中富集（GO：0005524），28个DEGs富集在葡糖基转移酶活性（GO：0035251），26个DEGs在营养活动(GO：0045735)中富集，20个DEGs在生物碱代谢过程（GO：0009820）富集。进一步分析发现，新品系春绿中66条DEGs参与苯丙素的生物合成（ko00940），其中参与苯丙氨酸与酪氨酸代谢中的大部分基因表达都增加，其中苯基丙酸合成基因（LOC114255969）上调2.49倍、甘露醇脱氢酶（LOC114256196）上调3.2倍、β-葡萄糖苷酶（LOC114256477）上调4.1倍。说明新品系春绿在氨基酸合成、香气合成等方面明显强于福云6号。这也进一步说明了新品系春绿制绿茶品质明显优于福云6号。

脂肪酸代谢由涉及脂肪酸或与脂肪酸密切相关的各种代谢过程组成，脂肪酸是属于脂质常量营养素类别的分子家族。这些过程主要可分为：（1）产生能量的分解代谢过程；（2）合成代谢过程，它们作为其他化合物的构建块。本研究发现，29条DEGs参与脂肪酸代谢（ko01212），大部分都显著上调，其中酰基辅酶A氧化酶（LOC114262883）上调3.1倍；ACP脱氢酶（LOC114265912）上调3.6倍。结合GO富集分析，发现393个DEGs在ATP结合中富集（GO：0005524），说明新品系春绿能量代谢明显比福云6号强。KEGG分析还发现，14条DEGs参与角质、亚伯碱和蜡的生物合成（ko00073），其中蜡粉基因CER1（LOC114261760）上调3.2倍，有研究发现CER1基因编码醛脱羰酶，在烷烃生物合成通路中催化醛脱羰形成烷烃，其突变体中烷烃的含量显著减少，而过表达CER1基因时，烷烃的含量会增加，器官呈现蜡粉合成减少[13]。还发现羟基棕榈酸酯-O-阿魏酰转移酶（LOC114285233、LOC11 4287670）表达下调，羟基棕榈酸酯-O-阿魏酰转移酶主要与萌芽相关，福云6号萌发期相对较早，这与表观性状观测结果基本一致。

本研究通过差异转录因子的聚类分析，有助于从分子水平进一步研究新品系春绿与福云6号的品质、性状差异，建立自己优选品种的基因表达信息，为茶树的本土育种和改良建立分子信息库，这对提高茶叶品质和产量具有重要意义。