基于高通量测序的都匀地区福鼎大白种茶树根茎叶分析

王芬 裴会敏 文狄 陈志 刘荣 姚玉仙 马媛

摘 要：为探究茶树中茶多酚等产物代谢途径的相关基因，该研究以贵州都匀地区福鼎大白种茶树的根茎叶为对象，利用高通量测序技术构建茶的转录组数据库并筛选其根茎叶差异表达基因。结果表明：共获得70.88 Gb Clean Data，各样品Clean Data均达到6.33 Gb，Q30碱基百分比在93.22%以上。将Clean Reads与中国种茶树参考基因组进行序列比对，比对效率从87.83%到91.14%。基于比对结果，进行可变剪接预测分析和基因结构优化分析，发掘新基因13 531个，其中10 244个得到功能注释。利用FPKM进行基因表达量分析，根据基因在不同样品中表达量识别差异表达基因。叶与茎的差异基因有5 595个，其中2 769个在茎中上调，2 826个下调，叶与根有9 650个差异基因，5 056个上调，4 594个下调，茎与根中有5 644个差异基因，2 938个上调，2 706个下调，并通过GO和KEGG分析，将差异基因进行功能注释和富集分析。上述结果为揭示都匀地区福鼎大白种茶参与类黄酮、茶氨酸和咖啡碱等代谢途径相关的基因提供了参考，为选育优良品种等提供了理论依据。

关键词：都匀地区，福鼎大白种，根茎叶，高通量测序，差异基因

中图分类号：Q945

文献标识码：A

文章编号：1000-3142（2020）09-1269-12

Abstract：Tea tree is rich in catechins，theanine，caffeine and other metobolite of health fuction. In order to study the related genes of the metabolisms of the polyphenols. We use high-throughput sequencing technology to study the root，stem and leaf of Fudingdabai tea and find differential expression genes （DEGs）. The results showed that 70.88 Gb Clean Data was obtained，6.33 Gb Clean Data is in each sample and Q30 is more than 93.22%. We map the Clean Reads to reference genome，the blast result is from 87.83% to 91.14%. Then，alternative splicing and gene structure optimization was analyzed. There are 13 531 new genes，in which，10 244 genes were annotated. GO and KEGG functional annotation and enrichment analysis were carried out in differential expression genes，which were identified according to gene expression level in different samples. There were 5 595 DEGs between leaf and stem，2 769 genes were up-regulated and 2 826 genes were down-regulated. 9 650 DEGs were found beween leaf and root，5 056 genes were up-regulated and 4 594 genes were down-regulated. 5 644 DEGs between stem and root，2 938 genes were up-regulated and 2 706 genes were down-regulated. The results are expected to provide reference for recognizing genes of catechins，theanine，caffeine pathways，provide the theoretical basis for breeding improved seeds.

Key words：Duyun，Fudingdabai，root，stem，leaf，high-throughput sequencing，differential expression gene

都匀地处云贵高原苗岭山脉南侧，具有低纬度、高海拔、寡日照、多云雾的自然生态环境，适宜茶树生长，资源丰富（陈世军等，2017）。不同品种的适制性、抗逆性和生长势不同，从而导致茶品质和产量的不同。20世纪90年代以来，都匀市等地开始大量引种单产高、品质好、适应性强的福鼎大白种茶树，增加了茶树的遗传多样性。目前，关于都匀地区的茶种质资源保护和抗性研究薄弱（张丽娟等，2016），特别是在转录组层面上的研究较少，在一定程度上限制了茶产业的应用发展。福鼎大白种茶园占全市茶园面积85%以上，因此开展都匀地区福鼎大白种根、茎、叶的转录组研究对茶树在品质、质量及抗逆相关的遗传改良上具有重要意义。

由于高通量测序技术能快速全面获取研究对象在某一状态下基因转录信息，因此该技术广泛应用于生物体转录组基因表达分析，能准确发掘重要功能基因。Shi et al.（2011）用RNA-seq技术对龙井43的嫩叶、成熟叶、茎、幼根、花蕾以及成熟种子进行测序，总共获得127 094个基因，其中包括与茶特有香气相关的代谢产物如咖啡碱、茶氨酸、黄酮等合成相关的基因。王君雅等（2019）利用转录组技术对龙井43和中茶126进行差异基因表达分析，结果表明99个表达模式完全相反的基因参与的生物过程主要有MAPK信号通路、谷胱甘肽和苯丙烷等。利用转录组测序技术可识別编码次生代谢的调节基因和结构基因，也可用来预测未知基因的功能，是进行基因组功能研究的重要工具（Novaes et al.，2008）。本研究利用高通量测序技术，对都匀地区福鼎大白种茶树根茎叶中类黄酮相关代谢途径进行研究，探索根茎叶中差异表达基因参与的代谢途径，为探讨都匀地区福鼎大白种茶树生长发育及组织间功能差异的分子机制提供理论依据。

1 材料与方法

1.1 材料

供试材料选用黔南州都匀地区种植的茶树-福鼎大白种。茶苗选用黔南民族师范学院茶园的扦插苗。选取9株培养条件相同长势相同的茶苗分为3组，每组3株，在第1组中取嫩根、嫩茎和嫩叶分别放入3个离心管作为第1组生物学重复的样本，分别标为根的第1个生物学重复，茎的第1个生物学重复，叶的第1个生物学重复。每个样本3个生物学重复，剩余的两组生物学重复分别在剩下的两组茶苗中取样，共9个样本，分别为叶1、叶2、叶3、茎1、茎2、茎3、根1、根2和根3。采样时在茶树上剪取发育阶段相似、部位相同、充分伸展的叶片、茎和根，样品随后迅速放入液氮罐中速冻后放入-80 ℃冰箱中储存备用。

1.2 转录组测序

首先将采集的福鼎大白种茶苗根茎叶样品放在干冰中送往北京百迈客生物科技有限公司进行二代Illumina高通量测序，对测序数据进行质量控制，并将二代得到的转录本与已测序的中国种茶树基因组通过序列比对的方法进行比较。然后进行转录组文库质量评估、SNP/InDel分析、可变剪接事件预测、差异基因分析和差异外显子分析等。

Illumina二代测序数据提交至NCBI的SRA数据库，BioProject的编号为PRJNA562747。

2 结果与分析

2.1 测序数据质量分析

从茶苗上分别取福鼎大白的根茎叶作为实验材料进行转录组测序，碱基质量值达到Q30的在93.22%以上，GC含量为44.56%～45.85%（表1），结果表明转录组文库质量符合分析要求。

2.2 转录组与中国种茶树基因组比对

2.2.1 比对效率 本研究利用Hisat2（Kim et al.，2015）对测序数据与中国种茶树基因组（Wei et al.，2018）进行比对，利用StringTie（Pertea et al.，2015）对比对上的Reads进行组装和定量，比对到参考基因组上的Reads占Clean Reads的百分比为87.83%～91.14%（表2）。

2.2.2 SNP/InDel分析 利用GATK（Mckenna et al.，2010）软件对Hisat2比对结果中的SNP位点和InDel进行识别，进而分析基因表达水平和基因功能，SNP分析表明基因区SNP位点数多于基因间区，转换型SNP多于颠换型SNP（表3）。InDel主要存在于内含子区和基因间区（图1）。

2.2.3 可变剪接预测 基因通过转录生成前体mRNA，再经过不同的剪接，产生不同的成熟mRNA，翻译为不同的蛋白质。利用Asprofil（Florea1. 基因间; 2. 基因内; 3. 内含子; 4. 上游; 5. 下游; 6. 受体剪切位点; 7. 供体剪切位点; 8. 剪切位点; 9. 起点缺失; 10. 移码; 11. 密码子缺失; 12. 密码子插入; 13. 密码子改变和缺失; 14. 密码子改变和插入; 15. 同义编码区; 16. 非同义编码区; 17. 同义终止区; 18. 终止区; 19. 终止缺失; 20. 其他。

1. Intergenic; 2. Intragenic; 3. Intron; 4. Upstream; 5. Downstream; 6. Splice site acceptor; 7. Splice site donor; 8. Splice site region; 9. Start lost; 10. Frame shift; 11. Condon deletion; 12. Condon insertion; 13. Condon change and condon deletion; 14. Condon change and condon insertion; 15. Synonymous coding; 16. Non synonymous; 17. Synonymous stop; 18. Stop gained; 19. Stop lost; 20. Other.et al.，2013）軟件对可变剪接类型和表达量进行分析，结果表明在9个样品中5′端和3′端外显子可变剪切最多（图2）。

2.2.4 基因结构优化 由于处理数据的软件不同或数据本身的局限性，本文对中国种茶树基因组的261个基因结构进行了优化并列出了部分优化的基因（表4）。

2.3 新基因

2.3.1 新基因功能注释 我们对利用StringTie拼接出来的转录本与中国种茶树基因组的注释信息进行比较，发掘出了13 531个新基因，并将它们与Swiss-Prot，COG，Pfam，KEGG，GO和NR数据库进行BLAST，获得新基因的注释信息（表5）。

2.3.2 基因表达量 本文通过最大流量算法，采用FPKM （Fragments Per Kilobase of transcript per1. 可变5′或3′端剪切; 2. 单内含子滞留; 3. 多内含子滞留; 4. 多外显子跳跃; 5. 单外显子跳跃; 6. 第一个外显子可变剪切; 7. 最后一个外显子可变剪切; 8. 可变5′或3′端剪切（模糊边界）; 9. 单内含子滞留（模糊边界）; 10. 多内含子滞留（模糊边界）; 11. 多外显子跳跃（模糊边界）; 12. 单外显子跳跃（模糊边界）。

1. Ae; 2. Ir; 3. Mir; 4. Mskip; 5. Skip; 6. Tss; 7. Tts; 8. Xae; 9. Xir; 10. Xmir; 11. Xmskip; 12. Xskip.

福鼎大白种茶根茎叶三个部位的基因表达量的密度分布对比（图3）表明，能够测序到的编码蛋白的基因FPKM值跨越10-2到104六个数量级。

2.4 差异表达分析

2.4.1 差异表达筛选 利用Deseq（Wang et al.，2010）进行样品间的差异表达分析，并且将Fold Change≥2，FDR<0.01作为筛选标准，该文列出了部分的差异表达基因并且统计了根茎叶两两样品间的差异基因的数目（表6，表7）。从叶和茎、叶和根、茎和根的差异表达火山图（图4）中可以看出，差异基因的统计学显著性以及上调基因和下调基因，红点代表上调，绿点代表下调。此外，我们还计算了三个组织差异表达基因的聚类热图（图5）。

2.4.2 差异基因功能注释和富集分析 对差异基因进行COG，GO，KEGG，KOG，NR，Pfam，Swiss-Prot和eggNOG功能注释（表8），并且绘制了叶和茎、 叶和根、 茎和根差异基因的GO功能分类图（图6）。应用超几何检验，对差异表达基因进行KEGG通路富集分析，统计了显著性q值最小的前20个通路（图7），图中每个圆圈代表一个pathway通路，圆圈大小代表通路总富集的基因数目，圆圈越大，表示基因越多。圆圈颜色代表q值，q值越小，代表差异基因在该通路中的富集性越可靠，所以说该图越靠近右下角的圆圈，参考价值越大。

2.5 DEU分析

本文应用DEXSeq（Anders et al.，2012）进行外显子水平的差异分析，即DEU（differential exon usage）分析，設置FDR<0.01，该文列出了叶和茎的部分DEU（表9）。

3 讨论

福鼎大白种茶树是异花授粉的经济作物，基因组杂合度高，具有优良的发芽率、抗寒性和耐旱性（刘本英等，2008）。本研究对都匀地区福鼎大白种茶根茎叶进行了转录组测序，Q30碱基百分比不小于93.22%，样品与参考基因组的比对效率在87.83%～91.14%之间。李明玺等（2018）对静安白茶芽和叶的转录组研究中，Q30值都为87.32%，Liu et al.（2017）对瑞雪的转录组研究中，Q20为97%左右，Wei et al.（2018）对龙井43的转录组研究中，Q30为90%左右，因此，福鼎大白种1. 代谢过程; 2. 细胞过程; 3. 单组织过程; 4. 生物调节; 5. 应激反应; 6. 定位; 7. 细胞成分组织或生物合成; 8. 发育过程; 9. 多细胞生物过程; 10. 信号; 11. 生殖; 12. 生殖过程; 13. 多组织过程; 14. 生长; 15. 解毒作用; 16. 免疫系统过程; 17. 生物附着; 18. 生物相; 19. 节律性过程; 20. 细胞杀伤性; 21. 运动; 22. 细胞; 23. 细胞组件; 24. 细胞器; 25. 细胞膜; 26. 细胞膜组件; 27. 细胞器组件; 28. 大分子复合物; 29. 胞外区; 30. 细胞膜内控; 31. 细胞连接; 32. 共质体; 33. 超分子复合物; 34. 拟核; 35. 病毒粒子; 36. 病毒粒子组件; 37. 胞外区组件; 38. 催化活性; 39. 结合; 40. 转运活性; 41. 结构分子活性; 42. 核酸结合转录因子活性; 43. 电子载体活性; 44. 信号转导活性; 45. 分子功能调节因子; 46. 抗氧化活性; 47. 分子转导活性; 48. 营养库活性; 49. 转录因子活性，蛋白结合; 50. 蛋白标签; 51. 金属伴侣活性; 52. 翻译调控活性。

1. Metabolic process; 2. Cellular process; 3. Single-organism process; 4. Biological regulation; 5. Response to stimulus; 6. Localization; 7. Cellular component organization or biogenesis; 8. Developmental process; 9. Multicellular organismal process; 10. Signaling; 11. Reproduction; 12. Reproductive process; 13. Multi-organism process; 14. Growth; 15. Detoxification; 16. Immune system process; 17. Biological adhesion; 18. Biological phase; 19. Rhythmic process; 20. Cell killing; 21. Locomotion; 22. Cell; 23. Cell part; 24. Organelle; 25. Membrane; 26. Membrane part; 27. Organelle part; 28. Macromolecular; 29. Extracellular region; 30. Membrane-enclosed lumen; 31. Cell junction; 32. Symplast; 33. Supramolecular complex; 34. Nucleoid; 35. Virion; 36. Virion part; 37. Extracellular region part; 38. Catalytic activity; 39. Binding; 40. Transporter activity; 41. Structural molecule activity; 42. Nucleic acid binding transcription factor activity; 43. Electron carrier activity; 44. Signal transducer activity; 45. Molecular function regulator; 46. Antioxidant activity; 47. Molecular transducer activity; 48. Nutrient reservoir activity; 49. Transcription factor activity，protein binding; 50. Protein tag; 51. Metallochaperone activity; 52. Translation regulator activity.

1. 有机含硒化合物代谢; 2. 类胡萝卜素生物合成; 3. 类黄酮生物合成; 4. 糖尿病并发症中的年龄信号通路; 5. 淀粉和蔗糖代谢; 6. 碳代谢; 7. 叶酸碳池; 8. 半乳糖代谢; 9. 角质、亚氨酸和蜡的生物合成; 10. 单萜类生物合成; 11. 甘氨酸、丝氨酸和苏氨酸代谢; 12. 光合生物固碳作用; 13. 乙醛酸和二羧酸代谢; 14. 氮代谢; 15. 氰基氨基酸代谢; 16. 植物激素信号转导; 17. 卟啉和叶绿素代谢; 18. 苯丙素的生物合成; 19. 光合作用-天线蛋白; 20. 光合作用; 21. 抗坏血酸和醛酸代谢; 22. 硫代谢; 23. 双醌和其它萜烯醌生物合成; 24. 谷胱甘肽代谢; 25. 萜类骨架生物合成; 26. 磷酸戊糖途径; 27. 脂肪酸降解; 28. 糖酵解/糖新生; 29. 亚麻酸代谢; 30. 二萜生物合成; 31. 倍半萜和三萜生物合成。

1. Selenocompound metabolism; 2. Carotenoid biosynthesis; 3. Flavonoid biosynthesis; 4. AGE-RAGE signalling pathway in diabetic complications; 5. Starch and sucrose metabolism; 6. Carbon metabolism; 7. One carbon pool by folate; 8. Galactose metabolism; 9. Cutin，suberine and wax biosynthesis; 10. Monoterpenoid biosynthesis; 11. Glycine，serine and threonine metabolism; 12. Carbon fixation in photosynthetic organism; 13. Glyoxylate and dicarboxylate metabolism; 14. Nitrogen metabolism; 15. Cyanoamino acid metabolism; 16. Plant hormone signal transduction; 17. Porphyrin and chlorophyll metabolism; 18. Phenylpropanoid biosynthesis; 19. Photosynthesis-antenna proteins; 20. Photosynthesis; 21. Ascorbate and aldarate metabolism; 22. Sulfur metabolism; 23. Biquinone and other terpenoid-quinone biosynthesis; 24. Glutathione metabolism; 25. Terpenoid backbone biosynthesis; 26. Pentose phosphate pathway; 27. Fatty acid degradation; 28. Glycolysis/Gluconeogenesis; 29. Alpha-Linolenic acid metabolism; 30. Diterpenoid biosynthesis; 31. Sesquiterpenoid and triterpenoid biosynthesis.

茶文庫构建成功且测序质量良好，可进行下一步分析。SNP位点数目为489 642～638 706，转换类型比例为65.46%～66.33%，颠换类型为33.67%～34.54%，杂合型SNP位点比例为34.70%～51.38%。对261个基因进行了结构优化，发现了13 531个新基因，其中10 244个新基因被注释了，注释到COG数据库中有2 363个，注释到GO的有6 218个，注释到KEGG的有3 948个，注释到KOG的有5 833个，注释到Pfam的有5 383个，注释到Swiss-Prot的有6 900个，注释到eggNOG的有9 243个，注释到NR的有10 107个。

福鼎大白种茶根茎叶中有12 595个基因参与代谢过程，12 507个基因参与催化活性，8 176个基因参与细胞组成成分，是基因数最多的三个生物过程。陈琳波等（2015）对“紫鹃”茶树的转录组分析中，Unigene的GO分析表明，代谢过程、细胞、细胞部分、催化活性等富集程度较高，与本研究结果相一致。三个组织中有 1 615个基因参与发育过程（278，405，191），即叶和茎参与此过程的有278个差异基因，叶和根中有405个差异基因，茎和根中有191个。根茎叶有1 401个基因参与运输活动（197，356，231），288个基因参与生长过程（47，65，20），262个基因参与信号转导活动（44，70，30），200个基因参与免疫反应过程（35，56，36）。

都匀地区茶独特的风味是由叶片中的茶氨酸、萜烯类和类黄酮的代谢产物赋予的。Wu et al.（2013）从茶树叶片的转录组研究中得出大多数基因参与编码合成黄酮类、咖啡碱和茶氨酸等次生代谢途径的重要酶。该文通过KEGG功能显著性富集分析，发现叶和茎、叶和根、茎和根差异基因参与的部分极其显著的通路有类胡萝卜素生物合成，泛醌和其他萜烯类醌的生物合成，单萜生物合成，甘氨酸、丝氨酸和苏氨酸代谢，乙醛酸和二羧酸代谢，植物激素信号转导，卟啉和叶绿素代谢，苯丙素的生物合成，光合作用天线蛋白和光合作用，类黄酮生物合成，萜类化合物的生物合成，亚麻酸代谢等显著富集。叶和茎的差异表达基因在KEGG通路中，有90个基因参与碳代谢（8.65%，21，69），占所有差异基因的8.65%，其中有21个差异基因在茎中表达上调，69个差异基因在叶中表达上调。79个差异基因参与植物激素信号转导（7.6%，55，24），75个差异基因参与苯丙素的生物合成 （7.21%，49，26），70个差异基因参与淀粉和蔗糖代谢（6.73%，46，24，），64个差异基因参与氨基酸生物合成（6.15%，22，42），43个差异基因参与植物病原体互作通路（4.13%，30，13），19个差异基因参与类黄酮的生物合成（1.83%，13，6），20个差异基因参与萜类化合物生物合成（1.92%，4，16），24个差异基因参与卟啉和叶绿素代谢（2.31%，1，23），44个差异基因参与光合作用（4.23%，0，44）。以上研究表明都匀地区福鼎大白种茶独特的香味和各种健康功能主要归因于萜类物质合成、光合作用、 类黄酮的生物合成等生物过程，类黄酮主要赋予茶的苦味，二萜类化合物的水解产物具有令人愉悦的香味和特色风味。为进一步研究都匀地区福鼎大白种茶生长发育和遗传育种提供了分子基础，并为揭示与类黄酮相关代谢通路有关基因提供了重要理论依据。