刘溶荣,陈炜曦,尹济容,李姿燕,王季春,3,4,荐红举,3,4*,吕典秋,3,4*
(1.西部(重庆)科学城种质创制大科学中心,西南大学,重庆 400715;2.西南大学农学与生物科技学院,重庆 400715;3.西南大学南方山地农业教育部工程研究中心,重庆 400715;4.西南大学薯类生物学与遗传育种重庆市重点实验室,重庆 400715)
近年来,极端天气频发,干旱作为全球变暖带来的十大危害之一,严重影响了农业生产[1]。马铃薯(Solanum tuberosumL.)作为浅根系作物,对水分十分敏感,干旱在马铃薯各个生育时期都有可能发生,且在各时期对马铃薯不同品种会造成不同程度的影响,最终影响其产量和品质[2-4]。
miRNA是一类长度为18~25 nt的内源性非编码单链小分子RNA,越来越多的证据表明,miRNA 参与调控植物对非生物胁迫的响应[5]。干旱胁迫上调了拟南芥中miR156、miR159、miR167、miR168、miR171、 miR172、 miR319、 miR393、 miR394a、miR395c、miR395e、miR396 和miR397 的表达水平,降低了miR161、miR168a、miR168b、miR169、miR171a 和miR319c 的表达水平[6,7]。但在某些情况下,miRNA 对干旱胁迫的响应取决于物种。例如,干旱胁迫下miR156 表达丰度在水稻[8]和玉米[9]中下调,而在拟南芥[6]、大麦[10]和野生二粒小麦[11]中上调。在拟南芥中,miR319-TCP4通过介导茉莉酸和生长素的合成途径参与干旱胁迫[12]。在抗旱能力不同的烟草品种中,miR398 在干旱胁迫后呈现出不同的表达模式,但是均可抑制其靶标铜锌超氧化物歧化酶基因(Cu/Zn-SOD,CSD)的表达[13]。Yang等[14]预测了干旱胁迫下马铃薯中脯氨酸积累相关的miRNAs,发现miR172、miR396a、miR396c和miR4233 调控吡咯啉-5-羧酸合成酶基因(P5CS),miR2673 和miR6461 分别调控吡咯啉-5-羧酸还原酶基因(P5CR)和脯氨酸脱氢酶基因(ProDH);Zhang 等[15]通过深度测序对马铃薯干旱胁迫下保守和新型的miRNAs 进行了鉴定,确定了4 个调控干旱相关基因(miR811、miR814、miR835、miR4398)。鉴定到3 个stu-miR159 成员(stu-miR159a、 stu-miR159b 和stu-miR159c)均在马铃薯干旱胁迫处理25 d后表达量显著降低,预测的3 个靶基因GAMyb-like 家族成员(StGAMyb1、StGAMyb-like2.1和StGAMyb-like2.2)表达量显著增加[16]。因此,关于响应马铃薯干旱相关的miRNA及其靶基因还有待深入挖掘。
本研究旨在通过20% PEG6000 模拟干旱处理,对0,1,3,6,12,24和48 h共7个时间点的马铃薯组培苗进行sRNA文库构建和测序,利用生物信息学分析的方法,鉴定响应干旱胁迫的miRNA及其靶基因,探究miRNA 及其靶基因在马铃薯干旱胁迫下的调控作用,为后期miRNA 的功能验证奠定基础,同时为阐明马铃薯参与干旱响应的分子机制提供理论依据。
以马铃薯品种‘克新18 号’组培苗为试验材料,供试材料为重庆市薯类生物学与遗传育种重点实验室保存。
在Murashige & Skoog(MS)培养基(3%蔗糖)进行组培苗的生长,14 d后移出,炼苗2 d后使用1/2 Hoagland 营养液水培5 d(16 h 光照/8 h 黑暗;22℃/18℃)。
使用20% PEG6000 替换营养液,在0,1,3,6,12,24 和48 h 分别取相同部位叶片,各处理2次生物学重复,每个生物学重复使用3 株叶片混合,液氮速冻后于-80℃保存。
将上述14 个样本送至百迈客生物科技有限公司,提取RNA,质检合格后,使用Illumina HiSeq 2500平台进行sRNA文库构建和测序。对获得的原始测序数据进行质量控制,去除未知碱基N(N 为无法识别的碱基)含量大于等于10%的和没有3′接头序列的Reads,剪切掉3′接头序列,获得18~30个核苷酸的高质量值序列(即Clean reads)。
利用Bowtie 软件,将Clean reads 分别与Silva数据库、GtRNAdb数据库、Rfam数据库和Repbase数据库进行序列比对,过滤核糖体RNA(rRNA)、转运RNA(tRNA)、核内小RNA(snRNA)、核仁小RNA(snoRNA)等非编码RNA(ncRNA)以及重复序列,获得包含miRNA 的Unannotated reads。利用Bowtie 软件将Unannotated reads 与马铃薯杂合二倍体参考基因组(RH89_039_16_v2)进行序列比对,获取在参考基因组上的位置信息(即Mapped reads)。
将Mapped reads与miRBase(v22)数据库中已知miRNA的成熟序列及其上游2 nt与下游5 nt的范围进行比对,最多允许一个错配,鉴定到已知miRNA(Known-miRNAs);利用miRDeep2 软件包,通过Reads 比对到基因组上的位置信息得到可能的前体序列,基于Reads 在前体序列上的分布信息(基于miRNA产生特点、Mature、Star、Loop)及前体结构能量信息(RNAfold randfold),采用贝叶斯模型经打分最终实现新miRNA(Novel-miRNAs)的预测。基于序列的相似性使用miRDeep2软件对检测到的已知miRNA和新miRNA进行miRNA家族分析。
对各样本中miRNA 进行表达量的统计,并用TPM 算法对表达量进行归一化处理。TPM 归一化处理公式为:
式中:Read count 表示比对到某一miRNA 的Reads数目。
使用皮尔森相关系数对样品间相关性进行计算,计算公式为:
使用|log2(FC)|≥0.58;P-value≤0.05 作为筛选标准,采用Benjamini-Hochberg 校正方法对原有假设检验得到的显著性P值(Pvalue)进行校正,使用DESeq2 进行样品组间的差异表达分析,获得差异表达miRNA集。
根据Known-miRNAs 和Novel-miRNAs 与RH89_039_16_v2的基因序列信息,使用TargetFinder软件进行靶基因预测。使用BLAST 软件将预测靶基因序列与NR、Swiss-Prot、GO、COG、KEGG、KOG和Pfam数据库比对,获得靶基因的注释信息。
利用富集因子(Enrichment factor)分析差异表达miRNA 靶基因在某一通路上的富集程度,并利用Fisher 精确检验方法计算P值,经FDR 矫正P值后得到Q值(Qvalue)。其中富集因子的计算公式为:
富集因子=(单个Pathway中的差异表达miRNA靶基因数/单个Pathway 中的所有靶基因数)/(所有Pathway中的差异表达miRNA靶基因数/所有Pathway中的所有靶基因数)
选取富集显著性最可靠(即Qvalue 最小)的前20 个通路进行结果展示,得到差异表达miRNA 靶基因的KEGG通路富集结果,利用联川生物云平台(https://www.omicstudio.cn/home)进行作图展示。
使用EASY spin 植物microRNA 快速提取试剂盒(艾德莱,中国)对20%PEG6000 处理0,3,12和48 h 的叶片提取包含miRNA 的总RNA;使用该公司的增强型miRNA 反转录试剂盒进行miRNA 3′末端的PolyA加尾和逆转录反应。
参照加PloyA 尾法[17],对筛选出的差异表达miRNA 设计正向荧光定量PCR 引物,以马铃薯18S rRNA 为内参基因。miRNA 上游引物见表1,下游为北京艾德莱生物科技公司荧光定量检测试剂盒中提供的通用引物。反应体系及反应程序参照试剂盒说明书。
表1 miRNA荧光定量PCR上游引物序列Table 1 Forward primer sequences for miRNA qPCR
每个miRNA 选择两个靶基因,通过Spud DB(http://spuddb.uga.edu/)获得靶基因的CDS 序列,使用Beacon Designer 7.9 软件进行定量引物设计(表2),以马铃薯EF1α作为内参基因。每个样品设3 次生物学重复、3 次技术重复,采用2-ΔΔCt算法计算基因的相对表达量。
利用Excel Microsoft 365 软件进行数据整理分析,使用GraphPad Prism 9软件进行作图展示。
对14 个样品的sRNA 测序原始数据进行质控,从各样品中获得了高质量值序列(表3)。14个样品中功能未注释的序列数在12 628 838(63.86%)~16 227 826 个(77.15%),rRNA、tRNA、snRNA、snoRNA 等ncRNA 以及重复序列占比较少。一般在植物中rRNA 低于60%被认为是合格品,本测序结果的14 个样本中rRNA 数目在4 314 282(20.51%)~6 763 846个(34.20%),说明sRNA文库的构建是合格的,可以进行下一步分析。
表3 小RNA的种类和数量Table 3 Categorization of unique reads and total reads of small RNA
将未注释的序列与马铃薯杂合二倍体参考基因组(Solanum_tuberosumRH89_039_16_v2)进行序列比对,获取在参考基因组上的位置信息(表4)。各样本比对到参考基因组上的Reads数在6 651 941(56.47%)~8 612 998个(60.93%),比对到正链上的Reads数占所有Unannotated reads的比例在39.54%~42.48%,比对到负链上的Reads 数占比在16.93%~18.53%。总体来看,14 个样本中一半以上的未注释序列可以比对到参考基因组上,且大多数序列比对到基因组正链上。
表4 参考基因组比对信息Table 4 Distribution of sequence reads mapped to genome
14 个样品共得到621 个miRNAs,其中已知miRNAs(Known-miRNAs)308 个,新预测miRNAs(Novel- miRNAs)313 个。 CK、 D1、 D3、 D6、D12、D24 和D48 样品中分别发现了286、291、275、272、269、272和265个Known-miRNAs,CK和D1 中Known-miRNAs 的数量明显多于其他处理,随着处理时间的延长,Known-miRNAs的数量逐渐减少;各处理中发现的Novel-miRNAs 数量相同,均为313个(图1A)。不同大小的miRNAs在各处理中分布相似,大都在20~24 nt 分布。且Known-miRNAs 和Novel-miRNAs 的长度都集中在21 和24 nt,但Known-miRNAs 在长度为21 nt 的序列含量最高,其次为24 nt 的序列,并且干旱胁迫处理3 h 后,21 nt 的Known-miRNAs 数目明显减少;Novel-miRNAs 则相反,在长度为24 nt的序列含量最高,其次为21 nt;在22 和23 nt 的KnownmiRNAs和Novel-miRNAs的序列数目相似(图1B)。
图1 miRNAs鉴定结果Figure 1 Results of identified miRNAs
miRNA在物种间具有高度保守性,基于序列的相似性对检测到的Known-miRNA 和Novel-miRNA进行miRNA 家族分析,研究miRNA 在进化中的保守性。250 个已知miRNAs 属于69 个家族,166 个新预测miRNAs 属于59个家族。总共鉴定到104个miRNA家族,其中有24个家族是Known-miRNA和Novel-miRNA所共有的(表5和图2A)。
图2 miRNA家族分析Figure 2 Analysis of miRNA family
表5 miRNA家族统计Table 5 Summary of miRNA family
miR399家族含有的Known-miRNA最多,有30个,其次为miR398 含有15 个Novel-miRNA。在处理后,共有16 个家族的miRNA 发生了变化,其中miR399 家族的miRNA 在干旱处理3 h 时明显减少,miR7982 家族的miRNA 在干旱处理6 h 后未检测到,miR7993和miR8005家族的miRNA也在干旱处理后出现不同程度的减少,miR169_1 家族的miRNA 在干旱3 h 前数目增加,3 h 后开始减少,另外,miR158 和miR6425 家族的miRNA 在干旱处理后才被检测到(图2B)。
对鉴定到的621 个miRNAs 进行差异表达分析,共鉴定到40 个差异表达miRNA,其中上调表达的12个,下调表达的28个(图3A)。干旱处理前6 h,上调的miRNA 逐渐减少,6 h 时全为下调表达的miRNA,6 h 后逐渐增加;干旱处理3 h 时,下调的miRNA减少,之后增加,干旱处理3 h时的差异miRNA数目最少(图3B)。
图3 响应干旱胁迫miRNAs差异表达分析Figure 3 Analysis of differentially expressed miRNAs in response to drought stress
对筛选出的40 个差异表达miRNA 做层次聚类分析,将具有相同或相似表达行为的miRNA 进行聚类,总共分为7 大类。Type Ⅰ为下调表达的miRNA,其中2 个Known-miRNA(stu-miR482a-5p和stu-miR391-3p),7 个为Novel-miRNA,他们大多在干旱胁迫6 h 及以后表达量下降;Type Ⅱ为先下调后上调表达的miRNA,其中有3 个KnownmiRNA,1 个Novel-miRNA,他们在干旱胁迫3~12 h 被下调,12 h 后开始上调表达;Type Ⅲ为先下调后上调再下调表达的miRNA,2 个都为Known-miRNA(stu-miR162a-5p 和stu-miR162b-5p),他们在干旱胁迫3~6 h 下调表达,12 h 被上调,12 h 后又被下调;Type Ⅳ为表达无明显规律的miRNA,总体呈上升趋势,包括5 个KnownmiRNA 和2 个Novel-miRNA;TypeⅤ为先下调后上调再下调再上调的miRNA,共13 个,大多数都是Novel-miRNA,只有1 个Known-miRNA(stumiR8016);TypeⅥ为上调表达的miRNA,2 个都是 Known- miRNA(stu- miR477b- 3p 和 stumiR477b-5p),他们在干旱胁迫6 h 开始逐渐被上调表达;TypeⅦ为先上调后下调表达的miRNA,包括1 个Known-miRNA(stu-miR477a-5p)和2 个Novel-miRNA,他们在24 h 后逐渐开始下调表达(图3C)。
靶基因的预测对于理解miRNA 的功能至关重要。鉴定到的621个miRNAs中有555个miRNAs被预测到了共9 222 个靶基因,其中分别有280 个Known-miRNAs 和275 个Novel-miRNAs 被预测到了6 889个和3 010个靶基因(表6)。对靶标基因只有1 个和2 个的miRNA 分析发现,Known-miRNA中同一家族的miRNA 靶基因相同,多个不同的Novel-miRNA 靶标同一基因,推测靶标同一基因的不同Novel-miRNA可能属于同一家族。
表6 miRNA靶基因数目预测统计Table 6 Prediction statistics of the number of miRNA target genes
对差异表达的40个miRNAs的靶基因进行分析,共有16个Known-miRNAs和21个Novel-miRNAs分别预测到558 和631 个靶基因。联合靶基因预测结果和转录组数据,对上述miRNA 所对应的靶基因与转录组中差异表达的基因取交集,最终分别对16 个已知(表7)和17 个新miRNA(表8)的81 个和70个靶基因进行了功能注释。
表8 差异表达新的miRNAs及其靶基因Table 8 Differential expression of novel-miRNAs and their target genes
由表7可以看出,靶基因中大多是与植物生长代谢有关的蛋白基因,例如:细胞色素P450、GDSL 脂酶、NADH 依赖性谷氨酸合成酶和UDP 糖基转移酶超家族蛋白。参与马铃薯干旱胁迫响应及氧化还原的一些蛋白和酶类有漆酶、2-氧戊二酸和铁(II)依赖性加氧酶超家族蛋白、氧化还原酶家族蛋白、锌指(C3HC4型环指)家族蛋白和CDPK相关激酶等。还有一些与激素相关基因,例如:赤霉素生物合成所需基因、类生长素抗性基因和乙烯反应/乙烯调节核蛋白(ERT2)。同时,预测到一些参与植物生长发育及胁迫响应的转录因子,例如:TCP 家族转录因子、转录因子bHLH61和转录调节因子。另外,还有与植物昼夜节律相关蛋白TIMLESS 和一些保守的未被表征的假定蛋白。从注释结果来看,一个miRNA 可靶标多个基因,不同miRNA也可能调控同一基因。
新的miRNA 预测出的靶基因参与马铃薯干旱胁迫的大多为各种蛋白和酶类。例如:ARM 重复序列超家族蛋白、含有P-环的核苷三磷酸水解酶超家族蛋白、2-氧戊二酸和铁(II)依赖性加氧酶超家族蛋白、受体凝集素激酶(RLK)、真核天冬氨酰蛋白酶家族蛋白、转录因子B3 家族蛋白、蛋白磷酸酶2C 家族蛋白、bHLH DNA 结合超家族蛋白以及EIN3结合F-box蛋白等(表8)。没有预测到参与胁迫响应的直接靶标转录因子。但从预测结果来看,属于同一家族的新miRNA 的靶基因同源性较高。
在生物体内,不同的基因相互协调来行使生物学功能。将差异表达miRNA靶基因的通路注释进行富集分析,选取富集显著性最可靠(即Qvalue最小)的前20个通路进行结果展示(图4)。差异表达miRNA的靶基因主要富集的前五个KEGG通路有谷胱甘肽代谢、倍半萜和三萜生物合成、其他类型的O-聚糖生物合成、脂肪酸伸长率和类固醇生物合成。其中,谷胱甘肽具有抗氧化和整合解毒作用,可防止自由基、过氧化物、脂质过氧化物和重金属等活性氧对重要细胞成分的损伤,在植物应对干旱胁迫中发挥重要作用。差异表达miRNA靶基因在苯丙烷代谢通路富集,苯丙烷代谢是重要的植物次生代谢途径之一,其代谢产物,如木质素、类黄酮、花青素和有机酸等,在调控植物生长发育及对非生物胁迫响应中发挥重要功能。在干旱胁迫下,植物通过释放α-亚麻酸重塑细胞膜的流动性,进而应答干旱胁迫。本研究中,差异表达miRNA靶基因通路也在α-亚麻酸代谢途径富集。另外,还在细胞信号转导重要组分鞘脂代谢通路、内质网中的蛋白加工通路、氨基酸生物合成以及基础转录因子等通路中富集。
图4 差异表达miRNA靶基因KEGG通路富集Figure 4 Enrichment of KEGG pathway of target genes with differential expressed miRNA
为了检测测序结果的可靠性,选取0,3,12和48 h 4 个处理时间点,使用qRT-PCR 对6 个miRNA(包括3个已知和3个新的miRNA)进行表达模式分析。定量表达结果与测序数据的miRNA 表达趋势一致(图5),表明测序结果真实可靠。
图5 miRNA测序及qRT-PCR结果Figure 5 miRNA sequencing and qRT-PCR results
miRNA对其靶标转录后调控是通过对靶标剪切来实现的,因此为了验证miRNA对靶标表达的调控作用,本研究对上述miRNA靶标基因的表达变化进行了qRT-PCR的验证及分析。大多数靶标的表达表现出与miRNA 相反的趋势,例如:stu-miR482a-5p、novel_miR_159和novel_miR_110在干旱胁迫后表现为先下调后上调,其靶基因StCLC-B和StAST、StALPL和StCYP以及StPP2C的相对表达量则表现出先升高后下降的趋势,novel_miR_110 的靶基因StACL5在干旱处理48 h 内的表达量呈一直上升趋势;stu-miR393-3p在干旱胁迫后表现为下调,而其靶基因StCHP1和StCHP2表现为上调;stu-miR391-3p 在干旱胁迫前12 h 表现为下调,在48 h 略有上升,其靶基因StELGP和StGBF的相对表达量则在前12 h 上升,之后逐渐下降;novel_miR_283 的相对表达量在前12 h下降,之后升高,其靶基因StPrp4的相对表达量在干旱48 h 内呈上升趋势,而另一靶基因StpsbY的相对表达量则呈下降趋势(图6)。
图6 miRNA及其靶基因的相对表达量Figure 6 Relative expression of miRNA and their target genes
马铃薯为浅根系作物,其整个生育时期都有可能受到干旱胁迫的影响。目前,植物对干旱胁迫的响应多集中在离子响应、脱落酸(Abscisic acid,ABA)依赖和非依赖途径以及转录因子等研究[3,18,19],关于miRNA 调控马铃薯干旱胁迫响应的报道较少[14-16,20]。因此,挖掘马铃薯中响应干旱胁迫的miRNA,进一步解析miRNA 及其靶基因之间的调控作用,为后期miRNA 的功能验证奠定基础,同时为阐明马铃薯参与干旱响应的分子机制提供理论依据。基于此,本研究通过20%PEG6000模拟干旱处理0,1,3,6,12,24 和48 h 共7 个时间点的14个马铃薯组培幼苗样品进行small RNA测序。测序结果中,Known-miRNA与Novel-miRNA的长度分别在21 和24 nt 的序列最多(图1),这与Zhang 等[15]、李莹等[21]、Xu 等[22]和Wu 等[23]在马铃薯、枣树、拟南芥和紫色马铃薯块茎中的报道一致。推测miRNA 的长度分布与物种、取材部位以及处理条件有关。此外,还发现长度为21 nt 的Known-miRNA 在干旱胁迫3 h 后明显减少(图1),说明马铃薯在干旱胁迫3 h时即对miRNA进行大量调控。
miR319 被报道其参与植物叶片发育、花器官形成、花粉发育、衰老、昼夜节律和激素信号转导等多种生长发育过程的调控[24]。本研究中,stumiR319-3p 在干旱胁迫24 h 时显著上调(图3),与拟南芥[7]以及紫花苜蓿[25]中在根和叶中都上调表达结果一致。靶基因预测有TCP家族转录因子、蛋白质丝氨酸/苏氨酸激酶和碳酸氢盐转运蛋白家族等,这与前人报道的miR159/319 家族的miRNA 通过调控MYB/TCP 转录因子来响应干旱胁迫一致[24,26]。但是,关于马铃薯中miR319是否及如何与其他激酶和转运蛋白等调控还有待进一步研究。
miR398和miR408家族在植物中较为保守,参与植物生长、发育和胁迫响应的调节,在植物遭到胁迫时对植物的生长发育过程发挥重要作用。在本研究中,stu-miR398a-3p 和stu-miR408a-3p 的表达量在PEG 处理1 h 后即出现显著上升(图3),这与蒺藜苜蓿、拟南芥中miR398和miR408上调表达一致[27,28]。stu-miR408a-5p 的表达量在PEG 处理24 h后出现显著下调(图3),这与水稻干旱胁迫后miR408 下调表达一致[8]。因此,干旱胁迫下miR398和miR408的表达模式变化与物种、品种之间抗性以及成员有关。对stu-miR398a-3p、stumiR408a-3p 和stu-miR408a-5p 的靶基因预测结果显示(表7),包括细胞色素P450(CYP450)、漆酶(LAC)、GDSL 脂酶以及四三肽重复序列样超家族蛋白等,这与前人已报道的拟南芥miR408 的靶基因有Laccase和蒺藜苜蓿中miR398 的靶基因有COX5b有相似之处。漆酶是一种含铜的多酚氧化酶,其在防御反应以及木质素和植物次生壁合成过程起重要作用。已有报道表明,AtLAC4过量表达植株在干旱处理后耐旱能力比野生型明显增强[29]。因此,在后续的研究中,可以进一步关注miRNA及其靶基因LAC在马铃薯干旱胁迫响应中的调控机理。也表明在干旱胁迫过程中,同一靶基因可能受到多个miRNA同时调控以应答胁迫环境。
miR482是植物中一个古老且相对保守的miRNA超家族,在植物的生长发育和抗病抗逆过程中发挥重要的作用[30]。在本研究中,stu-miR482a-5p在干旱处理6 h后被显著下调(图3),与miR482在大豆中下调一致[31];但stu-miR482e-3p在胁迫处理48 h后被显著上调(图3),与盐胁迫下miR482.2 在杨树[32]和干旱胁迫下miR482 在大豆[33]中被上调一致。在靶基因预测中,stu-miR482a-5p 靶向基因有转录调节因子(Transcription regulators,TRs)、类萌发素蛋白(Germin-like proteins,GLPs)以及氧化还原酶家族蛋白(Oxidoreductase family protein)等;而stu-miR482e-3p则主要靶向含NB-ARC结构域的抗病蛋白(NB-ARC domain-containing disease resistance protein)(表7)。关于该家族成员是如何协同植物各类生长发育以及代谢相关蛋白和转录因子,去抵御干旱及其他非生物胁迫,还有待进一步研究。
尽管,本研究鉴定出一些响应马铃薯干旱胁迫的miRNA,他们主要通过调节生物合成、生长代谢、氧化还原反应等过程来响应胁迫环境。但是,挖掘出的这些miRNA 是通过何种调控途径来响应干旱胁迫的分子机制研究,还需进一步通过转基因等手段验证。本研究结果为后续深入研究miRNA调控马铃薯干旱胁迫响应的分子机制提供线索,为抗旱种质资源的筛选和抗旱品种的选育提供理论依据。