不同浓度氮处理毛竹水通道蛋白基因表达模式及其调控网络

朱成磊 林泽铭 李天阔 高志民

(国际竹藤中心竹藤资源基因科学与基因产业化研究所，国家林业和草原局/北京竹藤科学与技术重点开放实验室 北京 100102)

氮是影响植物生长的重要营养元素，水分平衡是植物生长不可或缺的基础条件。植物体内水和氮运输之间的关系错综复杂，从蒸腾驱动的土壤中的物质流动，到通过膜转运器和通道被根系吸收并运输到地上器官，涉及水分调节和氮运输的信号级联、根系构型和栓化的疏水根屏障、质外体水和氮向维管系统迁移控制等许多方面。植物对氮等大量营养元素的可利用性在很大程度上受土壤中可供水量的调节，在干旱时植物生长同时受水和氮的限制[1]。土壤氮添加能够增加叶片和叶肉细胞厚度、单位叶面积细胞间隙暴露的叶肉表面积和气孔开放状态，提高叶片氮含量、水通道蛋白(AQP)和碳酸酐酶活性，显著促进叶肉和气孔导度的增加[2]。研究显示，当大麦(Hordeumvulgare)叶片中氮含量减少时，AQP基因(HvTIPs)表达上调，表明它们参与了氨和尿素的液泡卸载[3]；氮的吸收和同化同时发生在单穗果子蔓凤梨(Guzmaniamonostachia)的单个叶片中，在叶基部通过脲酶和硝酸还原酶吸收氮和生产铵，硝酸盐和铵的高亲和转运蛋白(NRT2.5和AMT1.2)、尿素转运蛋白(DUR3)和水通道蛋白(TIP2.1)基因上调，而在叶尖中硝酸盐和铵的低亲和转运蛋白(NRT1.2和AMT3.1)、水通道蛋白(PIP1.2)基因呈高表达[4]。由此表明，在植物生长中氮转运蛋白基因和水通道蛋白基因的表达存在密切的关联性。

水和氮对作物生产力的影响超过了其他大多数环境因子[1]。在全球气候变化日趋明显的背景下，植物获取充分的氮素营养并维持细胞水分平衡，是维持植物正常生长的关键。研究表明，植物AQP能够促进尿素通过膜的扩散，如高大鹦哥凤梨(Vrieseagigantea)的VgTIP2；1除了运输水外还能运输尿素，而VgPIP1；2可能促进铵/氨的扩散[5]。竹子具有生长速度快、材质优良和适应范围广的突出特点，通过施肥可以实现竹子的高产。目前，对竹子中氮吸收、运转的基因[6-8]和水通道蛋白基因[9-10]虽然已有一些报道，但对于在氮添加处理下AQP基因的表达模式以及氮运转蛋白和AQP基因表达的转录调控尚未见报道。本研究利用前期项目组获得的转录组数据，对不同氮浓度处理下毛竹(Phyllostachysedulis)AQP基因的表达模式进行分析，针对差异表达AQP基因分析其启动子中的转录因子调控元件，进一步鉴定差异表达的转录因子(TF)，构建差异表达AQP基因和TF的共表达网络，筛选关键的基因，并用转录组数据验证，以期为后续开展基因功能提供参考依据。

1 材料与方法

1.1 实验材料

选取在人工气候室[光周期为16 h光/8 h暗、光照强度为300 μmol/(m2·s)、温度为28 ℃、湿度为60%]培养2个月左右、且生长状态一致的毛竹幼苗用于实验，将其根部清洗干净，然后放入含有不同浓度KNO3(N0=0 mmol/L，N6=6 mmol/L，N18=18 mmol/L)的改性木村B溶液中培养2周[8,11]。随后采收幼苗根系，经液氮处理后置于-80 ℃下保存。每个氮处理组包含3个生物重复。利用试剂盒从9个样本中提取总RNA，进行转录组和小RNA测序。所有小RNA和转录组数据序列已存入NCBI数据库(SRA)，项目编号为PRJNA797724 (https://www.ncbi.nlm.nih.gov/bioproject/PRJNA797724)(SRR17650178—SRR17650186)和PRJNA797734 (https://www.ncbi.nlm.nih.gov/bioproject/PRJNA797734)(SRR17635191—SRR17635199)[8]。

1.2 差异表达基因筛选和启动子元件分析

为了鉴定N0、N6和N18处理之间的差异基因，利用DESeq2_EBSeq软件对转录组数据中的差异基因(DEGs)进行筛选，筛选标准为Fold Change (|FC|)≥1.5和错误发现率(FDR)<0.05的DEGs。基于转录因子结合位点预测网站(http://plantregmap.gao-lab.org/binding_site_prediction.php)[12]，以水稻(Oryzasativa)为参考物种，筛选标准阈值p≤e-4、q≤0.05，对PeAQPs的1 000 bp启动子区域的转录因子(TF)结合元件进行筛查。

1.3 相关性分析和共表达网络构建

根据PeAQPs和TFs在氮处理条件下的表达模式，利用基迪奥在线分析平台(https://www.omicshare.com/tools/home/report/reporticabg.html)，基于皮尔森相关系数(PCC)对PeAQPs和TFs间的相关性进行分析，筛选条件为：|cor|≥0.9和P≤0.05。利用Cytoscape (v.3.7.1)软件构建TFs和PeAQPs的调控网络并进行可视化展示。

2 结果与分析

2.1 不同浓度氮处理下PeAQPs的表达

利用不同浓度氮处理的毛竹根转录组数据，共鉴定到30个差异表达的PeAQPs，包括14个PePIPs、8个PeTIPs、7个PeNIPs和1个PeSIP2-1，分别占各亚家族成员的66.7%、42.1%、50.0%和33.3%(图1A)。这些PeAQPs在不同浓度氮条件下表现出不同的表达模式，主要分为2类。第1类，随着氮浓度的提高，基因上调表达；第2类，随着氮浓度的提高，基因下调表达(图1B)。在PePIPs亚家族的14个基因中，PePIP2-3、PePIP1-6、PePIP2-2、PePIP2-4、PePIP2-13、PePIP2-14和PePIP2-15为第1类，且PePIP2-3的表达量随着氮浓度提高而上升，其余6个基因表达量为先上升后下降，但N18组仍然高于N0组；而PePIP1-2、PePIP1-3、PePIP1-4、PePIP2-6、PePIP2-8、PePIP2-11和PePIP2-12为第2类。在PeTIPs亚家族的8个差异表达基因中，PeTIP1-2、PeTIP4-2和PeTIP4-3为第1类，且PeTIP4-2的表达量随着氮浓度提高而上升，PeTIP1-2和PeTIP4-3表达量为先上升后下降；PeTIP1-1、PeTIP2-1、PeTIP2-2、PeTIP2-3和PeTIP2-4为第2类。在PeNIPs亚家族的7个差异表达基因中，PeNIP1-1和PeNIP1-2为第1类；PeNIP1-3、PeNIP1-5、PeNIP2-2、PeNIP2-4和PeNIP3-3为第2类。PeSIP2-1基因的表达模式为第2类，随着氮浓度的提高，基因持续下调表达。

注：A：基因表达数量；B：基因表达谱。图1 PeAQPs在氮处理下的表达分析Fig.1 Expression analysis of PeAQPs under nitrogen treatments

2.2 PeAQPs启动子的转录因子调控元件

为探究上述30个差异表达PeAQPs的调控机制，对基因启动子中的TF调控元件进行分析。结果表明，在这些PeAQPs上游启动子的1 000 bp序列中，除了4个基因(PeNIP1-5、PeTIP2-3、PeTIP2-4和PeTIP4-2)启动子中未鉴定到TF调控元件外，其余26个PeAQPs启动子中共鉴定出14种TF调控元件(图2)。其中Dof转录因子结合位点在PeAQPs的启动子中分布最广，在20个基因启动子中均鉴定到，包括11个PePIPs、5个PeNIPs和4个PeTIPs。其中PeNIP1-1和PeTIP4-1拥有最多的Dof结合元件，分别为12个和10个，其他基因启动子中均拥有1～9个Dof的结合位点。BBR-BPC转录因子结合位点在PeAQPs的启动子中分布最多，在12个基因(6个PePIPs、4个PeNIPs和2个PeTIPs)启动子中共鉴定出117个结合位点，其中PePIP2-6、PePIP1-6和PeNIP1-3拥有最多的BBR-BPC结合元件，分别为35个、20个和23个，其他基因启动子中均拥有1～9个BBR-BPC的结合位点(图2A)。在5个PeAQPs的启动子中鉴定出C2H2转录因子调控元件，其中PeTIP2-2启动子中含有9个该元件。另外，在这些基因启动子中还鉴定出bZIP、ERF、GATA、B3、MICK-MADS、MYB-Related、NAC、SBP、HSF、HD-ZIP和AP2等转录因子调控元件(图2B)。

注：A：高峰度TF结合元件；B：低峰度TF结合元件。图2 PeAQPs启动子中的TF结合元件数量统计Fig.2 Statistics of the number of TF binding elements in the promoters of PeAQPs

在PePIPs亚家族中，PePIP1-6的启动子鉴定出的转录因子结合元件最多，包括20个BBR-BPC元件、4个bZIP元件、4个MICK-MADS元件和2个Dof元件，以及2个HSF元件。另外，在PePIP2-6启动子中，除了发现35个BBR-BPC元件外，还有3个C2H2元件以及1个Dof元件。在PePIP1-2、PePIP2-3、PePIP2-11和PePIP2-13启动子中均只鉴定到1种TF元件，分别为4个BBR-BPC元件、3个Dof元件、1个ERF元件和4个B3元件。在其余8个PePIPs的启动子中，均存在2～4种转录因子调控元件。在PeTIPs亚家族中，PeTIP1-2的启动子中鉴定出4种TF元件，包括6个ERF元件、3个Dof元件、1个GATA元件和1个C2H2元件。在PeTIP2-1、PeTIP2-2和PeTIP4-1启动子中也鉴定到Dof元件。而BBR-BPC元件在PeTIPs启动子中分布并不广泛，只在PeTIP1-1和PeTIP2-2中鉴定到6个和8个，在PeTIP1-1中还鉴定到NAC元件和B3元件。在PeTIP4-1的启动子中只鉴定到10个Dof元件。在PeNIPs亚家族中，Dof元件和BBR-BPC元件的分布较多，在5个PeNIPs启动子存在Dof元件，4个PeNIPs启动子中存在BBR-BPC元件。其中PeNIP1-3启动子中鉴定到23个BBR-BPC元件、4个ERF元件、1个MICK-MADS元件和1个C2H2元件，而PeNIP1-1的启动子中只鉴定到12个Dof元件。PeNIP2-4中不仅含有5个Dof元件和1个BBR-BPC元件外，还有1个HD-ZIP元件。这些结果表明，PeAQPs不同亚家族成员含有的元件种类和数量都不尽相同，表明它们的表达可能受到不同转录因子的调控，表达的差异进而导致功能方面的不同。

2.3 差异表达TF的筛选和PeAQPs的共表达网络

通过对转录组数据中差异表达基因的筛选和基于各个不同数据库的注释结果，在不同浓度氮处理下的差异表达基因中共筛选出979个差异表达转录因子(DETF)，在这些DETF中，包括AP2/ERF、Homebox、MYB、NAC以及bHLH等重要的TF家族。结合PeAQPs启动子的转录因子元件分析结果进行进一步筛选分析，从中筛选出403个可能与PeAQPs有关的DETFs，它们大致可属于10个基因家族，其中83个AP2/ERF包含14个AP2、64个ERF和5个RAV亚家族成员。30个C2C2包含20个Dof和10个GATA亚家族成员，在30个B3中还有20个ARF亚家族成员。此外，这403个DETF还包含74个NAC、48个C2H2、42个bZIP、30个HD-ZIP、23个HSF、22个MYB-related、13个SBP、7个MADs-MICK和1个BBR-BPC成员。推测这些DETFs在PeAQPs转录表达过程中发挥着重要的调控作用。相关性分析结果显示，310个DETFs和30个PeAQPs之间形成3 249个共表达基因对，结合PeAQPs的启动子TF元件结果，最终筛选出由21个PeAQPs和158个DETF形成的315对共表达关系的网络。

2.3.1PePIPs共表达网络分析

在PePIPs亚家族中，除了PePIP1-2、PePIP2-2和PePIP2-15未鉴定到共表达TF外，其余11个PePIPs与10个TF家族的83个DETF形成153个共表达关系对，包括91个正向共表达关系对和62个负向共表达关系对(图3)。分析发现，除了PePIP2-11外，其余10个PePIPs均与14个Dof存在28对正向共表达和21对负向共表达关系。PePIP1-3和PePIP2-11分别与13个和10个ERF存在共表达关系，PePIP1-4与3个AP2正向共表达，与2个MYB-related负向共表达。在PePIP1-6、PePIP2-13和PePIP2-14启动子中均鉴定到bZIP结合元件，基于转录组数据，也筛选到24个bZIP家族成员与这3个基因分别存在2个、24个和23个共表达关系。PePIP2-6与11个C2H2共表达，PePIP2-8与5个SBP和6个GATA共表达，这3个转录因子的调控元件在PePIPs的其他成员启动子中均未鉴定到，推测这3个转录因子可能参与PePIP2-6和PePIP2-8的表达调控。

2.3.2PeTIPs共表达网络分析

在PeTIPs亚家族中，PeTIP1-1、PeTIP1-2、PeTIP2-1和PeTIP4-1与49个DETF形成53个共表达关系对，包括34个正向共表达关系对和15个负向共表达关系对(图4)。PeTIP1-1与1个B3和9个NAC正向共表达，与3个NAC负向共表达，而这2类TF与PeTIPs其他成员均无共表达关系。PeTIP1-2只与3个ERF正向共表达，且PeTIP1-2启动子中有6个ERF结合位点，推测这3个ERF可能参与PeAQPs的氮响应过程。PeTIP2-1与6个Dof和12个bZIP正向共表达，与3个Dof和10个bZIP负向共表达。PeTIP4-1只与6个Dof存在共表达关系，且PeTIP4-1的启动子中鉴定到10个Dof结合元件，推测Dof与PeTIP4-1的表达存在着密切关系。

图3 PePIPs与TFs的共表达分析Fig.3 Co-expression analysis of PePIPs and TFs

图4 PeTIPs与TFs的共表达分析Fig.4 Co-expression analysis of PeTIPs and TFs

2.3.3PeNIPs和PeSIPs共表达网络分析

共表达结果显示，5个PeNIPs与6个TF家族的86个成员形成107个共表达关系对，包括57个正向共表达关系对和50个负向共表达关系对(图5)。PeNIP1-1和PeNIP1-2均只与Dof存在共表达关系，其中PeNIP1-1和PeNIP1-2均与7个Dof共表达，除此之外，PeNIP1-2还与1个Dof存在共表达关系。在这2个基因启动子中，分别鉴定到12个和3个Dof结合元件。PeNIP1-3与48个TF共表达，包括22个ERF、4个MICK-MADS和22个C2H2，而在其启动子中，也均存在这3种TF的结合元件。PeNIP2-2与8个Dof和8个MYB-related共表达，其中包括与3个Dof和2个MYB-related正向共表达。PeNIP2-4与7个Dof和21个HD-ZIP共表达。另外，PeSIP2-1的启动子中只鉴定到GATA元件，且筛选获得2个与其具有正向共表达关系的GATA基因，推测GATA基因对于PeSIP2-1的表达发挥着重要作用。

图5 PeNIPs和PeSIPs与TFs的共表达分析Fig.5 Co-expression analysis of PeNIPs and PeSIPs with TFs

2.4 关键PeDofs与PeAQPs的相关性

分析发现，Dof转录因子与PeAQPs的3个亚家族成员(8个PePIPs、2个PeTIPs和4个PeNIPs)存在共表达关系，7个PeDofs与14个PeAQPs存在65对共表达关系，推测氮处理下PeDofs可能在毛竹根系水分运输过程中发挥重要的调控作用(图3—图5)。结合TF结合元件的分析结果，PePIP1-3、PePIP2-13、PePIP2-14、PeTIP4-1、PeNIP1-1和PeNIP2-4等6个基因的启动子中均具有5个及以上的Dof结合位点，推测这7个PeDofs通过调控6个PeAQPs的表达从而影响其水分和氮的吸收与转运。相关性分析和共表达结果显示(图6)，PeDof_17788、PeDof_21000、PeDof_41006和PeDof_12766与PePIP1-3和PeNIP2-4正向共表达(cor≥0.92)，与PePIP2-13、PePIP2-14、PeTIP4-1和PeNIP1-1等基因负向共表达(cor≤-0.87)；PeDof_38007、PeDof_43753和PeDof_16228则与前4个Dof转录因子相反。另外，PeDof_21128只与PePIP1-4存在负向共表达关系。因此，推测这些PeDofs通过不同类型的调控方式影响PeAQPs的表达，从而影响水分和氮的协调运输。

图6 关键PeDofs与PeAQPs的相关性热图(A)及共表达网络(B)Fig.6 Correlation heat map and co-expression network of key PeDofs and PeAQPs

3 讨论

氮是一种重要的大量营养元素，是构成植物核酸、蛋白质、辅酶因子和植物激素的重要组成部分。氮对植物的生长发育起着决定性作用，它在植物叶片伸展以及开花过程中的作用已经得到了广泛证实，在植物根系生长和发育过程中的功能也得到了深入研究。AQP是重要的跨膜通道蛋白，通过细胞和细胞器膜进行水分子的跨细胞运动，同时还参与输送甘油、尿素、金属类、二氧化碳、氮气、过氧化氢和氨等小中性分子。因此，AQP在维持细胞内碳氮水平、转运关键调控分子等方面均发挥关键作用[13]。不同浓度KNO3处理的毛竹中AQP基因表达差异变化表明，不同的AQP基因可能以不同的方式参与氮的运输，另外差异表达的TF及其与AQP基因的共表达结果表明，这些TF在受不同氮处理影响的同时，还参与AQP基因表达的调控，进而参与毛竹体内的水分和氮的运输，以维持水分平衡和氮分配。

研究表明，AQP基因的表达除受TF调控外，还受到miRNA的调节以及蛋白修饰的影响。如在小麦中tae-miR1131与tae-miR408剪切抑制TaPIP1的表达，热激蛋白TaHSP90可与TaPIP1相互作用[14]。水稻叶片中水和氮平衡的机制引起参与膜转运的蛋白质的广泛变化，包括水通道蛋白OsPIP2-6的磷酸化[15]。在高浓度氮处理下玉米根系导水率的快速变化与质膜中的PIP蛋白丰度相关，表明PIP翻译后机制在根系水动力学参数的短期调节中发挥着重要作用[16]。另外，激素信号也能调节AQP基因的表达。在大麦所有HvTIPs基因启动子区域中均含有茉莉酸甲酯(MeJA)诱导的顺式作用基序，MeJA处理导致大麦叶片中氮含量下降[3]。外源脱落酸(ABA)处理，能够增加大麦茎的导水率和HvPIP2的表达丰度[17]。在大豆GmTIP1;6基因启动子中含有ABA和MeJA顺式作用元件，其表达受ABA和MeJA的强烈诱导[18]。因此，要深入解析竹子快速生长中AQP的作用机制，还需要考虑miRNA的调节、AQP翻译后修饰以及不同激素的影响。