一、研究背景

然而，目前干旱胁迫下毛果杨(Populus trichocarpa)的表观转录组与蛋白质组间的关系尚不清楚。在该项研究中，作者利用纳米孔直接RNA测序技术（Nanopore Direct RNA-seq ,DRS）和基于串联质谱标签(Tandem mass tag ,TMT)的蛋白质组学来分析研究干旱胁迫下茎分化木质部(SDX)的表观转录组和蛋白质组学调控,揭示了干旱胁迫下毛果杨SDX中N6-甲基腺苷（m6A）、多聚腺苷酸化、poly(A) tail length（PAL）和翻译之间的相互作用。

二、研究手段

研究对象：3个月的毛果杨，取对照组为无干旱处理（ND）毛果杨苗，实验组为5天(D5)和7天(D7)干旱处理后的SDX进行DRS和TMT蛋白组学分析。DRS做2个生物学重复，TMT蛋白组学分析做3个生物学重复。

研究方法：Nanopore DRS；TMT蛋白组学分析；RT-qPCR；LC-MS/MS

三、研究结果

1.干旱处理后，全长转录本的比例下降

DRS检测全长和非全长reads，结果显示，非全长reads显示出极端的3'偏好 (图1C)；在干旱处理后，全长reads的百分比下降(图1D、E); D5 vs ND和D7 vs ND分别有546个和2004个基因全长reads降低，而全长reads增加的基因很少(图1F、G)，干旱处理后全长reads比增加的基因在ND中的全长reads比非常低(20%)，而全长reads比较低的基因在ND中显示出较高的全长reads比(40%)(图1H)。

图1. 干旱胁迫下SDX的纳米孔DRS检测结果。A, 利用DRS和TMT标记分析干旱胁迫引发的表转录组和蛋白质组变化流程图。Nanopolish、PRAPI和Nanom6A是三种用于长 read分析的管道。B, 热图显示ND, D5和D7样本在全长和非全长读数之间的Spearman相关性。C，转录本全长和非全长读段的覆盖率。D,柱状图显示了干旱胁迫下全长read比的动态变化。E,Wiggle plot, histogram, and RT-qPCR显示干旱处理降低了DUF239的全长读取率。差异显著性采用t检验和利用星号评级 (*P < 0.05; **P < 0.01). Total_F/Total_R和Full_F/Full_R是引物位置。F,散点图显示干旱处理后不同的全长比。G，维恩图显示了D5与ND和D7与ND比较中具有增加的全长读取比的基因重叠。H，箱线图显示ND、D5和D7样品中不同全长读比基因的比例分布。box代表第25和75个百分位数;中位置线，第50个百分位数;晶须，第10百分位和第90百分位。圆圈表示离群值。

2. 干旱胁迫下，差异表达蛋白（DEPs）与差异表达基因(DEGs)存在显著重叠

ND、D5和D7的蛋白质组分析共鉴定出9872个蛋白对应的86904个肽段，与ND相比，D5与D7分别鉴定出802个DEPs(491个上调，311个下调)和821个DEPs(467个上调，354个下调)，然而D5和D7之间的DEP很少(图2A、B)，这表明干旱处理引起的蛋白质丰度变化在D5阶段已经完成。对DEPs进行GO分析，发现，在干旱处理后，上调DEPs (D5 vs ND)显著富集的GO项是“抗氧化活性”、“谷胱甘肽过氧化物酶活性”和“胚胎发育”，而上调的DEPs属于与纤维素生物合成相关的蛋白质(图2C)。进一步研究干旱胁迫下蛋白质丰度与mRNA水平之间关系的一致性，发现大约61.3%(301/491)和68.2%(212/311)的上调和下调蛋白在D5与ND的mRNA水平上被鉴定为差异表达，大约64.9%(303/467)和62.1%(220/354)的上调和下调蛋白在D7与ND的mRNA水平上被鉴定为差异表达(图2D、E)。

图2. 干旱胁迫下SDX的质谱蛋白质组学研究。A，火山图分别显示了D5与ND、D7与ND、D5与D7的DEPs。B, DEPs热图。深蓝色和浅黄色分别代表高和低蛋白质丰度。C, D5与ND的两两比较分别显示，分布在上调(左)和下调(右)蛋白中富集了GO term。D和E，维恩图显示D5与ND (D)和D7与ND (E)中DEPs和DEGs间的重叠。

3. 干旱胁迫下，蛋白质丰度和mRNA表达的相关性下降

文中通过对DRS进行TMT标签技术研究了干旱胁迫下全长和非全长RNA的mRNA表达和蛋白丰度（Rp）的Pearson相关系数，并对得到的数据做归一化处理，结果显示，随着标准基因表达阈值的增加，全长RNA和非全长RNA的Rp均下降(图3, A和B)，这与之前其发表的PacBio Iso-Seq data的数据趋势相同，此外，作者还观察到干旱处理下RNA和蛋白质之间的相关性降低，这表明在干旱胁迫下RNA表达和蛋白质丰度的变化是不同的(图3, A和B)。

为了比较全长RNA和非全长RNA的Rp 差异，作者利用GO term对全长RNA和非全长RNA进行下游分析，结果显示D5(图3D)和D7(图3E)的ND（图3C）非全长的RNA相比全长RNA均具有更高的Rp，为了进一步探讨干旱胁迫下Rp 的变化,作者对ND的D5和D7的GO terms重叠区进行分析，结果发现全长RNA中的Rp 在6个GO term中有提高，而在17个GO term中有下降(图3F)，Rp在涉及半胱氨酸型肽酶活性的term中增加，Rp 在与氧胁迫关联的term中下降。经过干旱胁迫后的非全长的RNA，Rp在5个GO term上升，在34个GO term中下降(图3G)。Rp 在具有蛋糖基化功能term中升高，在氧化应激和其他应激反应中降低(图3G)。

图3. 干旱胁迫下蛋白质丰度与mRNA表达的相关性。A\ -B，全长(A)和非全长RNA (B)不同RPM阈值水平的Rp 线形图。C - E, ND (C)， D5 (D)和D7 (E)中基于GO term的全长RNA和非全长RNA间不同Rp的散点图。F，热图显示干旱处理下全长RNA Rp 的变化。heatpap的颜色(深蓝色到红色)表示相关性(从低到高)。G，热图显示干旱处理下非全长reads Rp变化。

4. 干旱胁迫下m⁶A的动态变化

对m⁶A谱的转录组分析显示，所有样本的m⁶A峰主要分布在编码序列和3'-UTR中(图4A、B)。在干旱胁迫下，平均m⁶A比值增加（图4C），其中木材形成基因的m⁶ A比例增加，并伴有mRNA(图4E)和蛋白(图3F)的表达下降，表明m⁶A可能在胁迫下抑制木材形成中发挥作用。研究者进一步鉴定了差异表达的m⁶A (DE m⁶A)位点。与ND相比，D5共有1048个m⁶A降低，775个m⁶A升高，而D7与ND相比，1295个位点的m⁶A下降，717个位点的m⁶A增加(图3G)，结合转录组数据分析DE m⁶A位点与基因表达之间的关系，在D5和ND中鉴定出457个具有DE m⁶A位点的DEG，其中，在mUP-gDN（m⁶A比例增加，基因表达减少）类别中(图4I，象限II)，干旱胁迫降低了13个成木基因的表达，m⁶A比例增加，这与m⁶A甲基化通常与转录丰度呈负相关是一致的。此外，干旱胁迫下SDX中的m⁶A比率与蛋白质丰度和H3K9ac水平之间的关系表现出与基因表达相似的趋势(图4J、K)。

图4. 干旱胁迫下m⁶A甲基化的动态响应。A, m⁶A甲基化位点在转录本的分布。B，柱状图显示ND、D5和D7中包含UTR和CDS的注释基因区域中m⁶A位点的百分比。C，箱线图显示ND, D5和D7的m⁶A比率。D - F，箱形图显示ND、D5和D7木质形成基因的m⁶A比率(D)、基因表达(E)和蛋白水平(F)。框限代表第25和75个百分位数;中线，第50个百分位数;晶须，第10百分位和第90百分位。G，散点图显示D5(左)和D7(右)干旱处理后m⁶A变化差异(mod)。H, 转录因子WLIM1的差异m⁶A修饰。I - K, D5与ND两两比较的散点图显示了m⁶A比例与基因表达(I)、蛋白水平(J)和H3K9ac水平(K)的关系。

5. 干旱胁迫促进远端poly(A)位点的使用

CPSF30控制选择性聚腺苷酸化（APA）位点的选择。在干旱胁迫下，Ptr-CPSF30的表达上调，并且Ptr-CPSF30的近端poly(A)位点利用率(近端/远端和近端，PPR)下调(图5A)。发现在干旱胁迫下，AAUAAA和1-nt变异的poly(A)信号比略有增加，而UCGUUU和1-nt变异的poly(A)信号比略有下降(图5B、C)。与ND相比，D5和D7的PPR降低的基因更多，揭示了干旱胁迫下远端poly(A)位点的使用增加(图5D、E)。对显示PPR升高和PPR降低的基因分别进行GO富集分析，发现PPR升高的基因参与“纤维素合成酶活性”和“纤维素生物合成过程”，PPR降低的基因参与“蛋白去磷酸化”和“蛋白酪氨酸/丝氨酸/苏氨酸磷酸酶活性”(图5F-I)。

图5.在响应干旱胁迫中，poly(A)位点使用的动态变化。A，柱状图显示经RT-qPCR验证的干旱处理下CPSF30的表达增加和PPR降低。差异显著性采用t检验，采用星号评分系统(*P < 0.05;** p < 0.01)。B和C，注状图显示干旱胁迫下AAUAAA(和1-nt变化)(B)和UCGUUU(和1-nt变化)(C)六聚体基序的百分比。D，散点图显示PPR响应干旱胁迫变化。E, Wiggle图和柱状图显示IBR5的PPR下降(Potri.014G160500)。F/H，GO富集分析中基因分别的PPR的增多（F）和PPR的下降展示。G/I, Wiggle图和柱状图显示CESA7和VAB2作为PPR增加(G)和PPR降低(I)的基因的例子。

6.干旱胁迫诱导不同的PAL

利用DRS产生的全长poly(A)尾序列，确定ND、D5和D7样本的PAL中位数分别为79、80和75 nt(图6A)。PAL中位数最长的10个GO具有“转录共调节活性”和“激素应答”功能，PAL中位数最短的10个GO具有翻译和核糖体组成功能(图6B)。干旱胁迫下PAL增加最多的5个基因功能组(D5>ND)均与蛋白质加工有关，而PAL减少最多的5个基因功能组(D5<nd)的主要功能是酶活性和对激素的响应(图< span="">6C)。内含子保留（IR）reads具有更长的poly(A)（图6D）。

图6. 差异PAL对干旱胁迫的响应。A，小提琴图显示ND, D5, D7的PAL中位数。水平虚线表示PAL为80 nt。B，箱形图显示了ND中PAL最长的10个GOterm(橙色)和PAL最短的10个GO项(蓝色)。box代表第25和75个百分位数;中线，第50个百分位数;晶须，第10百分位和第90百分位。虚线表示PAL为80 nt。C，箱形图显示基因功能组在D5与ND(上)和D7与ND(下)中PAL(左)增加最多，PAL(右)减少最多。A，表示“激酶活性”;B，“核苷酸切除修复”;C,“蛋白质去泛素化”;D，“蛋白质代谢过程”;E, G蛋白偶联受体信号通路;F，“激素反应”;G，“n -乙酰转移酶活性”;H,“类异戊二烯生物合成过程”;1、“蛋白丝氨酸/苏氨酸激酶活性”;J，“含核碱基化合物代谢过程”;K，“金属离子跨膜转运体活性”;L,“大核糖体亚基”;M，“蛋白磷酸酶2A型调节剂活性”;N，“蛋白磷酸酶2A型复合体”;O代表“微管运动”;P，“微管运动活性”;Q，“驱动蛋白复合物”。水平虚线表示80 nt的PAL。D,箱形图表示完全拼接(no_IR)和IR转录本的PAL分布。虚线表示no_IR和IR转录本中位PAL。

7. PAL与转录水平和翻译相关

根据转录物丰度将基因分为三类:高、中、低丰度，发现高表达基因倾向于较短的PAL，而低表达基因倾向于较长的PAL(图7,A-C)。高表达基因倾向于富集最优密码子，并在ND、D5和D7中分别表现出对短PAL的强烈偏好(图7,D-F)。通过累积分数分析，观察到具有低、中、高水平最佳密码子的转录本在PAL分布、mRNA丰度和蛋白质丰度方面存在显著差异(图7,G-I)。然而，无论是干旱处理还是干旱处理，PAL和蛋白质丰度之间都没有实质性的关联(图7,J-L)。综上所述，这些结果表明高表达的转录本对短PAL有强烈的偏向。

图7. 高表达转录本倾向于较短的PAL。A - C为PAL分布于ND (A)、D5 (B)、D7 (C)基因中的不同表达水平。最高、中、最低分别代表500个表达量最高的基因、500个表达量最接近中位数的基因和500个表达量最低的基因。P-value采用曼-惠特尼U检验计算。D - F，热图表示ND (D)、D5 (E)和D7 (F)中最优密码子频率(Fop)和PAL之间的关系。G - I，累积总和图显示翻译效率(密码子优化)与PAL (G)、mRNA丰度(H)和蛋白丰度间(I)之间的关系。J - L，PAL在ND (J)、D5 (K)、D7 (L)不同蛋白水平池中的分布情况。p值采用曼-惠特尼U检验计算。

四、结论

本研究揭示了干旱胁迫引发的毛果杨SDX表观转录组和蛋白质组变化的全面概况。在干旱处理下，全长读比降低，m⁶A刺激纤维素和木质素相关基因的RNA降解。此外，干旱胁迫促进了包括CPSF30在内的许多基因向远端poly(A)位点的转移。最后，PAL与全长读取率和m⁶A修饰率有关。总之，这些发现极大地扩展了关于干旱对基因表达转录后影响的理解。

编辑 ：qingqin H审核：Tao Li

参考文献

Yubang Gao, Xuqing Liu, Yandong Jin, Ji Wu, Shuang Li, Yaxing Li, Binqing Chen, Yaxin Zhang, Linxiao Wei, Wei Li, Ruili Li, Chentao Lin, Anireddy S N Reddy, Pankaj Jaiswal, Lianfeng Gu, Drought induces epitranscriptome and proteome changes in stem-differentiating xylem of Populus trichocarpa. Plant Physiology. 2022 Sep, 190(1):459–479.