茶树新梢中香叶醇樱草糖苷含量的全基因组关联分析

王让剑杨 军张力岚高香凤

1福建省农业科学院茶叶研究所, 福建福州 355012;2国家茶树改良中心福建分中心, 福建福州 355012

香叶醇是茶树中含有的一种重要的挥发性单萜醇, 在茶树与环境互作中起重要作用, 如抗菌[1]、驱虫[2]、吸引昆虫传粉[3]、提高相邻茶树对冷胁迫的抗性[4]等。香叶醇具有温和微甜的玫瑰花香, 也是茶叶中关键呈香成分之一, 其在茶鲜叶中主要以香叶醇樱草糖苷形式存在[5]。茶鲜叶中的糖苷通常位于液泡中, 糖苷酶存在于细胞壁和细胞间的空腔区域[6],正常细胞中二者相互隔离, 在茶树遭受环境胁迫或茶叶付制时, 茶鲜叶细胞遭受破坏, 二者发生接触,液泡内积累的糖苷会在糖苷酶作用下水解, 释放出香叶醇等挥发性化合物, 从而使茶树提升自身抵抗环境胁迫能力以及使成茶表现出愉悦的花香。

香叶醇等挥发性化合物与糖结合形成糖苷后,水溶性增大, 稳定性增强[7]。糖苷形式有利于挥发性物质在植物体内储存运输、信号转导[8], 还能以低毒或脱毒的形式减少或避免其对植物自身造成毒害[9]。茶叶中糖苷分离提取与定性定量的报道较多, 热水[10]、乙醇[11]、甲醇[12]浸提等方法均能分离提取茶叶糖苷, 水解[13]、衍生[14]、液相色谱质谱联用[15]等方法均可对茶叶糖苷定性定量。嫩叶和成熟叶中均有萜烯醇糖苷分布, 嫩叶含量高于成熟叶, 春季含量相对较高[16]。红茶加工过程中, 揉捻后萜烯醇糖苷含量明显降低, 发酵末期樱草糖苷基本消失, 但葡萄糖苷变化不大, 说明樱草糖苷是红茶香气重要来源物质[17]。研究表明, 茶树体内樱草糖苷可以在糖基转移酶的作用下通过葡萄糖基化和木糖基化顺序合成[18]。糖苷形成还可能与糖苷酶相关, 糖苷酶除了能水解糖苷键, 还能发生转糖基反应生成键合态糖苷[19]。

全基因组关联分析(genome-wide association study, GWAS)是一种重要的遗传解析方法, 已广泛应用于水稻[20]、小麦[21]、玉米[22]、大豆[23]等作物重要农艺性状候选基因的鉴定研究之中, 茶树中基于该方法对新梢萌发期[24]、风味代谢物[25]等重要性状亦进行了研究。虽然香叶醇樱草糖苷等萜醇类糖苷在茶树中的生物合成、茶鲜叶中的分布以及茶叶付制过程中的变化已见报道, 但鲜见基于群体在全基因组水平上发掘与香叶醇樱草糖苷含量相关的标记位点与候选基因的报道, 鲜叶中香叶醇樱草糖苷含量的遗传调控机制尚不清楚。本研究以169个茶树种质构成的自然群体为研究对象, 利用SLAF-seq技术(specific locus amplified fragment sequencing,SLAF-seq)开发覆盖茶树全基因组的SNP, 然后在3个环境下对茶树新梢香叶醇樱草糖苷含量性状进行GWAS分析, 在全基因组水平上发掘与香叶醇樱草糖苷含量性状显著关联的SNP位点与候选基因, 以期为香叶醇樱草糖苷含量的遗传调控机制研究及茶树分子标记辅助育种提供参考。

1 材料与方法

1.1 供试材料

茶树种质FT516为黄观音(♀)与铁观音(♂)人工杂交创新种质, 该种质制乌龙茶香气馥郁幽雅, 品质优异, 保存于福建省乌龙茶种质资源圃内(27°13′40″N,119°32′11″E)。供试材料为FT516的自然杂交后代,含169个茶树种质。2014年11月份收集FT516母树上成熟的自然杂交种子, 温室条件下沙培至发根状态。2015年3月种植到田间, 单株种植, 株距0.5 m,行距1.5 m, 田间管理措施相对一致。2018、2019和2020连续3年分别采摘各参试种质春季新梢无损伤无病虫害的单芽, 用于香叶醇樱草糖苷含量测定、SLAF-seq, 其中香叶醇樱草糖苷极端高含量种质409、极端低含量种质1605单芽用于转录组测序、基因组DNA重测序。另取相同栽培环境下的已应用品种悦茗香、福云6号单芽, 进行转录组测序进一步验证候选基因, 其中悦茗香为香叶醇樱草糖苷高含量品种, 福云6号为香叶醇樱草糖苷低含量品种[26]。

1.2 香叶醇樱草糖苷含量表型性状鉴定

按照文献方法测定参试茶树种质新梢单芽香叶醇樱草糖苷含量[26]。使用SAS8.0进行表型数据统计分析, 利用PROC VARCOMP程序进行方差分量评估, 并计算香叶醇樱草糖苷含量性状广义遗传力(h2),计算公式:h2=σg2/(σg2+σe2), 其中,σg2为遗传方差,σe2为环境方差。

1.3 群体SNP分子标记开发

利用CTAB法提取茶树基因组DNA, 选用HaeIII酶切茶树基因组DNA, 酶切片段长度在314～364 bp的序列定义为SLAF标签, 然后对检测合格的各样品基因组DNA分别进行酶切。对得到的酶切片段进行3′端加A处理、连接Dual-index测序接头、PCR扩增、纯化、混样、切胶选取目的片段, 文库质检合格后基于百迈客Illumina HiSeq平台进行测序。利用BWA[27]将测序reads比对到参考基因组上(参考基因组地址http://tpdb.shengxin.ren/), 使用GATK[28]和samtools[29]开发SNP, 以2种方法得到的SNP交集作为可靠SNP数据集。按照完整度大于0.8, 次要基因型频率大于0.05标准, 对SNP数据集进行过滤,使用过滤的SNP进行后续遗传分析。

1.4 群体遗传结构与连锁不平衡分析

使用admixture软件[30]进行群体遗传结构分析,分别假设样品的分群数(K值)为1～10, 并对聚类结果进行交叉验证, 根据交叉验证错误率的谷值确定最优分群数。使用EIGENSOFT软件[31]进行PCA主成分分析(principal components analysis, PCA)。使用plink[32]软件计算两两SNP间的连锁不平衡强度(R2),将SNP间距与R2进行拟合, 绘制R2随距离变化的LD衰减图(linkage disequilibrium, LD)。

1.5 全基因组关联分析

基于Tassle软件Glm模型进行GWAS全基因组关联分析, 使用admixture软件的Q矩阵结果作为固定效应[33]。显著性阈值P设置为:P＜1.48×10-6[P=(10n)-1,n为使用标记的总数]。根据GWAS结果获得显著关联的SNP位点, 基于群体连锁不平衡分析LD区间和茶树参考基因组注释信息, 获得与香叶醇樱草糖苷含量相关基因。

1.6 极端表型材料RNA-Seq与候选基因筛选

分别取极端高香叶醇樱草糖苷含量种质409和极端低香叶醇樱草糖苷含量种质1605春梢单芽, 每个样品取3个生物学重复, 在北京百迈克生物科技有限公司有参全长转录组平台完成转录组测序, 流程包括样品质量检测、文库构建、文库质量检测和文库测序(参考基因组http://tpdb.shengxin.ren/)。样本之间差异表达基因GO和KEGG富集分析在百迈客云平台完成(https://international.biocloud.net/)。另取栽培品种悦茗香、福云6号春梢单芽样品同时进行转录组测序, 进一步筛选验证候选基因。筛选差异表达基因条件为: 采用CPM (counts per million)作为衡量基因表达水平的指标(CPM=比对到某一转录本上的reads数目×106/比对到参考转录组的reads数目)[34]。将GWAS分析得到的与香叶醇樱草糖苷含量相关的基因在2个极端种质以及栽培品种悦茗香和福云6号差异表达基因中搜索, 重叠基因视为候选基因。

1.7 候选基因编码区序列差异及上游作用元件分析

在百迈客全基因组重测序分析平台完成409和1605基因组重测序, 流程包括样品质量检测、文库构建、文库质量检测和文库测序, 然后将重测序获得的无冗余reads重新定位到参考基因组上, 进行编码区SNP非同义突变分析(参考基因组http://tpdb.shengxin.ren/)。将候选基因起始密码子上游2 kb的序列作为启动子区域, 提交到在线工具Plant CARE(http://bioinformatics.psb.ugent.be/webtools/plantcare/html/), 默认参数, 进行顺式作用元件的预测。

2 结果与分析

2.1 香叶醇樱草糖苷含量表型变异分析

对参试种质香叶醇樱草糖苷含量性状进行鉴定发现, 各环境下群体香叶醇樱草糖苷含量变异系数在77.6%～81.8%之间, 说明该表型性状变异较丰富。正态分布检验结果表明(Ryan-Joiner检验), 各环境下群体香叶醇樱草糖苷含量性状r值在0.855～0.917之间, 接近正态分布, 说明该性状是由多基因控制的数量性状(表1)。方差分析结果表明, 香叶醇樱草糖苷含量性状在基因、环境间均呈现显著差异,依据方差组分估算广义遗传力(h2)为62.6%, 说明香叶醇樱草糖苷含量性状变异主要受遗传影响(表2)。

表1 香叶醇樱草糖苷含量表型变异Table 1 Phenotypic variation of geranyl β-primeveroside content

表2 香叶醇樱草糖苷含量方差分析Table 2 Variance analysis of geranyl β-primeveroside content

2.2 群体SNP开发

基于SLAF-seq简化基因组测序技术, 共开发群体SNP标记10,802,426个。根据SNP变异位点在基因组发生区域进一步细分类型, 大部分检测到的SNP变异发生在基因间区、内含子区域、基因上下游区域, 基因编码区的同义突变和非同义突变SNP数目较少(表3)。

表3 群体SNP标记统计Table 3 Statistics of SNP marker types in the population

2.3 群体遗传结构与连锁不平衡分析

经质量控制, 筛选出675,245个SNP用于后续遗传分析, 每个标记间的平均距离为4.47 kb (3.02 Gb/675,245)。利用admixture软件分析, 当K=2时(图1-A), ΔK值最大, 参试群体被划分为2个亚群(图1-B), 利用 EIGENSOFT软件进行主成分分析(principal components analysis, PCA), 参试群体也被划分为2个亚群(图1-C), 2种方法划分结果一致。利用plink软件绘制R2随距离变化图发现,R2从0.2下降到0.1时, LD衰减距离约为100 kb (图1-D)。标记间的平均距离远小于LD衰减距离, 说明标记密度满足GWAS分析要求。

图1 遗传结构及连锁不平衡分析Fig. 1 Analysis of genetic structure and linkage disequilibriumA: 遗传结构分析; B: 主成分分析; C: 连锁不平衡分析。A: the genetic structure analysis; B: the principal component analysis; C: the linkage disequilibrium analysis.

2.4 全基因组关联分析

基于TASSEL软件GLM模型对香叶醇樱草糖苷含量性状进行GWAS分析, 3个环境下共检测到340个SNP标记与香叶醇樱草糖苷浓度显著相关,其中2个环境下均检测到的SNP标记共65个, 各标记贡献率在17.1%～23.9%之间, 分布在29个Scaffold上, 均位于基因间区、基因内含子以及基因下游区域(表4)。

表4 香叶醇樱草糖苷含量相关的SNP位点Table 4 SNP sites related to geranyl β-primeveroside content

2.5 香叶醇樱草糖苷含量候选基因预测

基于茶树参考基因组信息(http://tpdb.shengxin.ren/)和LD衰减距离(100 kb), 获得与香叶醇樱草糖苷含量相关65个SNP标记上、下游100 kb范围内的基因共88个, 涵盖了信号蛋白、激酶、磷酸酶、离子转运蛋白、转录因子、热激蛋白、激素相关蛋白、抗性蛋白、萜类代谢酶、糖基转移酶、糖苷酶类等(表5)。

表5 候选基因信息Table 5 Candidate gene information

(续表5)

(续表5)

2.6 极端表型材料RNA-Seq分析

利用RNA-Seq获得2个香叶醇樱草糖苷极端含量种质差异表达基因共4393个, 其中显著上调表达基因2049个(图2-A), 显著下调表达基因2344个(图2-B)。GO功能分析结果表明, 差异表达基因生物过程最多的2类均是细胞过程(cellular process)、代谢过程(metabolic process), 细胞结构最多的2类均是膜(membrane)、膜构件(membrane part), 分子功能最多的2类均是结合活性(binding)、催化活性(catalytic activity) (图3-A, B)。KEGG功能分析结果表明, 上调差异表达基因主要富集在油菜素甾醇、类黄酮、核糖体蛋白、氨基酸等代谢途径(图4-A), 下调差异表达基因主要富集在单萜类、倍半萜类、花生四烯酸类、玉米素类等代谢途径(图4-B)。

图2 差异表达基因Fig. 2 Differentially expressed genesA: 409 vs 1605显著上调表达; B: 409 vs 1605显著下调表达。A: 409 vs 1605 significantly up regulated; B: 409 vs 1605 significantly down regulated.

(图3)

图3 差异表达基因GO功能分析Fig. 3 GO functional analysis of differentially expressed genesA: 409 vs 1605显著上调表达; B: 409 vs 1605显著下调表达。A: 409 vs 1605 significantly up regulated; B: 409 vs 1605 significantly down regulated.

图4 差异表达基因KEGG功能分析Fig. 4 KEGG functional analysis of differentially expressed genesA: 409 vs 1605显著上调表达; B: 409 vs 1605显著下调表达。A: 409 vs 1605 significantly up regulated; B: 409 vs 1605 significantly down regulated.

2.7 香叶醇樱草糖苷含量候选基因筛选

将GWAS分析得到的显著关联SNP标记上、下游100 kb范围内的相关基因在2个香叶醇樱草糖苷极端含量种质差异表达基因中搜索, 得到重叠基因10个, 均表现为高含量种质中表达量高于低含量种质。利用栽培品种悦茗香、福云6号进一步验证10个重叠基因的表达水平发现, 10个重叠基因在香叶醇樱草糖苷高含量品种悦茗香中表达量亦高于香叶醇樱草糖苷低含量品种福云6号(表6), 故将10个重叠基因初步视为候选基因。

表6 香叶醇樱草糖苷含量候选基因预测Table 6 Prediction of candidate genes of geranyl β-primeveroside contents

2.8 候选基因编码区序列差异及上游作用元件分析

对2个极端香叶醇樱草糖苷极端含量种质进行基因组DNA重测序, 10个候选基因编码区序列在2个极端种质之间均存在SNP非同义突变, 造成了不同种质之间编码蛋白氨基酸序列的差异(表7)。利用Plant CARE在线工具对所有10个候选基因上游2 kb的区域进行了预测(表8)。在所有类型的元件中, 与激素响应相关的元件及与环境胁迫相关的元件是代表植物核心生理过程的调控元件。候选基因中具有与激素响应相关的顺式作用元件的基因占8个(除TEA023801.1、TEA009864.1外), 其中与赤霉素、脱落酸和茉莉酸甲酯有关的顺式作用元件含量较丰富。候选基因中具有与环境胁迫相关的顺式作用元件的基因占9个(除TEA009864.1外), 其中与防御和应激反应性、干旱诱导、无氧诱导有关的元件含量较丰富。

表7 2个极端种质中差异基因编码区序列变异Table 7 Sequence variation of differential gene coding region in two extreme germplasms

(续表7)

表8 茶树香叶醇樱草糖苷候选基因上游与环境相关的顺式作用元件的预测Table 8 Prediction of environmental relatedcis-elements in the upstream regions of candidate genes for tea geraniol primrose glycoside

(续表8)

3 讨论

分析了参试茶树种质在3个环境下新梢中香叶醇樱草糖苷含量的遗传差异, 其广义遗传力为62.6%, 说明香叶醇樱草糖苷含量变异主要受遗传控制, 其候选基因编码区序列均存在SNP非同义突变, 预示可以通过遗传改良手段定向调控茶树新梢中香叶醇樱草糖苷含量。候选基因中多数基因上游2 kb区域存在不同数量的与激素信号转导和环境胁迫相关的顺式作用元件, 预示激素水平和环境变化会影响茶树香叶醇樱草糖苷的积累。候选基因涉及到信号转导、基因表达、代谢变化等多个过程(图5),说明茶树新梢香叶醇樱草糖苷积累的遗传调控机制是一个复杂的体系, 由基因和环境共同决定, 它们之间的相互作用仍有待深入研究。

图5 候选基因示意图Fig. 5 Schematic diagram of candidate genes括号中的数字表示基因个数。信号转导: TEA023908.1, TEA006075.1, TEA029081.1, TEA027247.1; 基因转录: TEA023801.1; 蛋白翻译:TEA023811.1, TEA009864.1; 抗病蛋白: TEA006076.1; 糖代谢: TEA001888.1, TEA018925.1。Numbers in brackets indicate the number of genes. Signal transduction: TEA023908.1, TEA006075.1, TEA029081.1, TEA027247.1; Gene transcription: TEA023801.1; Protein translation: TEA023811.1, TEA009864.1; Resistant protein: TEA006076.1; Glucose metabolism:TEA001888.1, TEA018925.1.

富含亮氨酸重复类受体蛋白激酶(LRR-RLK)包含1个胞外的富含亮氨酸(LRR)结构域、1个跨膜结构域和1个胞内丝氨酸/苏氨酸(Ser/Thr)激酶结构域[35]。LRR-RLK定位在细胞质膜上, 能够感知外界环境信号, 然后通过磷酸化将信号传递到胞内相应部位, 可作为信号感知元件参与植物体应答生物和非生物胁迫过程[36]。不同的LRR-RLK中LRR数目和排布有较大的区别, LRR的多样性使LRR-RLK能够特异的识别不同的信号分子, 从而准确传递不同的胞外信号[37]。水稻接种稻瘟病菌后,OsBRR1在叶片中的表达受强烈诱导, 同时过量表达OsBRR1转基因植株可以提高稻瘟病抗性[38]。水稻中的一个LRR-RLK蛋白激酶还能正调控与生物胁迫相关的MAPK及WRKY表达水平, 从而调控植物对昆虫取食的防御机制[39]。烟草中LRR-RLK基因参与了盐胁迫反应的调控过程[40]。辣椒中LRR类受体蛋白激酶基因CaLRR1的表达受病原菌侵染的诱导[41]。丝氨酸/苏氨酸蛋白激酶在接受到生物和非生物胁迫信号后迅速在其丝氨酸、苏氨酸部位发生磷酸化而被激活, 并通过对下游分子的磷酸化作用, 最终将外界信号传递给细胞核, 调节特异基因的表达[42]。本研究发现2个富含亮氨酸重复类受体蛋白激酶基因(TEA023908.1、TEA006075.1)和2个丝氨酸/苏氨酸蛋白激酶(TEA029081.1、TEA027247.1)与香叶醇樱草糖苷的含量有关, 推测这些蛋白激酶作为信号转导分子会将环境胁迫信号传递到代谢网络下游, 导致茶树香叶醇樱草糖苷发生代谢变化, 从而引发香叶醇在茶树与环境互作中发挥重要作用。

茶树体内樱草糖苷可在糖基转移酶作用下通过葡萄糖基化和木糖基化顺序合成, 虽然鼠李糖、半乳糖、木糖、葡萄糖醛酸均可作为糖基转移酶的活性糖供体, 但一般认为葡萄糖是最常见的糖基供体[18]。β-1,3-葡聚糖酶能够水解β-1,3-葡聚糖成为葡萄糖和小分子寡糖, 在正常情况下β-1,3-葡聚糖酶表达水平较低, 当植物遭受生物因子以及非生物因子诱导后, 其表达活性会显著增强[43]。研究表明, 棉花β-1,3-葡聚糖酶基因参与了棉花对黄萎病的防卫反应[44]。烟草β-1,3-葡聚糖酶基因转化番茄可使番茄发病率降低36%～58%[45]。本研究鉴定到2个β-1,3-葡聚糖酶基因(TEA001888.1、TEA018925.1)与香叶醇樱草糖苷累积相关, 推测β-1,3-葡聚糖酶通过水解作用产生葡萄糖为香叶醇樱草糖苷生物合成提供糖配基并参与调节香叶醇的累积和释放, 从而抵御生物和非生物胁迫。

植物在生长发育过程中为了应对复杂的环境变化, 进化形成了2个不同层次的防御反应, 初级免疫反应(PTI)和次级免疫反应(ETI), 其中ETI比PTI更为强烈和迅速, 可以通过局部程序性细胞死亡引起植物超敏反应, 能更为有效地遏制病原菌的进一步扩散[46]。RPM1是一种植物抗病蛋白, 可以触发ETI使植物获得抗病性[47]。本研究发现1个抗病蛋白RPM1基因(TEA006076.1)与香叶醇樱草糖苷含量相关, 推断香叶醇樱草糖苷以及香叶醇在与环境互作过程中会协同产生抗病蛋白RPM1, 从而提高茶树抗病性。

植物萜类化合物的合成途径分为甲羟戊酸(MVA)和2-C-甲基-D-赤藓糖醇-4-磷酸(MEP)2种途径, 前者位于细胞质中, 后者位于质体中, 二者虽有所不同, 但共同点都以异戊烯基焦磷酸(IPP)为主要中间产物。IPP与其异构体二甲丙烯焦磷酸(DMAPP)在异戊烯基二磷酸合酶(IPPS)催化下生成香叶基焦磷酸(GPP), 后者在萜类合成酶(TPS)的作用下生成香叶醇[48]。本研究在Scaffold349上280 bp范围内鉴定到6个SNP, 分别位于3,542,803 bp、3,543,009 bp、3,543,014 bp、3,543,021 bp、3,543,061 bp、3,543,082 bp位置, 在该区段10.9 kb和85.8 kb附近区域分别存在 1个核糖体蛋白基因(TEA023811.1)和1个RNA聚合酶II转录亚基基因(TEA023801.1), 转录组测序表明这2个基因在香叶醇樱草糖苷极端含量种质中均表现显著差异表达,说明他们可能与香叶醇樱草糖苷含量相关。有趣的是在该区段8.6 kb区域附近还存在1个异戊烯基二磷酸合酶基因(TEA023826.1), 推测核糖体蛋白基因(TEA023811.1)和RNA聚合酶II转录亚基基因(TEA023801.1)参与异戊烯基二磷酸合酶基因(TEA023826.1)的转录和翻译, 进而导致香叶醇及其樱草糖苷代谢和累积变化。

萜类化合物与单糖(或寡糖)一起作为底物, 可在糖基转移酶的催化下生成不易挥发的糖苷化合物。虽然本研究筛选的候选基因中未包括糖基转移酶基因, 但是在Scaffold3611的1,622,376 bp位置鉴定到1个SNP, 其附近14.9 kb区域存在一个糖基转移酶基因(TEA013410.1)和一个WRKY转录因子基因(TEA013393.1), 他们是否参与香叶醇樱草糖苷的生物合成, 仍有待证实。

4 结论

对茶树新梢中香叶醇樱草糖苷含量进行了GWAS分析, 初步鉴定到10个与香叶醇樱草糖苷含量相关的候选基因, 其编码区序列在香叶醇樱草糖苷极端含量种质之间均存在非同义SNP突变, 并且多数基因上游2 kb区域存在不同数量的与环境胁迫和激素相关的顺式作用元件, 说明香叶醇樱草糖苷的含量变异与遗传及环境胁迫相关。该结果对茶树香叶醇樱草糖苷含量的遗传调控机制研究及育种实践具有较好的理论和应用参考价值。