茶树叶型性状的全基因组关联分析

王小萍,王 云,熊元元,刘 晓,李明红,张 厅,赖 谦,李春华

(1.四川省农业科学院茶叶研究所,成都 610066;2.泸州市经济作物站,四川 泸州 646000;3.古蔺县农业农村局,四川 泸州 646500)

【研究意义】茶树叶型与茶树抗性强弱、产量和品质密切相关。叶长、叶宽、长宽比、叶身、叶色等是茶树重要的叶型性状,易受环境和遗传等多因素的共同影响,属于典型的数量性状[1]。研究茶树叶型性状的遗传特性,选育具有合理叶型性状的新品种,对促进茶树高效优质育种具有积极意义。【前人研究进展】关联分析技术可直接对基因型变异和表型变异进行分析,从而发掘特定条件下与性状相关联的候选基因[2]。随着基因组测序技术的快速发展,全基因组关联分析(GWAS)已在玉米[3-4]、水稻[5-6]、小麦[7]、油菜[8]、甜瓜等[9]的相关性状研究和基因发掘中得到广泛应用,而在茶树上的应用起步晚且较少。 Lu等[10]利用全基因组关联分析对来自贵州、云南的8个古茶树群体进行遗传结构及树型与叶型相关性状进行分析,挖掘出6个候选基因,与4个叶型性状和2个树型性状显著相关。Hazra等[11]利用全基因组关联技术对23份大吉岭茶树种质的农艺与品质性状进行研究,最终检测到57个SNP显著位点,多与茶树品质性状相关联。Huang等[12]对“龙井43”“白鸡冠”的198份F1代植株开展游离氨基酸含量的关联分析,挖掘出4个关联候选基因。【本研究切入点】虽然前人基于茶树叶型性状研究发掘出部分候选基因,但茶树叶型受多基因控制与影响,有必要发掘更多的关联位点和候选基因辅助分子标记育种。【拟解决的关键问题】对收集于四川本地的91份茶树自然群体开展简化基因组测序,基于测序获得的SNP开展叶型性状的全基因组关联分析,筛选与叶型性状关联的SNP位点,并挖掘相关的候选基因,为今后茶树分子标记辅助育种提供基因资源和参考依据。

1 材料与方法

1.1 试验材料

91份茶树资源(表1)包括70份野生乔木型茶树资源[古蔺52份(以G1,G2......表示)、叙永5份(以X1,X2......表示)、合江1份(以H1表示)、雷波4份(以LB1,LB2......表示)、崇州6份(以CP1,CP2......表示)、云南大叶种2份(以YD1、YD2表示)],8份北川野生茶树(以BC1,BC2......表示)和13份四川中小叶群体种材料(以CCXX表示,CC表示川茶首字母,XX表示具体来源地首字母,如CCJA表示源自江安县的四川中小叶群体种,GYGC表示“高阳贡茶”)。

表1 91份茶树资源

续表1 Continued table 1

1.2 茶树资源表型性状统计

根据《茶树种质资源描述规范和数据标准》于2020年10—11月调查并测量茶树的叶长(YC)、叶宽(YK)、长宽比(CKB)、叶身(YS)、叶色(YSZ)、叶尖(YD)、叶脉对数(YM)和叶缘(YY)共8个性状。其中叶身(YS)、叶色(YSZ)、叶尖(YD)和叶缘(YY)4个描述型性状参照《茶树种质资源描述规范和数据标准》转化为数值数据,使用Excel整理调查统计的数据,并利用R软件(V3.5.0)对表型数据进行统计分析,计算均值、标准差、标准误、偏度、峰度等。

1.3 茶树资源群体结构与亲缘关系分析

于2021年5—6月取参试资源夏梢幼的嫩芽叶,采用改良CTAB法提取基因组DNA,质检合格的基因组DNA由广东基迪奥生物科技公司进行30×的全基因组简化测序(GBS),将测序数据进行过滤和比对,获得高质量的SNP标记。基于高质量的SNP标记,利用Admixture软件(V1.3)[13]进行供试材料群体结构推断,根据交叉验证错误率(CV)确定群体最优分群数(交叉验证错误最小值对应最优的分群数),即K值;同时利用gcta软件(v1.92.2)[14]进行亲缘关系分析,获得供试材料两两间亲缘关系矩阵。

1.4 茶树资源全基因组关联分析

根据最优K值计算各材料相应Q值,作为固定效应协变量矩阵(Q矩阵),以亲缘关系矩阵作为随机效应协变量矩阵(K矩阵),应用Tassel软件(v5.2.54)[15]的广义线性模型(Q模型),基于SNP标记对供试91份茶树资源叶型8个性状的表型数据进行全基因组关联分析(GWAS)。取P<0.01时与性状显著关联的SNP位点,并在显著位点上下游50 kbp区域内进行候选基因挖掘[7]。根据候选基因GO与KEGG 注释,进一步推断预测与叶型性状关联的功能基因。

2 结果与分析

2.1 茶树资源表型性状分析

对91份茶树资源叶型8个性状进行分析(表2)显示,8个叶型性状均表现较大变异,其中叶尖的变异系数最大(47.18%),长宽比变异系数最小(11.34%)。 4个数值型性状中叶长、叶宽、长宽比和叶脉对数均值分别是12.44 cm、5.26 cm、2.38和8.49对,波动范围较大。4个描述型性状中叶身、叶色、叶尖和叶缘的均值分别是1.51、3.14、1.49和2.18,且4个性状在所有赋值等级均有分布。8个性状中,除性状长宽比和叶色的峰度数值较大外,其余性状偏度、峰度接近1,基本符合正态分布。将各性状进行相关性分析(表3),叶长与叶宽相关性最大(r=0.86),其次是叶长与叶脉对数(r=0.57),分别呈显著正相关;叶长与叶尖(r=-0.59)、叶宽与叶尖(r=-0.52)分别呈显著负相关;除叶缘、叶色2个性状与其他性状相关系数较低以外,其余两两性状间均存在一定的相关性,表明这8个叶型性状间存在协同变化的关系。

2.2 茶树资源群体结构与亲缘关系分析

对91份茶树资源基因组DNA进行分析,构建GBS文库及酶切测序,共获得97.79 G数据,过滤后获得95.69 G数据,734 642 116条高质量reads,平均每份茶树资源获得1.05 G数据,8072 990条高质量reads。91份茶树资源Q20≥97.27%,Q30≥91.84%,表明测序质量较好,GC含量在40.02%～41.76%,分布正常。将高质量reads比对到茶树参考基因组上,参试茶树共获得5569 787个变异,其中5374 140个SNPs,195 647个Indels。

基于筛选获得863 468个高质量的SNP标记,利用Admixture软件(V1.3)对参试茶树资源的遗传结构进行分析,当K=2时,CV值最小(图1-A),参试茶树资源被分为2个亚群(图1-B)。其中亚群Ⅰ包含22份茶树资源,分别为8份北川茶树、11份四川中小叶群体种、2份云南茶树和1份古蔺茶树;亚群Ⅱ包含69份茶树,以古蔺、叙永、崇州、合江、雷波的野生茶树为主,同时还包含2份四川中小叶群体种茶树。

表2 91份茶树资源叶型相关性状基本统计

表3 叶型相关性状的相关性分析

图1 91份茶树资源的群体结构分析Fig.1 Population structure analysis of 91 tea resources

进一步利用gcta软件(v1.92.2)对茶树资源进行亲缘关系分析,获得91份茶树资源两两间亲缘关系。由图2可知, 38份茶树资源(占比41.75%)两两材料间亲缘关系小于0, 48份茶树资源(占比52.75%)两两材料间亲缘关系介于0～0.5,仅5份资源(占比5.49%)的两两材料间亲缘关系介于0.5～1.5。

图2 91份茶树资源的亲缘关系Fig.2 Kinship of 91 tea resources

2.3 茶树资源叶型相关性状的关联分析

基于863 468个SNP标记对91份茶树资源叶型8个性状的表型数据进行关联分析,由表4可知,在P<0.01时,8个性状中的叶长(YC)、长宽比(CKB)、叶色泽(YSZ)、叶尖(YD)、叶缘(YY)、叶脉(YM) 6个性状共检测到38个显著SNP位点,叶宽和叶身2个性状未检测到显著SNP位点。38个显著位点中：叶长(YC)关联到1个显著位点,长宽比(CKB)关联到的显著位点最多,为24个,叶色泽(YSZ)关联到6个显著位点,叶尖(YD)关联到2个显著位点,叶脉对数关联到1个显著位点,叶缘则关联到4个显著位点。这些SNP显著位点对关联性状的解释率在31.37%～60.70%,从SNP显著位点的位置信息来看,其中34个(占比89.5%)关联位点位于基因间区,3个(占比7.9%)关联位点位于内含子区域,1个(占比2.6%)关联位点位于转录终止下游1 kb区域内。

表4 基于SNP标记的关联分析检测到的叶型相关性状显著位点

续表4 Continued table 4

2.4 茶树资源候选基因的注释与预测

基于参考基因组的序列数据,对获得的SNP显著位点上下游50 kbp区域内进行扫描,共筛选出潜在的候选基因8个,其中CSA011953(Sec11a真核信号肽酶s26b)、CSA027293(CAD3肉桂醇脱氢酶3)、CSA020561(GAE2udp-葡萄糖醛酸4-外酯酶3)、CSA008530(AIM1过氧化物酶体脂肪酸-氧化多功能蛋白靶向1)、CSA008531(MFP过氧化物酶体脂肪酸-氧化多功能蛋白靶向1)、CSA008532(AIM1过氧化物酶体脂肪酸-氧化多功能蛋白靶向1)与性状长宽比(CKB)相关联,CSA008763(CYP72A154细胞色素p450)与叶片色泽(YSZ)相关联,CSA029592(RABG3Fras相关蛋白rab7亚型x1)与叶脉对数(YM)相关联(图3)。对8个候选基因进行GO和KEGG注释(表5),初步推测基因CSA027293(CAD3)与长宽比(CKB)紧密关联、CSA008763(CYP72A154)与叶片色泽(YSZ)紧密关联。

图3 8个候选基因在性状关联SNP位点中的位置Fig.3 Location of eight candidate genes in related SNP markers

表5 8个候选基因的GO和KEGG注释信息

3 讨 论

本研究对茶树成熟叶片的叶长(YC)、叶宽(YK)、长宽比(CKB)、叶身(YS)、叶色(YSZ)、叶尖(YD)、叶脉对数(YM)和叶缘(YY)8个性状进行表型分析,多样性较高,变异系数介于11.34%～47.18%,与其他以野生自然群体茶树为研究对象的相关研究相比(性状相同),高于贵州姑菁野生茶树资源的表型多样性[16],略低于云南野生茶树资源[17-18]。

在遗传结构分析中,91份茶树资源分为2个亚群,源自2个遗传血统。从茶树资源来源看,亚群Ⅰ中四川中小叶群体种取自德昌、广安、平武、旺苍、宣汉、天全、江安等地,北川茶树资源取自山林间,零星分布,当地称之为“野生茶”,云南大叶种则源自德宏州芒市当地大叶种群体自然杂交后代,这3个自然群体遗传结构一致,具有共同血统来源,初步印证四川中小叶群体种部分源自云南茶果,但还需进一步扩大样本数量证实;亚群Ⅱ中,16份茶树资源(含14份古蔺大树茶与2份雷波大叶种)遗传结构一致且单一,有54份茶树资源则存在亚群Ⅰ和亚群Ⅱ 2个血统杂合,表明2个亚群间存在基因渗透交流的趋势,但两者仍然独立且亚群Ⅱ(野生型)仍在遗传结构中占据主导。这也间接印证了钟渭基[19]认为四川野生大茶树系本地原产的观点,也与周斌等[20]对雷波、云南茶树遗传结构研究中,认为四川与云南的野生茶树群落相互独立观点一致。

叶片是茶树重要的生长器官之一,同时也是茶树的直接利用部位。叶片形态不仅直接影响光合效率(产量),还与茶树抗性强弱密切相关。Lu等[10]对茶树叶型和树型共11个性状开展GWAS研究中,定位发现2个基因TEA01277与TEA02816与茶树叶色紧密关联,基因TEA01277属编码钙调素结合转录激活因子2,基因TEA02816则属cop1/spa泛素连接酶复合物。本研究中,筛选出基因CSA008763(CYP72A154)与叶色(YSZ)紧密关联。该基因属细胞色素P450家族CYP72A亚家族,CYP72A154被证实与甘草酸合成途径中的C-30氧化有关[21],并在野生型甘草中高度表达[22]。另外,基因CSA027293(CAD3)与长宽比紧密关联。CAD(肉桂醇脱氢酶)是木质素单体合成最后一步还原反应关键酶,目前CAD基因已在小麦[23]、茶树[24]、杨树[25]、陆地棉[26]、棉花[27]、黑莓[28]等物种中克隆鉴定得到。研究表明,CAD基因的表达与植株木质素累积或木质素单体结构组成等有关。如在小麦[23]、高粱[29]中,CAD4基因与抗倒伏能力有关,而茶树中基因CAD3在新梢第5叶中高度表达,与茶树鲜叶的木质化程度密切相关[24]。后续将对进一步验证这2个候选基因的准确性开展研究。

4 结 论

基于简化基因组重测序(GBS)对四川本地的91份茶树自然群体进行分析,利用测序获得的863 468个高质量的SNP标记对叶型8个性状(叶长、叶宽、长宽比、叶身、叶色、叶尖、叶脉对数和叶缘)开展全基因组关联分析,获得38个显著SNP关联位点,并初步推测基因CSA008763(CYP72A154)与叶色、基因CSA027293(CAD3)与长宽比紧密关联,这为今后茶树分子标记辅助育种提供基因资源和理论参考。