袁婷婷，朱成磊，杨克彬，宋新章，高志民＊

(1. 国际竹藤中心竹藤资源基因科学与基因产业化研究所，国家林业和草原局/北京市共建竹藤科学与技术重点实验室，北京 100102；2. 浙江农林大学，浙江 杭州 311300)

氮是影响植物生长发育的重要矿质营养元素，也是叶绿体、核酸、蛋白质以及很多次生代谢产物的重要组成成分[1]。大多数陆生植物以硝态氮作为主要氮源，并以主动运输的方式从外部环境吸收硝态氮，运输到其它组织器官中，供植物生长发育[2]。植物体内的硝态氮转运体( NRTs)负责从土壤中吸收硝态氮、植物体内硝态氮的运输和细胞内硝态氮再分配，如NRT1、NRT2和NRT3等[3-4]。植物在进化过程中形成2类不同的转运蛋白，即高亲和力转运蛋白/系统( HATS)和低亲和力转运蛋白/系统(LATS)，它们又分别分为组成型(cHATS和cLATS)和诱导型(iHATS和iLATS)转运系统。拟南芥(Arabidopsis thaliana)AtNRT1.1(CHL1)是高等植物中第1个被克隆的硝态氮转运蛋白基因，它的突变导致硝态氮的吸收减少[5]，之后又在水稻(Oryza sativaL.)、玉米(Zea maysL.)和小麦(Triticum aestivumL.)中获得了多种同源基因[6-8]。由于NRT1和PTR (Peptide transporter)成员序列同源性高，且能转运种类广泛的底物，它们被统一命名为NRT1/PTR FAMILY，即NPF家族[9]。

植物NPF蛋白可运输硝态氮、多肽和脱落酸(ABA)等多种底物，是一个功能广泛的家族，可分为8个亚家族[9-10]。研究表明，在拟南芥中AtNPF4.6负责硝态氮的吸收，它是一个组成型表达的低亲和转运蛋白，也可转运ABA；AtNPF7.3介导硝态氮的木质部装载[11]；AtNPF1.2的作用是负责木质部到韧皮部的转移，将硝酸盐重新分配到正在发育的叶片中[12]。水稻NPF8.9负责硝酸盐的吸收，有2个不同剪接形式的转录本(NPF8.9a和NPF8.9b)，它们过表达可增加植株地上干质量；OsNPF6.5作用于硝态氮的吸收和地上地下硝态氮的转运[13]。玉米中ZmNPF6.6是一个高亲和硝酸盐特异转运蛋白，而ZmNPF6.4是一个低亲和转运蛋白，同时具有转入和转出硝酸盐的双向转运活性[14]。NPF广泛存在于多种植物中，且成员众多，如拟南芥NPF家族有53个成员，水稻NPF家族有93个成员[13]，之后在黄瓜(Cucumis sativusL.)、毛果杨(Populus trichocarpaTorr. & Gray)、苹果(Malus domesticaBorkh.)、甘蔗(Saccharum officinarumL.)中相继鉴定了多个NPF成员[15-18]，而竹子中NPFs的相关研究尚属空白。

竹子是禾本科(Gramineae)竹亚科(Bambusoideae)植物，其中，毛竹(Phyllostachys edulis(Carrière) J.Houz.)是最具特色的竹种之一，具有生长快、材性好、笋材两用、固碳能力强等特点[19-20]。速生的特点，使得毛竹必须从外界吸收大量的氮素供其快速生长所需[21]，因此，生产中毛竹林施用氮肥非常普遍。在一定范围内(250 kg·hm-2)，随着氮肥施用量增加，毛竹冬笋、春笋和竹材的产量均呈增加的趋势[22]；然而，毛竹如何吸收利用氮素仍不清楚，其分子调节机制及对外界胁迫的响应机制等均有待深入研究。NPF家族成员为氮素的重要运输者，因此，开展毛竹NPF基因的鉴定与表达模式研究对于揭示其功能具有重要意义。本研究以毛竹为研究对象，利用已有基因组数据[23]对其中NPF家族基因成员进行全面分析，包括基因启动子、基因结构、编码蛋白理化性质、保守结构域、系统进化等，并利用毛竹转录组数据分析NPF基因的组织表达特异性，以及非生物胁迫冷和干旱处理、激素GA3和NAA处理后的表达模式，以期为深入研究毛竹中NPF基因的功能奠定基础，为探索毛竹快速生长发育期的氮素吸收利用机制提供科学依据。

1 材料与方法

1.1 毛竹NPF基因序列获取

分别从 Rice Genome Annotation Project(http://rice.plantbiology.msu.edu/)和Tair (https://www.arabidopsis.org/index.jsp)数据库下载水稻和拟南芥NPF基因的CDS和氨基酸序列，作为诱饵序列在毛竹基因组数据库BambooGDB (http://bamboo.bamboogdb.org/)中进行BlastP、BlastN (e-value=1e-10)比对，查找并获取其中NPF同源基因的候选序列。对获得的候选序列基于SMART数据库(http://smart.embl-heidelberg.de)和PFAM数据库(http://pfam.xfam.org)逐条进行蛋白保守结构域的鉴定和分析，仅保留具有编码完整保守结构域的序列，并进行基因命名(PeNPFs)。

1.2 基于NPF同源序列的系统进化分析

为分析毛竹NPF成员的进化关系，预测其潜在的功能，使用MEGA 7.0软件的Clustal W功能对毛竹、水稻和拟南芥的NPF蛋白序列进行同源比对分析，使用邻接法(Neighbor-joining)构建系统进化树，bootstrap试验重复1 000次[24]，其它参数设置为默认值。

1.3 毛竹NPF基因结构及其编码蛋白序列分析

使用Gene Structure Display Server 2.0 (http://gsds.cbi.pku.edu.cn/)在线工具绘制PeNPFs的基因结构图，利用Plant CARE (http://bioinformatics.psb.ugent.be/webtools/plantcare/html/)在线平台分析其启动子序列所含作用元件，用ProtParam (http://web.expasy.org/prot-param/)对毛竹NPF基因编码蛋白质的基本理化特性进行分析，使用MEME Version 5.1.0 (http://meme-suite.org/tools/meme)分析毛竹NPF蛋白的保守基序，使用Plant-mPLoc算法(http://www.csbio.sjtu.edu.cn)预测PeNPFs的亚细胞位置。

1.4 组织表达特异性分析

依据本实验室前期毛竹7个不同组织(叶、早花期花序、盛花期花序、根、鞭、20 cm笋和50 cm笋)的转录组数据[23]，用不同组织中NPF基因的RPKM (Reads Per Kilo-bases Per Million Reads)值表示基因的表达丰度，为有效地可视化展示基因表达量，对RPKM的值取以2为底的对数即Log2(RPKM)作为基因的表达量[25]。采用TBtools的热图工具[26]制作热图展示基因的表达情况，颜色由蓝色到红色逐步变化表示基因的表达丰度数值递增。

1.5 获得不同胁迫与激素处理后的转录组数据

根据文献查找，从NCBI数据库(https://www.ncbi.nlm.nih.gov/sra/)下载获取毛竹非生物胁迫和激素处理的转录组数据。冷和干旱处理的转录组数据源自2年生毛竹，5株的离体枝条分别进行冷(0℃和无光条件)和干旱(20℃和相对湿度50%条件)处理2 h后，每个处理分别选取10片未展叶进行转录组测序，其RNA-seq数据登录号为SRS-1759772[27]；GA3处理的转录组数据均源自毛竹实生苗，1个月的毛竹实生苗分别用水和100 μmol·L-1的GA3喷洒处理后4 h取整株进行转录组测序，其对照和处理的RNA-seq数据登录号分别为SRS256-6465～SRS2566467、SRS2566468～SRS2566470[28]；NAA处理的转录组数据均源自1个月的毛竹实生苗，用5 μmol·L-1的NAA喷洒处理，间隔1 h处理1次，4 h后取根系进行转录组测序，其对照和处理的RNA-seq数据登录号分别为SRS2294012～SRS2294014、SRS2294015～SRS2294017[29]；上述处理取样时均采用3次生物学重复。从这些转录组数据的基因注释中比对查找，筛选出对应毛竹NPF基因的FPKM值，用于表达分析。

1.6 不同处理后基因的表达模式分析

利用获取的冷、干旱、GA3和NAA处理毛竹及其对照中NPF基因的FPKM值，比较分析处理前后的基因表达变化情况。为有效地可视化展示基因表达量，用FPKM+1的值取以2为底的对数即Log2(FPKM+1)作为基因的表达量[30]。对所有表达量数据用TBtools软件绘制表达谱热图，分析不同基因的表达情况，同时基于基因FPKM值对具有显著差异(p＜ 0.05)表达的基因绘制点线图。

2 结果分析

2.1 毛竹NPF家族基因成员鉴定与系统进化分析

植物NPF蛋白通常含有PTR2保守结构域，通过比对分析，在毛竹基因组中共鉴定到具有编码完整PTR2的候选基因27个，其编码蛋白均含有保守的跨膜结构域，符合MFS (Major Facilitator Superfamily)超家族的NPF家族的特征[31]。根据候选基因在毛竹Scaffold中的位置以及同源基因的名称，依次命名为：PeNPF1.1～PeNPF8.8。

为确定PeNPFs与其它物种中已鉴定的NPF蛋白之间的进化关系，预测其潜在功能，对毛竹、水稻和拟南芥NPF的氨基酸序列进行了多重序列比对，并利用MEGA7.0软件构建了系统进化树。结果表明：所有毛竹、水稻和拟南芥的NPFs被聚类到8个亚家族(S1～S8)，PeNPFs成员分布在除S3之外的其它7个亚家族中，每个亚家族中成员数量不等，其中，亚家族S1、S4、S6分别只有1个成员，S8中的成员最多为8个，其余亚家族分别有3～7个成员不等(图1)。

图1 基于毛竹、水稻和拟南芥NPF氨基酸序列构建的系统进化树Fig.1 Phylogenetic tree constructed based on amino acid sequence of NPFs from moso bamboo, rice and Arabidopsis

进一步分析发现，PeNPFs半数成员被聚类在S7 (7个)和S8 (8个)亚家族，而双子叶植物拟南芥在这2个亚家族的成员分别只有3个和5个，在S5亚家族中，水稻和拟南芥NPF成员分别达到17个和16个，但仅包含6个毛竹成员。虽然NPF家族成员在功能上具有相似性，但不同亚家族成员在功能上也存在一定的差异，如AtNPF1.1和AtNPF1.2负责硝酸盐向嫩叶的分配[32]，而OsNPF 2.4主要在表皮、木质部薄壁细胞和韧皮部的伴胞中表达[33]。另外，所有PeNPFs均优先与水稻NPF聚类在一起，说明在功能上PeNPFs或许与水稻等单子叶植物的同源基因相似。

2.2 PeNPFs基因结构及启动子分析

基因结构分析表明，PeNPFs均有内含子，数量为2～5个，单个内含子在基因组上的跨度为33 bp(PeNPF7.6的第5个内含子)～3 489 bp (PeNPF2.2的第2个内含子)，而单个外显子在基因组上的跨度为36 bp (PeNPF7.3的第1个外显子)～1 601 bp(PeNPF7.3第3个外显子) (图2)。PeNPFs各个成员之间的内含子和外显子的数量以及所在的位置存在一定的差异，这可能是影响其转录表达的主要原因之一。

图2 PeNPFs基因结构Fig.2 Gene structure of PeNPFs

为更好地了解PeNPFs对外界环境的响应和受内在调控因子调节的情况，对PeNPFs的起始密码子ATG上游2 000 bp的序列所含顺式作用元件和应答元件进行分析。结果表明：所有序列中均含有启动子必需的核心元件TATA-box和CAAT-box。27个PeNPFs的启动子序列中均包含多种参与基因表达调控的顺式作用元件和外界刺激的应答元件，如低温响应元件LTR，干旱响应元件MBS和TCrich repeats，赤霉素响应元件TATC-box、P-box和GARE-motif，生长素响应元件TGA-element、AuxRR-core和AuxRE，茉莉酸甲酯响应元件CGTCA-motif和TGACG-motif，脱落酸响应元件ABRE等(图3)。由此表明，PeNPFs的转录表达可能会受到非生物胁迫和激素的影响。

图3 毛竹PeNPFs基因启动子响应元件种类统计Fig.3 Statistics of response elements in the promoter of PeNPFs

2.3 PeNPFs的保守基序分析

保守基序分析表明：在毛竹PeNPFs序列中共获得10个保守基序，分别命名为motif 1～motif 10，其中，motif 1、motif 2和motif 4是27个PeNPFs共有的高度保守基序，而其它基序则在部分序列中缺失。27个PeNPFs中有15个成员含有10个motif，其它12个成员缺失1～3个motif，如PeNPF2.1、PeNPF4.1和PeNPF8.8缺少motif 3，PeNPF8.2缺少motif 5，PeNPF2.1缺少motif 6，PeNPF2.1、PeNPF7.5、PeNPF7.7、PeNPF8.2、PeNPF8.3和PeNPF8.8缺 少motif 7，PeNPF5.2、PeNPF5.5缺少motif 8，PeNPF7.6缺少motif 9，PeNPF5.1、PeNPF5.2、PeNPF5.3、PeNPF5.5缺 少motif 10(图4A)。

另外，各motif在氨基酸组成的偏好性方面有一定相似性，如motif 1、motif 2和motif 4高度保守，均富含疏水性氨基酸，但各motif在氨基酸组成的偏好性方面也存在一定的差异，如motif 1富含亮氨酸，motif 2富含丝氨酸，motif 4富含甘氨酸(图4B)。

2.4 编码蛋白理化性质和亚细胞定位分析

PeNPFs编码氨基酸序列最长的为761 aa(PeNPF7.4)，最短的为452 aa (PeNPF8.8)，蛋白分子量介于49.56～82.08 kDa之间，理论等电点介于5.17～9.85之间，大部分是中性或碱性蛋白。亲疏水性分析显示，PeNPFs的亲水性值(GRAVY)为0.063～0.607，推测它们属于亲水性蛋白。亚细胞定位预测显示，所有PeNPFs均定位于细胞膜或液泡中，其中，有8个成员在细胞膜和液泡中均有定位(表1)。

表1 PeNPFs编码蛋白的基本特征Table 1 Basic characteristics of the proteins encoded by PeNPFs

2.5 PeNPFs的组织表达模式分析

利用毛竹7个不同组织的转录组测序(RNASeq)表达谱数据，对PeNPFs的组织表达特异性进行分析。结果(图5)表明：在不同组织中，PeNPFs表达模式存在一定差异，其中，PeNPF5.3和PeNPF7.6在7个组织中均检测到表达，且在7个组织中的表达量变化不大，呈组成型表达模式，它们可能与毛竹整个生长发育过程密切相关；而其它PeNPFs成员至少在1个组织中检测到表达，但各成员的表达丰度差异较大，如在叶片中PeNPF2.3、PeNPF7.5、PeNPF7.7和PeNPF8.7的表达丰度较高；PeNPF-2.1、PeNPF5.6、PeNPF6.1、PeNPF7.1、PeNPF7.2和PeNPF8.2在早花期花序中表达丰度较高；PeNPF-1.1、PeNPF5.2、PeNPF5.6、PeNPF6.1和PeNPF8.1在盛花期花序中表达丰度较高。一些PeNPFs成员表现出组织特异性，如PeNPF4.1主要在竹鞭中表达；PeNPF5.1和PeNPF7.4主要在根中表达；PeNPF8.8主要在20 cm和50 cm高度笋中表达。PeNPFs在叶、花序、鞭和根等不同组织中的表达差异，表明它们在这些组织的生长发育中发挥着不同的功能。

图5 PeNPFs在毛竹7个组织中的表达分析Fig.5 Expression analysis of PeNPFs in seven tissues of moso bamboo

2.6 毛竹PeNPFs响应非生物胁迫和激素胁迫的表达分析

启动子调控元件分析表明：PeNPFs的表达可能会受到逆境胁迫的影响，因此，本文利用已发表的转录组数据，对冷和干旱处理下毛竹叶片中PeNPFs的表达进行分析。结果表明：在冷和干旱处理的毛竹叶片转录组数据中获得了18个PeNPFs的FPKM值，进一步分析发现它们的表达存在着一定的差异。与处理前(CK)相比，冷处理2 h后共有9个PeNPFs呈下调表达趋势，9个PeNPFs呈现上调表达趋势(图6A)；其中，下调表达的基因中有6个(PeNPF2.3、PeNPF5.3、PeNPF5.4、PeNPF7.2、PeNPF7.3和PeNPF8.7)呈显著下调(p＜ 0.05)，而在上调表达的基因中仅有3个基因(PeNPF5.2、PeNPF7.1和PeNPF8.1)表现为显著上调(p＜ 0.05) (图6B)。在干旱处理2 h后共有15个基因呈现下调表达趋势，仅有3个基因呈现上调表达趋势(图6C)；其中下调表达的基因中有6个(PeNPF5.3、PeNPF7.2、PeNPF7.3、PeNPF8.1、PeNPF8.7和PeNPF8.8)呈现显著下调(p＜ 0.05)，而上调基因的表达差异均不显著(图6D)。综合分析发现，4个基因(PeNPF5.3、PeNPF7.2、PeNPF7.3和PeNPF8.7)在冷和干旱处理后均表现出显著下调，这与其启动子序列中存在与低温和干旱响应调控元件(图2)相一致，推测它们可能在低温和干旱响应中发挥着重要作用。

图6 冷和干旱处理下PeNPFs的表达分析Fig.6 Expression analysis of PeNPFs in moso bamboo under cold and drought treatments

另外，NPF基因的表达还受到激素的影响[34]，本文对GA3和NAA处理毛竹幼苗的转录组数据中PeNPFs基因的表达模式进行分析。结果表明：在转录组数据中获得了18个PeNPFs的FPKM值，分析发现它们的表达在处理前后存在一定的差异。与处理前(CK)相比，18个PeNPFs在100 μmol·L-1GA3处理4 h后，各有9个基因呈现上调或下调表达趋势(图7A)，其中，上调表达的5个基因(PeNPF4.1、PeNPF5.3、PeNPF5.4、PeNPF7.1和PeNPF7.6)，下调表达的3个基因(PeNPF2.1、PeNPF5.1和PeNPF7.3)的表达差异显著(p＜ 0.05)(图7B)；在5 μmol·L-1NAA处理4 h后，7个基因呈上调表达趋势，11个基因呈现下调表达趋势(图7C)，其中，PeNPF7.1和PeNPF7.3表现为显著上调，而PeNPF2.2、PeNPF5.2和PeNPF7.6等基因表现为显著下调(p＜ 0.05) (图7D)。对比分析发现，PeNPF7.3和PeNPF7.6在2种激素处理后呈相反的表达变化，而PeNPF7.1的表达趋势相同，推测这些基因在响应激素方面可能发挥着不同的功能。

图7 GA3和NAA处理下PeNPFs的表达分析Fig.7 Expression analysis of PeNPFs in moso bamboo treated with GA3 and NAA

3 讨论

NPF蛋白是MFS超家族中成员数量最多的家族之一，这些蛋白质广泛分布于真核生物和原核生物中，它对于植物、动物、酵母和细菌的氮代谢具有重要作用。本研究从毛竹基因组数据库中，鉴定到27个NPF成员，这些成员均含有PTR2保守结构域及跨膜结构域，这可能是NPF家族成员具有吸收和转运硝态氮的基础。同一个亚家族内的部分NPF成员间的基因结构以及理化性质和亚细胞定位高度相似，如PeNPF8.6和PeNPF8.7的序列的相似性为86.83%，氨基酸长度相差3个，相对分子量、亚细胞位置、基因结构、所含保守基序以及在7个组织中的表达模式相似。另外，部分NPF成员的基因结构、理化性质和亚细胞定位存在较大差异，如毛竹PeNPF5.1和PeNPF5.4序列的同源性为31.19%，氨基酸长度相差136个，二者的分子量分别为56.96 kDa和71.33 kDa，且二者在7个组织中的表达模式相反，这可能是由于非保守区的差异造成的，同时各成员序列非保守区的差异可能导致各个成员的功能具有多样性[35]。

从毛竹基因组数据库中，鉴定了27个PeNPFs成员，分别属于7个亚家族，但没有鉴定到S3亚家族成员。系统进化分析发现，在其它物种(如水稻)中存在S3亚家族成员，推测本研究中没有鉴定到S3亚家族成员的原因可能是所用的是毛竹基因组草图，毛竹可能存在S3亚家族成员，作者将进一步用最新的毛竹基因组进行分析。同时，通过物种进化关系分析发现，聚类在同一个亚家族的成员可能具有相似的功能。在水稻和拟南芥中，S1、S2、S5、S6和S7亚家族成员的主要功能为转运硝态氮，如OsNPF2.2负责根到地上部分硝态氮的运输[36]，小麦TaNPF6.3负责硝态氮的吸收并调控非根系组织的硝态氮的分配[37]，AtNPF7.2和AtNPF 7.3分别负责硝态氮在木质部的卸载和装载[38]，推测毛竹在这5个分枝的成员可能也具有硝态氮的转运功能。另外，研究表明，AtNPF2.4和AtNPF2.5负责根中Cl-的转运[39]，AtNPF5.2在胁迫条件下可以转运二肽和三肽，AtNPF7.3还具有转运K+的功能[40]，所以，也不排除相同亚家族中毛竹PeNPFs也具有一些其它功能。基因在不同组织的表达模式是基因功能的重要体现，转录组表达谱表明，PeNPFs在毛竹各组织的功能存在差异，多数基因在叶片和花序中的表达量较高，在根和笋中的表达量相对较低。3个毛竹NPF基因具有组织特异性，PeNPF1.1只在盛花期花序中表达，PeNPF5.1只在根中表达，PeNPF8.8只在笋中表达，表明它们可能分别在盛花期花序、根以及笋中发挥着重要作用；而它们分别属于不同的亚家族，说明各亚家族之间功能有一定的差异。拟南芥S8是小肽转运体，AtNPF8.1和AtNPF8.2能够转运多肽，同时AtNPF8.3在开花和种子发育过程中将小肽转出液泡，而毛竹S8成员中多个成员在花序发育时期的表达量较高，说明它们也可能参与毛竹开花发育时期小肽的转运。

真核生物基因的表达受到多层次调控，转录因子能够通过与基因的启动子特定区域结合从而调控基因的转录表达，因此，对基因的启动子所含的调控元件进行分析对于基因功能的预测具有重要意义。本研究发现，在PeNPFs中发现大量与逆境胁迫及激素响应相关的调控元件，这些元件已经被证实在植物逆境胁迫[41]及植物生长与果实成熟等方面均发挥着重要作用。当植物受到生物或非生物胁迫时，NPF基因通过上调或下调表达发挥其功能。冷处理后，18个PeNPFs中有9个PeNPFs的表达量发生了显著变化(上调或下调)，说明毛竹在一定程度通过PeNPFs表达量变化来响应冷处理；干旱处理后，18个PeNPFs中有15个基因的表达量下调，其中，6个与处理前差异显著，这与干旱处理能显著诱导小麦NPF基因表达量下调的结果一致[34]。在GA3处理4 h后18个PeNPFs均发生变化，其中，PeNPF4.1、PeNPF5.3、PeNPF5.4、PeNPF7.1和PeNPF7.6的表达量显著上升，推测它们可能受GA3的诱导并通过表达量变化进行响应，这与施用GA3能够增加ZmNRT2.1/2.2的表达水平相似[42]；NAA处理后，毛竹有5个基因发生显著上调或下调，推测毛竹S7亚家族可能参与毛竹响应NAA过程。毛竹中PeNPFs不同成员对GA3和NAA处理呈现不同程度、甚至相反的表达变化，表明各成员功能的差异，这与水稻中OsPTRs在GA3和NAA处理后的表达类似[4]。PeNPFs在毛竹中的具体功能和作用机制需要进一步实验验证。

4 结论

本研究从毛竹基因组中共鉴定出27个PeNPFs基因，所有成员编码的蛋白均包含完整的PTR2保守结构域，根据其蛋白序列差异分为7个亚家族，各亚家族成员之间的分子特征和组织表达特异性均存在着一定的差异，部分成员的表达表现为组织特异性或组成型表达。在PeNPFs启动子序列中鉴定出多种与非生物胁迫和激素响应相关的调控元件，冷、干旱以及GA3和NAA处理后毛竹的转录组表达谱证明PeNPFs参与了不同程度的非生物胁迫和激素的应答，研究结果为揭示PeNPFs的生物学功能提供了重要依据。但对于PeNPFs的研究尚在起步阶段，其在毛竹中的具体功能尚需展开大量的实验进行验证。