PEG和NaCl胁迫下毛竹萌发种子中环状RNA特征及其表达研究

王晓静,王 涛,杨 凯,李潞滨

(1. 中国林业科学研究院林业研究所,林木遗传育种国家重点实验室,国家林业和草原局林木培育重点实验室, 北京 100091;2.运城学院,山西 运城 044011;3.北京植物园,北京市花卉园艺工程技术研究中心/ 城乡生态环境北京实验室,北京 100093;4.北京农学院植物科学与技术学院,北京 100096)

毛竹(Phyllostachysedulis)是禾本科(Gramineae)竹亚科(Bambusoideae)刚竹属(Phyllostachys)的多年生木本散生竹,是我国种植面积最广、应用及研究最为广泛和深入的竹种,具有重要的经济、生态和文化价值[1-2]。使用毛竹种子进行实生苗造林是毛竹造林的重要方式,具有竹鞭生长较快、萌笋多、竹材利用率高、运输方便、成本低等优点,有利于提高造林整齐度、成活率和遗传多样性[3-5]。但毛竹种子寿命短,其种子萌发容易受到环境的影响,限制了实生苗造林及其研究[6-12]。因此,开展毛竹种子萌发相关研究对毛竹种质资源的保存和应用具有重要意义。

随着全球气候的极端变化,干旱和盐胁迫等非生物胁迫变得越来越频繁,增加了植物生长和发育面临的挑战。前期研究结果表明,干旱和盐胁迫是抑制毛竹种子萌发的重要环境因素[13],研究者已经对干旱和盐胁迫下毛竹种子的萌发率、活力指数、氧化酶活性以及萌发后幼苗的生长开展了研究[6, 10, 12],但对分子水平上的相关调控机制仍然缺乏了解。因此,亟须挖掘调控毛竹种子萌发和萌发期抗旱抗盐的因子,为竹类资源的遗传改良、培育和应用提供依据。

circRNAs是一类由转录本反向可变剪接产生的内源非编码RNA,在真核生物中非常丰富[14]。circRNAs的环状结构相对稳定,能够参与调节线性转录本表达、选择性剪接以及与miRNA或蛋白质的相互作用等[15]。近年来对植物中与circRNAs的研究逐渐兴起,目前PlantcircBase数据库(Version 7.0,更新于2022-03-15)已经收录了拟南芥(Arabidopsisthaliana)、水稻(Oryzasativa)、玉米(Zeamays)、小麦(Triticumaestivum)等21种植物的171 118个circRNAs[16]。据报道circRNAs参与了植物生长、发育和胁迫响应[17-18]。Zhou等[19]通过CRISPR-Cas9技术发现水稻circRNAOs05circ02465的无效突变体在种子萌发期间表现出较高的耐盐性,表明其能够参与胁迫下种子萌发的调控。竹类植物中,Wang等[20]对毛竹茎秆快速生长阶段的circRNAs进行了鉴定和分析,Li等[21]报道了circRNAs在3月龄的毛竹幼苗中对干旱胁迫的响应。然而,竹类植物中与circRNAs相关的研究还存在大量空白,特别是在干旱或盐胁迫下circRNAs如何参与毛竹种子萌发调控研究鲜见报道。

本研究使用聚乙二醇(PEG)和NaCl模拟干旱和盐胁迫,通过构建链特异性文库首次对毛竹萌发种子中的circRNAs进行鉴定,对PEG和NaCl胁迫下毛竹种子萌发阶段circRNAs的表达进行分析,并对其源基因进行功能富集,以期为深入开展竹类植物种子在非生物胁迫条件下萌发的分子调控机制研究奠定基础。

1 材料与方法

1.1 试验材料

试验所用毛竹种子于2020年9月由广西壮族自治区林业科学研究院收集,剥除外稃后挑选健康饱满且大小均匀的种子备用(图1a)。种子千粒质量为(21.85±0.04)g,含水率为(12.05±0.35)%。将种子表面消毒后转移到底部有两层无菌滤纸的无菌培养皿中,每个培养皿放置20粒种子。试验分为5组(A1—E1):分别在培养皿中加入H2O(A1)、10%(质量分数,下同) PEG(B1)、15% PEG(C1)、50 mmol/L NaCl(D1)和100 mmol/L NaCl(E1),每组萌发200粒种子,重复3次。在(26 ± 2) ℃的持续黑暗环境下进行萌发,每2 d更换滤纸和培养液。观察和记录种子萌发状态,在第4天选择生长状态一致且到达种皮破裂阶段的种子用于建库(图1b)。对于每个处理,在3个重复中各取5粒种子并将15粒种子混样,使用液氮速冻后置于-80 ℃保存。

a)毛竹干种子dry seed;b)种皮破裂阶段取样的萌发种子sampled seed germinated at seed coat rupture stage。图1 实验用毛竹种子Fig. 1 Moso bamboo seeds used in the experiment

1.2 建库、测序与数据过滤

使用植物多糖多酚RNA提取试剂盒(DP441,天根生物科技有限公司,中国)提取样品的总RNA,RNA质检合格后首先使用Ribo-ZeroTM试剂盒(MRZSR116,Epicentre,中国)去除rRNA,随后将RNA片段化并合成cDNA第1条链,在cDNA第2链合成时用dUTP取代dTTP。使用尿嘧啶特异性切除酶消化含有U的cDNA第2链,然后对第1链连接接头后进行PCR扩增和富集得到最终的链特异性文库,并对其进行初步定量。使用Illumina高通量测序平台NovaSeq 6000(美国)对质检合格的链特异性文库进行双末端测序,测序长度为150 bp。对测序后的原始数据进行过滤,过滤掉接头序列、含有≥1%的未知碱基的序列、低质量(≤20)碱基数≥50%的序列。

1.3 circRNA鉴定

使用find_circ[22]和CIRI2[23]鉴定circRNAs,find_circ查找接头序列与线性序列位置相反,且内部存在GT-AG剪接位点的序列;CIRI2寻找具有目的信号和GT-AG剪切位点信息的连接点处的读序(junction reads)序列。然后根据与参考基因组的比对结果[24-25],对鉴定结果取交集,并保留至少1个样品中junction reads数≥5的circRNAs。

1.4 表达水平分析

使用具有circRNA剪接位点的junction reads数目估计circRNAs的表达量,将表达量进行SRPBM(spliced reads per billion mapping)归一化定量[26]。使用edgeR进行差异表达分析[27],差异表达circRNAs筛选条件为差异倍数(foldchange)以2为底的对数绝对值≥1且P<0.05。使用TBtools绘制表达量热图[28]。

1.5 源基因及其功能富集分析

使用HISAT2将circRNAs比对至毛竹染色体水平参考基因组,鉴定其源基因[24-25]。对源基因进行GO和KEGG注释,P<0.05为富集结果显著性阈值[29-30],使用ImageGP在线平台展示富集结果[31]。

2 结果与分析

2.1 测序数据产量

本研究得到lncRNA原始数据如表1所示,过滤后共产出14.66～15.28 Gb数据,clean reads的Q20值均≥98%,Q30值均≥95%,GC含量为52.4%～53.7%,将过滤后数据进行比对,77.70%～84.98%的数据可以比对到基因组,表明数据产量较高,且质量较好,能够满足后续分析。测序原始数据提交至NCBI数据库(PRJNA545401)。

表1 文库测序数据统计

2.2 circRNA的鉴定

将两个软件的鉴定结果取交集,结果在本研究不同处理下萌发的毛竹种子中共鉴定到1 446个circRNAs,其长度分布范围为173～99 647 bp,每个样本中长度≤2 000 bp的circRNAs频率均为最高(图2a)。1 446个circRNAs主要分布在毛竹染色体scaffold 1—24,分布在scaffold 21和scaffold 17的数目最多,而分布在scaffold 1、scaffold 2、scaffold 19染色体的circRNAs最少,此外有少量(11个)circRNAs分布在基因组未组装片段中(图2b)。1 446个circRNAs中有76.34%来源于基因组外显子区,19.23%来源于基因间隔区,来源于内含子区的circRNAs最少(4.42%)。每个样本分布在基因组不同来源区域circRNAs的数量见图2c,结果显示每个样本中来源于外显子区的circRNAs数量最多,只有少量circRNAs来源于内含子区。

a)5个样本circRNAs长度分布与频率the length and frequency of circRNAs in five samples;b)circRNAs在毛竹染色体的分布distribution of circRNAs in moso bamboo chromosomes;c)circRNAs在基因组的分布distribution of circRNAs in genome。图2 毛竹萌发种子中circRNAs的特征Fig. 2 Characteristics of circRNAs in germinated seeds of moso bamboo

2.3 源基因鉴定

对circRNAs的源基因进行分析,结果表明,产自外显子区和内含子区的1 168个circRNAs共来源于867个基因,有698个源基因产生了1个circRNA,176个源基因产生了2个及以上circRNAs,14个源基因产生了5～13个circRNAs。如胼胝质合成酶基因PH02Gene21577(基因编号,后同)(callose synthase 9,CALS9)产生了13个circRNAs,核孔锚蛋白基因PH02Gene01074(nuclear-pore anchor,NUA)和CHD3型染色质重塑因子PICKLE编码基因 PH02Gene11133(CHD3-type chromatin-remodeling factor PICKLE,PKL)产生了9个circRNAs,胼胝质合成酶基因PH02Gene00919(callose synthase 9,CALS9)产生了8个circRNAs。其中PKL能够参与渗透胁迫和外源对种子萌发的调控,pkl突变体在脱落酸诱导下ABI3和ABI5过量表达,其启动子区H3K9和H3K27甲基化水平比野生型大大降低,但不影响其在萌发后期幼苗的甲基化水平,同时PKL能够沉默转基因逆境诱导型RD29A启动子的转录沉默[32-33],但目前关于circRNAs在PKL调控中的功能仍是未知的,毛竹萌发期种子中PKL基因产生的circRNAs是否能够参与PEG和NaCl胁迫下种子萌发的调控仍需进一步研究。

2.4 表达特性分析

对照组(A1)毛竹萌发种子中表达量最高的circRNA为plant_circ_0001703,其源基因编号为PH02Gene18257,编码油质蛋白(olein)。10% PEG(B1)、15% PEG(C1)、50和100 mmol/L NaCl(E1)胁迫下毛竹萌发种子中表达量最高的circRNA分别为plant_circ_0001728、plant_circ_0000358、plant_circ_0001236、plant_circ_0001728,其源基因分别为球蛋白编码基因PH02Gene28357、水通道蛋白TIP3-1编码基因PH02Gene14831、功能未知基因PH02Gene04721和PH02Gene28357。对照组或胁迫处理下毛竹萌发种子中表达量前10的circRNAs热图及其源基因的注释见图3,其中包含了大量种子储藏蛋白编码基因,如球蛋白基因、vicilin-like种子储藏蛋白基因、油质蛋白基因等,表明在毛竹种子萌发过程中种子储藏蛋白基因产生并积累了大量的circRNAs,但这部分circRNAs是否能够调控毛竹种子萌发仍需要进一步研究。

差异表达circRNAs(图4)可见,与对照组相比,10%与15% PEG、50与100 mmol/L NaCl胁迫下差异表达circRNAs(differential expressed circRNAs,DECs)数量分别为524、505、467和474,其中有122个circRNAs在4种胁迫下均显著差异表达;此外,有503个circRNAs在PEG处理B1-C1比较组中差异表达,465个circRNAs在不同浓度NaCl胁迫下差异表达(D1—E1);6个比较组中共有1 056个circRNAs差异表达。上述circRNAs在不同胁迫条件下萌发的毛竹种子中存在表达量的显著变化,推测可能参与了胁迫下毛竹种子萌发的调控。

2.5 差异表达circRNAs源基因功能富集分析

图5 不同比较组差异表达circRNAs源基因GO显著富集前20结果Fig. 5 Top20 significantly enriched GO terms for original genes generating DECs in differential comparision groups

图6 不同比较组差异表达circRNAs源基因KEGG显著富集结果Fig. 6 Significantly enriched KEGG pathways of original genes generating DECs in different comparison groups

为了进一步探究PEG和NaCl胁迫下毛竹萌发期种子中circRNAs的响应,对差异表达circRNAs(DECs)的源基因进行GO和KEGG功能富集。6个比较组中显著富集前20的GO条目见图5,其中作用于酸酐的水解酶活性(hydrolase activity, acting on acid anhydrides,GO:0016817)在6个比较组中均被显著富集到前20,推测其对毛竹萌发期种子中circRNAs对PEG和NaCl的响应可能具有重要意义。KEGG富集结果(图6)表明,6个比较组中显著富集途径数目分别为4(A1-B1)、3(A1-C1)、2(A1-D1)、5(A1-E1)、5(B1-C1)、1(D1-E1)。KEGG途径中鞘脂代谢途径(sphingolipid metabolism)在A1-B1、A1-C1、A1-E1 3个比较组均被显著富集,推测鞘脂代谢途径是circRNAs响应PEG和NaCl的重要途径。此外,过氧化酶体途径,脂肪酸降解途径,基本转录因子途径在A1-C1和B1-C1比较组中被显著富集,淀粉与蔗糖代谢在A1-C1比较组中被显著富集,糖胺聚糖降解途径在A1-E1、D1-E1比较组中均被显著富集,ABC transporters在A1-E1比较组中被显著富集,这些途径与种子萌发密切相关,推测在PEG或NaCl胁迫下毛竹种子萌发时,上述途径可能通过相关circRNAs响应胁迫,并调控萌发进程。

3 讨 论

本研究在对照组、PEG和NaCl胁迫处理下萌发的毛竹种子中共鉴定1 446个circRNAs,对其特征分析表明这些circRNAs来源主要为外显子区,只有少部分circRNAs来源于内含子区,这与毛竹幼苗叶片中circRNAs的基因组来源区域是一致的,在拟南芥、玉米、水稻、小麦中也发现了相似的结果[34-37];同时本研究中circRNAs在染色体的分布是不均匀的,这与毛竹幼苗叶片中circRNAs的分布结果也是一致的。

已有研究[8,35,38-39]表明,植物中circRNAs参与了对干旱、盐胁迫的响应。在玉米368个自交系中,circRNAs的源基因富集在与circRNAs表达和耐旱性相关的SNP,表明了circRNAs在玉米干旱响应中的重要作用[8]。拟南芥中过表达由保卫细胞外向整流K+通道基因产生的circGORK导致转基因植株对脱落酸高度敏感,但对干旱不敏感,表明circGORK在耐旱性中发挥了积极作用[35]。在紫花苜蓿(Medicagosativa)、鹰嘴豆(Cicerarietinum)和大豆(Glycinemax)中,大量的circRNAs能够响应干旱或盐胁迫[38-39]。本研究在不同比较组中也鉴定了大量响应PEG和NaCl胁迫处理的差异表达circRNAs,这些circRNAs对于干旱和盐胁迫下的毛竹种子萌发可能具有重要调控作用。

对不同比较组差异表达circRNAs的源基因进行GO富集分析,结果表明作用于酸酐的水解酶活性的相关基因在6个比较组中均显著富集;涉及到的源基因包括CHD3型染色质重塑因子PICKLE编码基因(PKL)PH02Gene11692、PH02Gene11133、PH02Gene11133,ABC转运蛋白编码基因(ABC transporter C,ABCC)PH02Gene23726、PH02Gene32020,钙离子转运ATP酶基因(Calcium-transporting ATPase,ACA)PH02Gene33073、PH02Gene34105等,其中PH02Gene11133是本研究中产生circRNAs最多的源基因。前人的研究中PKL、ABCC和ACA均参与了种子萌发的调控[33, 40-41],推测作用于酸酐的水解酶活性相关基因可能对circRNAs调控干旱和盐胁迫下的毛竹种子萌发具有重要意义。

根据差异表达circRNAs源基因的KEGG富集分析,鞘脂代谢途径在A1-B1、A1-C1、A1-E1中均被显著富集,推测其可能参与了PEG和NaCl胁迫下circRNAs对毛竹种子萌发的调控。鞘脂代谢途径参与了玉米种子萌发过程中病原真菌诱导的细胞程序性死亡[42]和拟南芥种子萌发期对冷胁迫的响应[43],但目前关于鞘脂代谢途径参与种子萌发和萌发期抗逆的研究相对较少,鞘脂代谢途径是否能够通过circRNAs参与PEG和NaCl胁迫下的种子萌发仍需进一步研究。本研究中DECs源基因中还涉及基础转录因子,基础转录因子在A1-C1,B1-C1中显著富集,涉及的7个基因中6个属于转录起始因子TFIID(transcription initiation factor, TFIID),推测TFIID可能参与了circRNA对毛竹种子萌发的调控。前人研究表明TFIID在干旱胁迫下对水稻生长和光合作用的调节可能起着关键作用[44]。在拟南芥中TAF10参与构成TFIID复合物,atTAF10的过表达提高了种子萌发期间对盐胁迫的耐受性,而敲除突变体对盐胁迫更为敏感,表明转录起始因子是植物盐害毒性的生理靶标复合体[45]。但目前关于TFIID通过circRNAs进行基因调控的研究还未见报道,TFIID是否能够通过产生circRNAs调控种子萌发还需要进一步研究。

本研究首次对毛竹萌发期种子中的circRNAs进行了初步鉴定和胁迫响应分析,丰富了毛竹circRNAs数据,并为胁迫下毛竹种子萌发调控因子的挖掘及相关研究提供参考;后续将进一步开展circRNAs功能的验证工作,并对其在毛竹种子萌发中的调控机制进行深入研究。