通过植株表型分析与统计比较两种水稻材料的生长差异：成熟期野生型zyk639的株高显著高于zyk-lm（图1a、d）；zyk-lm的千粒重显著低于zyk639（图1c）；但zyk-lm的有效穗数和每穗粒数均高于野生型zyk639（图1e、f）。从苗期到成熟期，zyk-lm叶片均表现出类病斑表型：倒三叶上出现大量红棕色病斑，倒二叶也有少量病斑，而旗叶几乎无病斑（图1b）。
图1.野生型（zyk639）与突变体（zyk-lm）的植株形态及农艺性状比较。a，成株期植株表型。白色标尺=20厘米。b，成株期叶片表型。黑色标尺=2厘米。每组图像从左至右依次为：左侧3片叶分别为zyk639的旗叶、倒二叶、倒三叶；右侧3片叶分别为zyk-lm的旗叶、倒二叶、倒三叶。黑色标尺=2厘米。c，千粒重。d，株高。e，单株有效穗数。f，每穗粒数。数值为10株植株的平均值±标准差（SD），**表示在p<0.01水平上差异极显著。对zyk-lm和zyk639接种10个（P2、P3、P4、P5、P6、P7、P8、P9、P10、C4）xoo菌株（Xanthomonas oryzae pv. Oryzae，水稻白叶枯病的病菌），观察并统计各叶位病斑长度（图2a）。结果显示，易感材料zyk639在接种P6和P10菌株后病斑较长，且病斑长度按旗叶、倒一叶、倒三叶的顺序递增（图2b）；而zyk-lm在接种10个白叶枯病菌株后均未表现发病表型。
图2.野生型（zyk639）与突变体（zyk-lm）的抗病性比较。a，接种10个菌株后的白叶枯病病斑长度。b，3个叶片部位的发病表型。黑色标尺=2厘米。对三个叶位进行台盼蓝（TB）染色发现，与zyk639相比，zyk-lm叶片的蓝黑色斑点数量显著增加，且病斑周围叶肉区域也出现小蓝斑，表明zyk-lm病斑部位发生细胞死亡并影响周围细胞，且随着叶片衰老死亡现象加剧；而zyk639三个叶位几乎无蓝斑，说明其叶片未发生细胞死亡（图3a）。3,3’-二氨基联苯（DAB）染色结果显示，zyk639三个叶位均无明显褐色斑点，而zyk-lm叶片出现明显褐色斑点，且倒三叶的褐色标记大于倒二叶和旗叶（图3b），这与TB染色结果一致，推测可能是zyk-lm中ROS积累引发细胞PCD，进而导致类病斑表型。
图3.六叶期zyk639和zyk-lm植株叶片的DAB染色与台盼蓝（TB）染色结果。a，DAB染色结果（从左至右）：依次为zyk639的旗叶、倒二叶、倒三叶，以及zyk-lm的旗叶、倒二叶、倒三叶。b，台盼蓝（TB）染色结果（从左至右）：依次为zyk639的旗叶、倒二叶、倒三叶，以及zyk-lm的旗叶、倒二叶、倒三叶。
采用Illumina深度测序技术对18个叶片样本进行测序，主成分分析（PCA）显示重复样本间表达相关性较高，证实数据重复性良好，可支持后续分析。以P<0.05且FC≥1.5为筛选标准，共获得11,930个DEmRNAs（差异表达基因）。9个对比组（zyk639a vs zyk-lma、zyk639b vs zyk-lmb、zyk639c vs zyk-lmc、zyk639a vs zyk639b、zyk639a vs zyk639c、zyk639b vs zyk639c、zyk-lma vs zyk-lmb、zyk-lma vs zyk-lmc、zyk-lmb vs zyk-lmc）的DEmRNAs数量分别为4,478、7,199、7,170、1,566、2,307、981、6,640、7,416、1,779个，其中zyk-lma vs zyk-lmc组的DEmRNAs数量最多（图4a）。火山图显示，不同对比组中上调与下调基因的分布模式清晰（图4b）。
图4.各比较组中差异表达mRNA的分析。mRNA表达水平的火山图，红色点代表上调表达基因，黄色点代表下调表达基因，蓝色点代表无显著差异基因。
差异表达长链非编码RNA（DElncRNAs）鉴定
共鉴定出5,289个lncRNAs，包括95个已知lncRNAs和5,194个新预测lncRNAs。根据新lncRNAs在基因组上相对于蛋白质编码基因的位置，可将其分为5类：基因间区lncRNA（intergenic lncRNAs）、双向lncRNA（bidirectional lncRNAs）、内含子lncRNA（intronic lncRNAs）、反义lncRNA（antisense lncRNAs）、正义重叠lncRNA（sense overlapping lncRNAs）（图5a）。共获得1,873个DElncRNAs，在9个对比组中，zyk639c vs zyk-lmc组的DElncRNAs数量最多；除zyk-lmb vs zyk-lmc和zyk639a vs zyk639b组外，其余7个组的上调基因数量均显著多于下调基因（图5b）。
图5.差异表达长链非编码RNA（DElncRNAs）的鉴定。a，长链非编码RNA（lncRNAs）的类型。b，长链非编码RNA（lncRNAs）中显著差异调控基因的数量。采用RNAplex预测反义lncRNA与mRNA的互补关系，共获得883对lncRNA-mRNA靶基因对，包括874个lncRNAs和768个mRNAs，其中仅9个lncRNAs对应2个mRNA，其余均为一对一关系。对具有顺式调控功能的lncRNAs进行mRNA靶基因预测，共鉴定出12,739对lncRNA-mRNA靶基因对，包括4,569个lncRNAs和9,214个mRNAs，其中超过70%的预测lncRNAs可靶向2个以上mRNA（图6）。通过分析样本间lncRNAs与蛋白质编码基因的表达水平相关性（共表达分析）进行反式调控靶基因预测，从18个样本中获得39,213对lncRNA-mRNA靶基因对，包括549个lncRNAs和4,500个mRNAs，其中52%的lncRNAs可靶向10个以上mRNA，另有976个mRNA仅对应1个lncRNA。此外，发现WRKY家族转录因子Os11g0490900可靶向53个lncRNAs，而MSTRG.3292.1可靶向2个R2R3-MYB家族转录因子，提示水稻PCD相关lncRNAs可能通过调控转录因子影响PCD发生。
图6.长链非编码RNA（lncRNAs）的顺式（cis）和反式（trans）靶基因分析。a，与mRNA存在顺式（cis）关系的lncRNAs预测数量。b，与lncRNA存在顺式（cis）关系的mRNA预测数量。c，与mRNA存在反式（trans）关系的lncRNAs数量。d，与lncRNA存在反式（trans）关系的mRNA数量。
为深入探索lncRNA与水稻叶片PCD的潜在调控关系，对zyk639和zyk-lm的lncRNA转录本与mRNA转录本的基因表达矩阵进行筛选，共获得16,485个基因用于WGCNA分析。通过构建基因聚类关系将其分为18个模块，分析各模块中基因表达谱与所有样本的相关性（图7a），并绘制模块-样本矩阵热图（图7b）。结果显示，蓝色模块中的基因表达谱在PCD发生与PCD未发生样本间存在显著差异（图7c）。蓝色模块包含5,870个基因，主要与蛋白质合成相关。对蓝色模块中lncRNA靶向的共表达mRNA进行GO和KEGG富集分析：GO功能富集结果显示，这些mRNA显著富集于“细胞大分子代谢过程（GO0044260）”、“含氮化合物代谢过程（GO0006807）”、“结合（GO0005488）”、“细胞（GO0005739）”、“细胞内（GO0005622）”、“细胞器（GO0043226）”等条目（图8a）；KEGG通路富集结果显示，主要涉及“核糖体（ko03010）”、“剪接体（ko03040）”、“RNA转运（ko03013）”、“真核生物核糖体生物发生（ko03008）”、“DNA复制（ko03030）”等通路（图8b）。
图7.zyk639与zyk-lm中鉴定基因的加权基因共表达网络分析（WGCNA）。a，基于拓扑重叠差异层次聚类得到的差异表达基因（DEGs）树状图。b，展示各模块中样本表达模式的热图。c，蓝色模块（module blue）中所有共表达基因的表达谱。颜色梯度代表Z值，柱状图代表该模块中对应共表达基因的共有表达模式。图8.蓝色模块中长链非编码RNA（lncRNA）靶向的共表达mRNA的GO富集分类与KEGG通路富集分析。a，mRNA在GO中显著富集（P<0.05）。b，涉及前20个术语的KEGG通路。Q值范围为0-1，Q值越接近0，富集程度越显著；Rich Factor值越大，表明富集程度越高。
为进一步明确这些DEGs间的关系，构建基因共表达网络。选取蓝色模块中与核糖体合成和DNA复制相关的富集基因，绘制权重值前100的差异基因间关系图（图9）。该网络中，rpl27-1（Os08g0404300）、OsRPS27等关键基因在核糖体通路中显著富集。
图9.蓝色模块的共表达调控网络分析。
取蓝色模块中的mRNA和lncRNA与DElncRNAs、DEmRNAs取交集，结合miRNA构建ceRNA网络（采用Cytoscape软件可视化）。共鉴定出101对调控关系，包括9个lncRNAs、11个miRNAs和76个mRNAs，构建8个复杂竞争关系网络（图10）。其中3个关键lncRNAs（MSTRG.24301.2、MSTRG.17476.1、MSTRG.24898.1）可分别结合miR8175-z、miR11033-z、osa-miR5809；MSTRG.21792.1可与novelm0237-5p协同调控Os02g0557000、Os03g0850400、Os06g0157000、Os08g0512400；MSTRG.5369.2可与miR477-x协同调控Os01g0308300、Os02g0220500、Os03g0773300、Os06g0609700；MSTRG.31242.3可与miR5225-z、miR3627-z协同调控Os05g0393400、Os07g0666300；MSTRG.16622.1可与miR5257-z协同调控Os03g0749800。值得注意的是，Os01g0566100基因可同时被MSTRG.7176.1和MSTRG.28624.1调控（图10）。此外，发现WRKY家族转录因子Os09g0334500和JMJ家族转录因子Os09g0483600可在MSTRG.24898.1与osa-miR5809结合后调控靶基因（图10）。
图10.基于蓝色模块中筛选出的差异表达长链非编码RNA（DElncRNAs）、差异表达信使RNA（DEmRNAs），并结合信使RNA（mRNAs）、长链非编码RNA（lncRNAs）与微小RNA（miRNAs）构建的内源竞争RNA（CeRNA）网络。
随机选取18个基因进行qRT-PCR验证，结果显示，qRT-PCR检测的基因转录水平与RNA-seq估算的表达丰度高度一致（图11），表明本研究获得的RNA-seq数据可靠，可支持上述转录组分析结果。
图11.实时荧光定量聚合酶链式反应（qRT-PCR）的表达谱与转录组测序（RNA-Seq）分析结果的比较。
将预测的候选基因与GFP融合后转入水稻9311原生质体，通过激光共聚焦显微镜观察目标基因在细胞内的定位。构建4个差异基因（PAC，Os05g0393400；PTD，Os05g0588200；ELF3-2，Os01g0566100；Q2R432，Os11g0490900）的亚细胞定位载体，激光共聚焦显微镜观察显示：PTD定位于细胞质，ELF3-2定位于细胞核（图12），提示多个细胞器可能参与水稻PCD的触发与调控。
图12.选定基因的亚细胞定位分析。
已有研究表明，植物中ROS和H₂O₂的过量积累可引发PCD。因此，本研究通过检测zyk639和zyk-lm中的过氧化物含量、抗氧化酶活性及渗透调节物质含量，探究水稻PCD中的抗氧化系统变化。如表1和表2所示，zyk639与zyk-lm不同叶位的抗氧化酶活性和渗透调节物质含量存在显著差异（p<0.05）：zyk-lm的SOD、POD、CAT活性及MDA含量均高于zyk639，且zyk-lmb的抗氧化酶活性高于zyk-lma和zyk-lmc；O₂⁻和H₂O₂含量变化趋势与抗氧化酶活性基本一致；zyk-lma的脯氨酸含量为35.97μg/g，显著高于其他样本，分别是zyk-lmb的2.22倍、zyk639c的2.13倍；而zyk639的可溶性蛋白含量高于zyk-lm，其中zyk639a的可溶性蛋白含量是zyk-lmc的1.44倍；zyk-lmb的可溶性糖含量最低，为17.92mg/g（表2）。
表1.zyk-lm对抗氧化活性的影响。
不同小写字母表示差异显著（P<0.05）。
表2.zyk-lm对渗透调节物质含量的影响。
不同小写字母表示差异显著（P<0.05）。