大叶杨干旱和水淹胁迫的转录组分析

张佳琪,张新新,冯佳珺,施婷婷,张钰麟,王 婧

(四川大学 生命科学学院 生物资源与生态环境教育部重点实验室,成都 610000)

大叶杨(Populuslasiocarpa)是杨柳科(Salicaceae)杨属(Populus)大叶杨组的落叶乔木,主要分布于中国湖北、四川、陕西、云南、贵州等地,生长在海拔1 300～1 500 m的山坡、灌丛等水资源较为丰富的区域[1],是中国特有的杨属物种。相较于其他广域分布的杨属植物如山杨、银白杨、青杨等[2],大叶杨的分布范围较为局限,这可能与其喜湿的生长习性有关。近年来,全球气候变暖导致极端气候事件频发,这可能导致干燥地区变得更加干燥,潮湿地区则变得更加潮湿[3]。IPCC AR6报告指出,在过去的50余年,全世界范围内强降水事件频发以及干旱事件增加[4],而由此引发的干旱和洪水等自然灾害已经严重阻碍了中国林业建设与发展和生态环境的改善。因此在地区生境条件日益恶劣的背景下,培育抗逆性较强的树种并制定相应的保护和育种措施是恢复与重建生态系统平衡和提高林业生产力的重要方式。

长期的洪涝会使得土壤含水量处于饱和状态,植物的有氧呼吸逐渐替换为无氧呼吸[5],其形态结构、生长生理和新陈代谢等逐步受到影响[6-8]。如枫杨在水淹12 h后,体内类黄酮合成、亚麻酸、茉莉酸等代谢途径明显活跃[9]。而缺水造成的干旱胁迫可能导致农作物和林木的减产,超过其他环境因子胁迫所造成减产的总和[10]。另外,缺水环境下也能够引起植物的叶子或幼果释放大量的乙烯,从而脱叶、落果[11]。植物也可通过增加体内的脱落酸含量以提高抗旱能力[12]。因此为了维持自身的生长发育,植物在分子水平上已进化出一系列响应非生物胁迫的机制。而外界的非生物胁迫也能影响植物体内功能基因和调控基因的表达[13-14],其中就包括调控基因表达的众多转录因子家族成员[15]。AP2/ERF转录因子家族在植物的生长发育、代谢以及生物和非生物胁迫的响应中发挥了关键作用[16]。根据AP2/ERF结构域的数量以及是否含有其他结构域将其分为5个亚家族[17]:AP2、ERF(响应乙烯)、DREB(与脱水响应元件结合)、RAV、Soloist。AP2亚家族含有2个串联重复的AP2/ERF结构域;ERF和DREB亚家族中只含有1个AP2结构域,其共同特点是基因结构中内含子数量很少[18],但是二者具有不同的作用,ERF亚家族主要与乙烯响应元件结合[19],而DREB亚家族,主要通过DRE/CRT顺式作用元件来调控基因的表达[20];Soloist亚家族也只含有1个AP2/ERF结构域,但是核苷酸序列和基因结构与ERF等亚家族明显不同,并且核苷酸序列在很多植物中高度保守[17];RAV亚家族除了含有1个AP2/ERF结构域外,还带有1个B3结构域。AP2亚家族主要在植物生长发育阶段起重要作用,例如开花[21]、种子生长与萌发[22]等;ERF主要作用于激素信号传导、非生物胁迫响应等[23];RAV主要参与叶片衰老、生物与非生物胁迫响应调控[24];Soloist在拟南芥中仅包含1个基因,调控水杨酸的合成以及防御细菌病原体的入侵[25]。

目前对于大叶杨的研究多集中于种质资源调查和保存[26],而对于大叶杨如何响应非生物胁迫以及其胁迫耐受能力研究鲜有报道。本研究通过对水淹和干旱胁迫后的大叶杨进行转录组测序,结合生物信息学分析鉴定到一系列响应胁迫的差异表达基因,并发现大叶杨中AP2/ERF转录因子家族成员在胁迫响应基因中显著富集,最后筛选出可能参与调控干旱或水淹胁迫耐受的候选基因。该研究结果将为今后深入探究极端水分环境下大叶杨的分子响应过程和适应性进化奠定基础,为大叶杨分子遗传育种提供理论参考。

1 材料和方法

1.1 试验材料

大叶杨种子采集于四川省凉山彝族自治州美姑县瓦候乡(28°43′N,103°8′E),铺种萌发后移至温室进行培养,温室生长条件为:光照16 h,(25±1) ℃,黑暗8 h,(20±1) ℃,空气湿度在70%左右。待植株高度生长到30～40 cm时,挑选一部分进行预实验以此确定水分胁迫时间(以出现表型变化但未枯死不影响测序为准),最终确定干旱胁迫时间为5 d,水淹胁迫为15 d。挑选出长势均匀良好的植株用于正式试验。

1.2 胁迫处理

正式试验将材料分为3组:水淹(waterlogging stress,WS)处理组(水没过盆土3 cm左右)、干旱(drought stress,DS)处理组(胁迫开始前保持正常浇水,胁迫开始后停止浇水,干旱胁迫结束时土壤含水量约在56%)以及对照(control group,CK)控制组(保持正常浇水,土壤含水量在70%～80%)。为确保胁迫结束的时间一致,于2021年11月18日早上9:00先开始水淹胁迫处理,在水淹胁迫的第10天早上9:00 DS组材料停止浇水,开始干旱胁迫处理,期间保持所有处理组及对照组材料在相同环境,水淹处理组的水每3 d更换1次,以免微生物滋生干扰结果。

在水淹和干旱胁迫第3天和胁迫结束的最后一天采集样品,每组各取3株叶片材料用于转录组学测序,叶片统一取第5～6片新鲜的成熟叶片,并用液氮速冻后暂存于-80 ℃冰箱。将试验样品材料送至深圳华大基因股份有限公司(https://www.bgi.com/),使用T7 PE150平台进行cDNA文库构建和RNA-Seq测序,每组处理的3株材料作为生物学重复。

1.3 生物信息学分析

原始测序数据(raw reads)先用FastQC[27]检验数据质量,Trimmomatic[28]过滤接头及低质量的reads。过滤后的reads用HISAT2[29]比对到组装好的大叶杨参考基因组上,再用StringTie[30]统计基因表达量,得到标准化后的TPM(transcripts per kilobase of exon model per million mapped reads)值。根据脚本prepDE.py(https://github.com/gpertea/stringtie/blob/master/prepDE.py)生成用于DESeq2[31]差异表达分析的矩阵(count)文件,将错误发现率FDR(false discovery rate)≤0.05,差异倍数|log2αFD|≥1的基因定义为差异表达基因(differentially expressed genes,DEGs)。topGO[32]包对DEGs进行GO富集分析。

1.4 大叶杨AP2/ERF转录因子家族分析

在Pfam蛋白质家族数据库(https://pfam-legacy.xfam.org/)下载AP2保守结构域的蛋白序列(序列号:PF00847)HMM文件,HMMER v3.1b2[33](http://hmmer.org/)鉴定大叶杨基因组数据库中的AP2/ERF候选基因。SMART(https://smart.embl.de/)和NCBI CDD(https://www.ncbi.nlm.nih.gov/cdd/)在线数据库对基因蛋白序列进行检测,保留有AP2保守结构域的序列用于后续分析。

在PlantTFDB植物转录因子数据库(http://planttfdb.gao-lab.org/)中下载拟南芥AP2/ERF基因家族成员的蛋白质序列。Muscle v3.8.31[34]对大叶杨和拟南芥的目的基因蛋白序列进行多序列比对,trimAL v.1.4.1[35]快速修剪和对齐,MEGA11[36]中采用 NJ(neighbor-joining)法构建系统发育树,Bootstrap设置为1 000,在FigTree v1.4.4(http://tree.bio.ed.ac.uk/software/figtree/)对树进行可视化。

1.5 qRT-PCR验证

提取RNA后反转录得到cDNA原液用于实时荧光定量PCR[37]实验。NCBI Primer(https://www.ncbi.nlm.nih.gov/tools/primer-blast/)对12个大叶杨AP2/ERF家族成员在线设计候选基因和内参基因的特异性引物(表1),UBQ-10为内参基因[38]。采用南京诺唯赞生物科技股份有限公司的(https://www.vazyme.com/Home.html)试剂Taq Pro Universal SYBR qPCR Master Mix进行定量分析。每个孔总体系为20 μL:1 μL模板cDNA,0.4 μL上游引物,0.4 μL下游引物,10 μL 2×Tap Pro Universal SYBR qPCR Master Mix,8.2 μL无菌水。每个样品有3个技术重复。扩增反应程序设置为:95 ℃预变性30 s;95 ℃ 10 s,60 ℃ 30 s,循环40次;95 ℃ 15 s,60 ℃ 60 s,95 ℃ 15 s,循环1次。用2-ΔΔCT法[37]计算候选基因相对于内参基因的表达量。

表1 引物设计信息

2 结果与分析

2.1 转录组数据分析

2.1.1 转录组数据比对结果

将不同胁迫下的处理组:水淹3 d(WS-3 d)、水淹15 d(WS-15 d)、干旱3 d(DS-3 d)、干旱5 d(DS-5 d)以及对照组(CK)各3个重复的总共15个转录组数据的原始数据(raw reads)进行质控,得到Q20和Q30的均值分别为97.59%和93.67%和平均GC含量为43.71%和测序深度均在30×以上的clean reads。将其比对到组装好的大叶杨参考基因组上,平均比对率为93.82%(表2)。以上结果表明转录组测序数据质量较高,可用于后续分析。利用StringTie统计每个样本中基因的TPM值,为了检验每组数据3个个体间的重复性,基于皮尔逊相关系数对不同处理组中3个生物重复的表达量数据绘制相关性热图(图1)。结果显示每个处理组内的3个个体间,其相关性都在0.9及以上,表明各处理组内的3个样本具有良好一致性,可以进行下一步分析。

图1 大叶杨对照和胁迫样本相关性热图

表2 转录组样本数据概述

2.1.2 差异表达基因鉴定

利用DESeq2分别将处理组与对照组进行差异分析。干旱处理组:在DS-3d与CK组鉴定到5 747个DEGs,其中包含3 399个上调和2 348个下调的DEGs;在DS-5 d与CK组鉴定到8 486个DEGs,其中上调和下调的DEGs分别有4 035,4 451个(图2,A,C)。水淹处理组:在WS-3 d与CK组共获得2 114个DEGs,其中1 455个上调和659个下调;在WS-15 d与CK组共获得1 077个DEGs,上调和下调的DEGs分别有572,505个(图2,B,C)。随着水淹胁迫的天数增加,DEGs的数量由水淹3 d的2 114个逐渐减少到1 077个;而随干旱胁迫天数增加,DEGs数量明显增加。接下来筛选出在WS-3 d和WS-15 d共享的DEGs作为全程响应水淹胁迫的DEGs(命名WS-DEGs),同理也筛选出全程响应干旱胁迫的DEGs(命名为DS-DEGs),分别鉴定出385个WS-DEGs和3 986个DS-DEGs(图2,A,B)。对二者取交集,其中有237个WS-DEGs和DS-DEGs共享的DEGs(命名为Shared DEGs),还有148个WS-DEGs和3 749个DS-DEGs作为独立响应水淹和干旱胁迫的DEGs(分别命名为WS-uniq DEGs和DS-uniq DEGs),见图2,D。大叶杨中响应干旱胁迫的差异基因数量明显多于水淹胁迫时差异基因数量。

2.1.3 差异表达基因GO富集

为了探究大叶杨在水分胁迫下的生物学过程及其响应胁迫的大致机理,将上述鉴定到的Shared DEGs、WS-uniq DEGs和DS-uniq DEGs分别进行基于Fisher检验的GO富集。Shared DEGs富集到的生物过程主要包括原花青素生物合成途径、调控叶片衰老、对水分刺激的响应、光周期以及开花调控、获得系统抗性、免疫防御等(图3,A);在WS-uniq DEGs中,富集到色氨酸代谢、生长素介导的信号通路正调控、芳香族氨基酸代谢过程、吲哚类化合物的代谢过程等通路(图3,B);在DS-uniq DEGS中,主要富集到光合通路、叶绿素代谢、获得系统抗性、应激防御反应等(图3,C)。依次以不同类型DEGs(Shared、WS-uniq、DS-uniq)的数量作为随机取样值,对大叶杨基因组中全部基因进行20 000次随机抽样,统计每次抽样数据集中AP2/ERF家族成员的数量,比较后发现AP2/ERF家族成员显著富集于差异表达基因中。观测值表示不同类型DEGs中AP2/ERF基因的数量(图3)。

2.2 大叶杨AP2/ERF家族成员鉴定和系统发育关系

参考Nakano等[18]在拟南芥中对AP2/ERF基因家族成员分类的方法,将大叶杨中的候选基因分为5类,包括AP2、ERF、DREB、RAV和Soloist。在大叶杨中,有65个(31%)DREB亚家族基因,103个(50%)ERF亚家族基因,30个(14%)AP2亚家族基因,6个RAV亚家族基因和1个Soloist亚家族基因,其中ERF亚家族中的基因数量最多,Soloist亚家族中只有1个Polas10384基因。利用NJ法则对205个大叶杨AP2/ERF家族成员和138个拟南芥AP2/ERF家族成员进行蛋白序列的多重比对,构建系统发育树(图4)。在前人研究中[39],将ERF与DREB亚家族统称为ERF亚家族,二者主要的区别在于AP2/ERF结构域中第14位和第19位氨基酸残基位置不同。因此可以看出,大叶杨中DREB与ERF亚家族成员之间亲缘关系更近。

圆圈代表后续qPCR实验中挑选的大叶杨AP2/ERF基因;黄色和蓝色分别代表干旱和水淹胁迫。

2.3 qRT-PCR实验验证

根据大叶杨AP2/ERF基因家族成员在水淹和干旱胁迫下的转录水平变化,一共筛选出12个在水淹和干旱胁迫下基因表达存在明显差异的基因进行qRT-PCR验证(图5、图6)。

A.AP2/ERF基因在干旱胁迫下基因聚类热图; B.AP2/ERF基因在水淹胁迫下基因聚类热图。

前6个为AP2/ERF基因在干旱胁迫下相对表达量,后6个为AP2/ERF基因在水淹胁迫下相对表达量。

所有候选基因在系统发育树中的位置见图4。结合试验结果和基因热图发现Polas35282、Polas00804、Polas05763、Polas28115在干旱处理下具有明显下调趋势(图5,A和图6),其中Polas28115同时在水淹胁迫下也明显下调(图5,B,图6)。qPCR实验结果表明基因表达趋势与转录组数据分析的结果一致,这不仅进一步证实了转录组数据的可靠性,也表明筛选的AP2/ERF转录因子家族成员基因在受到水分胁迫时发生差异表达,可能在大叶杨水淹和干旱胁迫下发挥重要的调控作用。

3 讨 论

本研究运用高通量测序技术对水淹、干旱胁迫和对照的大叶杨样本进行转录组测序,得到测序质量较好且相关性系数均在0.9以上的转录组数据。通过转录组数据分析比较大叶杨在不同水分条件刺激下差异基因的变化,发现干旱胁迫下DEGs的数量远高于水淹胁迫,推测与大叶杨长期生长在水资源丰富地区有关,其体内可能已经进化出一系列应对水淹、洪涝等非生物胁迫的适应性机制。随着干旱胁迫天数增加,DEGs的数量明显增加。相反的是,持续水淹胁迫并没有导致DEGs数量的增加。WS-uniq DEGs GO显著富集在生长素介导的信号通路等,而生长素对植物新生根系的形成具有重要的调控作用,研究表明施加生长素与吲哚丁酸能够促进苹果幼苗产生新生根[40],因此推测大叶杨在水淹胁迫后期形成的不定根可能与生长素调控机制有关,而不定根的形成能够缓解植物根系因水淹胁迫而产生的不良影响[41],这表明大叶杨可能具有较高的耐水淹能力。DS-uniq DEGs的GO结果显著富集在光合作用通路、叶绿素代谢等过程,参与该过程的基因表达均呈下调趋势。因此进一步推测干旱胁迫降低了大叶杨光合作用的能力,这与先前研究提出的杨树在干旱胁迫下,光合速率和蒸腾速率显著下降的结果[42]一致。

AP2/ERF是植物中特有的转录因子家族,参与了许多植物的非生物胁迫响应,它主要通过激素调节方式应对非生物胁迫。在拟南芥[18]、毛果杨等[43]模式植物中已经对AP2/ERF基因家族进行了研究。发现在拟南芥中At4g13040是Soloist家族中唯一成员,正向调控水杨酸的积累且能够防御细菌病原体的入侵[25];玉米中的ZmEREB180属于ERF亚家族,通过乙烯的诱导表达,从而促进不定根的形成来增强玉米的抗涝性[44],而目前尚未对大叶杨中的AP2/ERF基因家族进行探究。本研究鉴定了大叶杨中的AP2/ERF家族成员,一共有205个基因属于该家族。根据基因表达聚类热图及qPCR实验结果表明,发现Polas28115(AtRAP2.4)在水淹和干旱胁迫下均显著下调,而RAP2.4在拟南芥中参与干旱胁迫应答,在提高植物抗旱性中有重要作用[45-46],但尚未有研究揭示其在水淹胁迫下的调控作用;Polas35282(ATRAP2.7)随着干旱胁迫的天数逐渐增加,基因表达量明显下调,在四倍体棉花的盐胁迫和干旱胁迫研究中,发现RAP2.7的同源基因GhirPA0801G001500为差异表达基因[47]。该结果说明大叶杨Polas35282可能参与了干旱胁迫的应答过程,至于其调控机理还需进一步深入研究。

目前极端气候事件频发正在导致动植物遭受越来越多的非生物胁迫压力,如气温升高、极端高温事件频繁发生、土壤干旱、土地沙漠化严重。干旱和水淹胁迫作为植物育种生产的一大难题,可能导致作物减产、林木等植物种群的地理分布和生存能力受限。因此探究植物胁迫响应的分子基础并筛选识别与抗逆性有关的转录因子是植物非生物胁迫研究和作物育种的关键。本研究初步探究了大叶杨在水淹和干旱胁迫下基因的差异表达情况,同时对广泛参与非生物胁迫的AP2/ERF转录因子家族成员进行了筛选与鉴定,未来可以联合基因组学、群体遗传学进一步对大叶杨的适应气候变化潜力和适应性进化机制进行探究,为大叶杨分子遗传育种和种质资源保护制定更为完善的方案,这也将为其他林木树种的遗传保护提供更多的思路及参考。