基于硝态氮或铵态氮条件下杨树根尖miRNAs 特征分析

阚东旭，逯岩，吴江婷，陈昕，石文广，周婧＊

(1.林木遗传育种国家重点实验室，国家林业和草原局森林培育重点实验室，中国林业科学研究院林业研究所，北京 100091；2.林木遗传育种国家重点实验室，东北林业大学，黑龙江 哈尔滨 150040)

杨树（PopulusL.）作为一种快速生长的木本植物，是我国最主要的速生丰产林树种之一[1]。杨树具有速生、易繁殖、适应性强和生产力高等特性，在木材加工、碳汇造林、制浆造纸和生物燃料等方面发挥着重要作用[2-4]。目前，我国许多杨树人工林种植园生长在氮素贫瘠的地区[5]，易导致杨树人工林生产力低下，原材料利用率不高，产品质量差等问题。因此，开展杨树根系响应氮素分子调控机制的研究，可以实现杨树人工林氮素高效利用的经济效益，意义重大。

氮素是植物生长发育所必需的大量矿质营养元素之一，是蛋白质、核酸、叶绿素等重要物质的组成成分[6]。植物根系是吸收硝态氮或铵态氮的主要器官，良好的根系形态对植物氮素高效吸收利用至关重要[7]。有研究表明，植物可通过对外界不同氮形态响应做出反应，从而改变根系形态并逐渐做出相应的适应[8-12]。Wei等研究表明，杨树根系形态对低硝态氮处理响应的初期表现是加快增加主根长度，随后加速侧根增殖，在低硝态氮处理2～4 d后，根系生物活性显著增加，其根长和侧根数显著高于对照[5]。类似的动态响应在玉米(Zea maysL.)中也有报道[13]。上述研究结果表明：植物根系有能够根据不同氮形态供应而改变其形态的能力[14]。

除植物根系形态发生变化外，分子水平调节在植物响应不同氮形态吸收同化过程中也发挥着重要作用，特别是来自具有调控功能的miRNAs 的调节[15-16]。miRNAs 是一类广泛存在于真核生物中，长约18～24个碱基(nt)的高度保守小分子RNAs[17]。利用小RNA 高通量测序技术并结合生物信息学分析，前人研究中已鉴定出大量参与氮素吸收同化响应相关的miRNAs[16,18-19]，其中，部分miRNAs 也参与了植物根生长和发育过程[16,20-21]。如在拟南芥（Arabidopsis thaliana(L.) Heynh.）根中柱鞘细胞中，5 mmol·L−1硝态氮可以抑制miR167 的表达，从而提高靶基因auxin response factor 8(ARF8)的表达水平，促进拟南芥侧根起始和随后的萌发[20,22]。在水稻（Oryza sativaL.）中，miR396 的靶基因growth-regulating factor 4(GRF4) 能够促使水稻根系吸收铵态氮从而改变水稻根系结构[23]；而miRNA 调控机制在不同氮形态处理杨树根尖上存在怎样的差异表达模式尚不清楚，值得深入研究。

综上所述，本研究以灰杨(Populus×canescens)根尖为试材，利用小RNAs 高通量测序技术，分析根尖差异表达的miRNAs，并结合降解组和转录组测序结果，分析miRNAs 靶基因的差异表达模式，阐明杨树根尖对硝态氮或铵态氮处理条件下miRNAs及其靶基因的表达调控机制。本研究成果可为培育氮素吸收利用效率高的林木奠定理论基础，具有理论与实践意义。

1 材料和方法

1.1 试验材料

将灰杨组培苗在组培室 (白天/黑夜温度25℃/18℃，相对湿度50%～60%，光照长度14 h，光照强度150 μmol·m−2·s−1) 中培养4 周。随后，将生根的灰杨植株种植在10 L 的花盆里进行砂培。每株植株隔天浇50 mL Long Ashton (LA)营养液(0.5 mmol·L−1NH4NO3、0.5 mmol·L−1KCl、0.9 mmol·L−1CaCl2、0.3 mmol·L−1MgSO4、0.6 mmol·L−1KH2PO4、42 μmol·L−1K2HPO4、10 μmol·L−1Fe-EDTA、2 μmol·L−1MnSO4、10 μmol·L−1H3BO3、7 μmol·L−1Na2MoO4、0.05 μmol·L−1CoSO4、0.2 μmol·L−1ZnSO4和0.2 μmol·L−1CuSO4,pH 5.5)。砂培植株在温室(与组培室气候条件相同)中培养14 d 后，具有相似高度(株高约30 cm)的植株被选择并分成2 组(每组24株植物)。水培1～2 周后，将这2 组植株移植到分别添加等量氮供应的0.5 mmol·L−1NaNO3(硝态氮)和0.5 mmol·L−1NH4Cl (铵态氮) 代替0.5 mmol·L−1NH4NO3的改良LA 营养液中进行水培，培养时间为10 d。在处理10 d后对其进行根系形态特征分析，然后选取根尖40 mm 进行收获。所有收获样品迅速置于液氮中，放于−80℃冰箱备用。为了获得足够的材料进行进一步分析，将8株植物的样本均匀混合作为一个重复，每种处理水平3个生物学重复。

1.2 试验方法

1.2.1 根系形态特征分析 试验苗进行硝态氮或铵态氮处理后，对灰杨根系形态构型进行监测，包括测量根系长度，二级可见侧根发生位置等指标。

1.2.2 不同处理水平根尖小RNA 测序 不同处理水平根尖小RNAs 测序按照Illumina 公司提供的标准步骤执行。简单来说，将2种不同氮形态处理下的样品，每个处理3个重复分别利用RNA 提取试剂盒提取总RNAs (TRK1001,联川生物技术公司,杭州,中国)，将提取的总RNAs 质量检测后，利用miRNA 特有的5′端磷酸基团和3′端羟基基团属性，将一个腺苷化单链DNA 3′接头和5′接头相继连接到小RNAs 上，通过与3′端互补的RT 引物进行反转录反应，对反转录产生的cDNA序列进行PCR 扩增。随后通过6% polyacrylamideTris-borate-EDTA 切胶回收的方式对140～160 bp 长度范围的PCR 产物进行胶回收，最后利用杭州联川生物技术公司Illumina HiSeqTM2000 高通量测序技术进行小RNA 测序分析。

1.2.3 miRNAs 测序结果分析 将测得的RNAs 序列使用联川生物技术公司开发的软件进行分析。首先，去除低质量的序列，保留干净的唯一的序列进行分析，将剩余序列比对RFam 数据库和重复序列数据库用以去除rRNA、tRNA、snRNA 和snoRNA等，过滤后的数据使用已有的灰杨基因组数据库(http://aspendb.uga.edu/index.php/databases/spta-717-genome)、miRNAs 数据库(http://www.mirbase.org/)进行比对，筛选出已知miRNAs。针对新的miRNAs 序列，利用miRNA 前体的标志性发夹结构，通过Mfold 软件(http://rna.tbi.univie.ac.at/cgibin/RNAfold.cgi)分析RNA 二级结构，鉴定新的miRNAs。

1.2.4 miRNAs 靶基因分析 把等量的每个处理水平的3个重复样品合并，共2个处理水平的合并样品利用联川生物技术公司Illumina HiSeq2500 技术平台进行降解组测序，采用CleaveLand 软件(3.0.1版本)预测miRNAs 的靶基因。为了进一步挖掘靶基因的表达水平，6个cDNA 文库被构建，利用转录组数据Hiseq 测序平台分析相关miRNAs 靶基因的表达模式，同时结合生物信息学Kyoto Encyclopedia of Genes and Genomes (KEGG)软件和MapMan 软件，分析差异表达miRNAs 靶基因的生物学功能。

1.2.5 荧光定量PCR 引物设计及合成 根据灰杨数据库，利用Primer Premier 3.0 软件设计引物并由上海生工生物工程有限公司合成相关miRNAs 和靶基因的引物。5.8s rRNA 和Actin 分别作为miRNAs和靶基因的内参基因。序列见表1。

表1 RT-qPCR 引物序列Table 1 Primers used for RT-qPCR

1.2.6 统计分析 使用Statgraphics 软件(STN,St Louis,MO,USA)进行数据统计分析。在进行统计分析之前，先对数据进行正态性检验。所有数据均采用单因素方差分析(ANOVA)，以P<0.05 作为统计意义上的显著水平。

利用transcripts per million(TPM)对miRNAs表达水平进行定量。使用miRNA 在铵态氮处理下的TPM 除以硝态氮处理下的TPM 来计算miRNA的差异倍数(FC)。差异表达miRNAs 筛选阈值为P<0.05。

采用Fragments Per Kilobase of exon per Million mapped (FPKM)对靶基因mRNAs 表达水平进行定量。基于FPKM 值，使用Ballgown package 计算mRNAs 的差异表达水平。使用mRNAs 在铵态氮处理下的FPKM 除以硝态氮处理下FPKM 来计算这个基因的差异倍数(FC)。差异表达mRNAs 筛选阈值为log2(FC) ≥1 或≤ −1，P<0.05。

对于RT-qPCR 的测定，将qPCR 得到的Ct 值进行归一化，计算miRNAs 及其靶基因的相对表达量。

2 结果与分析

2.1 根系形态特征分析

由图1 可知：在不同氮形态处理条件下灰杨根的形态特征具有显著差别。测量表明，硝态氮处理下的根长为27.5 cm，比铵态氮处理的14.3 cm 几乎长1 倍(p<0.001)，但2种处理水平下其二级侧根发生位置比较一致。这说明不同氮形态处理对根系生长发育的影响不同。选取二级侧根形成前的根尖部分(40 mm)进行后续分子实验。

图1 不同氮形态处理灰杨根系的表型Fig.1 Morphological parameters of P.×canescens roots with different nitrogen forms

2.2 miRNAs 深度测序数据分析

基于上述根尖形态特征分析的结果，进行高通量小RNAs 测序分析。测序结果显示：2个处理水平的测序分别得到的原始序列读数为5 752 458 和7 800 546，进行数据处理和杂质的过滤，去掉低质量序列和3′、5′缺失序列等，分别获得1 996 697和2 608 788 条干净序列(表2)。长度在18～25 nt的序列被保留下来做进一步研究。

表2 小RNA 文库测序数据统计Table 2 Distribution of small RNAs in different categories

2.3 已知和新的miRNAs 鉴定

在不同氮形态处理下，共鉴定灰杨523个已知的miRNAs 和42个新的miRNAs (表3)。研究表明，植物miRNAs 长度大多在21 nt 或22 nt，统计本实验2个处理水平miRNAs 文库中所有18～25 nt 的miRNAs 长度，21-nt miRNA 出现频率最高 (图2)，这说明该测序数据可靠。

图2 灰杨根尖已知和新miRNAs 长度分布Fig.2 Lengths of known and novel miRNAs in root tips of P.×canescens.

表3 基于高通量测序技术鉴定灰杨根尖已知和新的miRNAs 数量Table 3 The number of identified known and novel miRNAs in root tips of P.×canescens based on high-throughput sequencing

2.4 不同氮形态处理下灰杨根尖miRNAs 差异表达

为了解不同氮形态处理下，灰杨根尖对miRNAs的响应，对筛选出来的灰杨中已知miRNAs 和新miRNAs 的表达量进行计算。分析表明：灰杨根尖共有96个miRNAs (p<0.05)在不同氮形态处理下展示出不同的表达模式，其中，包括33个已知miRNAs 家族和9个新的miRNAs。相比较于硝态氮处理，铵态氮处理有44个上调表达的miRNAs和52个下调表达的miRNAs (图3a)，其中29个显著差异表达的miRNAs (p<0.001)见表4。

表4 29个显著差异表达miRNAsTable 4 twenty-nine significantly differentially expressed miRNAs

同时，随机挑选了10个显著差异表达miRNAs，利用RT-qPCR 进一步证实了测序结果的可靠性(图3b)。虽然miRNAs 表达水平的差异与测序得到的差异并不一致，但上调或下调表达的趋势是相似的。2种技术在比率方面的差异是由2种技术的算法和灵敏度的本质不同造成的。总体而言，本研究测序结果的可靠性较高。

图3 不同氮形态处理条件下灰杨根尖差异表达的miRNAsFig.3 Significantly differentially expressed miRNAs in root tips of P.×canescens under different nitrogen forms treatments

2.5 不同氮形态处理下灰杨根尖靶基因的功能

为了进一步了解miRNAs 的作用，分别使用等量的总RNA 样本对硝态氮和铵态氮2个混合降解组池进行测序，以识别miRNAs 靶基因。结果表明：分别有67.13%和65.37%去除冗余的reads 可以比对到灰杨数据库中。从降解组中共鉴定出2729个靶基因位点，其中，2104个靶基因显著差异表达。为了进一步了解上述靶基因的功能，对识别出的2104个显著差异表达靶基因进行KEGG 通路分析(图4)；同时MapMan 分析进一步揭示了这些差异靶基因在氮代谢过程中以及植物生长发育过程中参与的各种生物学过程，如ptc-miR166i-p5 有4个靶基因potri.015g017500.1/2/3/4均属于NADH依赖的谷氨酸合成酶(NADH-dependent glutamate synthase family) 蛋白参与氮代谢过程。gra-MIR8723bp3_2ss6TC21TC 有3个靶基因Potri.003G111500.1/2/3均作为硝态氮转运体基因(nitrate transmembrane transporter1.1/NRT1.1) 在硝态氮运输中起到关键作用。ptc-MIR6462a-p5_1ss14TC 的靶基因Potri.005 G079200.1、gma-miR6300_1ss5TG和gma-miR6300_R+1_1ss5TG的靶基因Potri.001G300900.1/2均 编码氨基酸运输或代谢相关蛋白，也对氮代谢有所响应 (表5)；同时，与植物生长发育相关的一些miRNA的靶基因也相继被发现。如miR164 家族的多个靶基因均属于NAC基因家族成员，其在植物生长发育过程中发挥重要功能 (表5)。

表5 灰杨根尖氮代谢以及生长发育相关差异表达靶基因Table 5 Significantly differentially expressed target genes related to nitrogen metabolism and growth and development of P.×canescens root tips

图4 KEGG 通路分析鉴定差异表达靶基因Fig.4 KEGG pathway analysis of significantly differentially expressed target genes

2.6 不同氮形态处理下灰杨根尖miRNAs-靶基因调控网络

为了研究不同氮形态处理下，杨树根尖差异表达miRNAs 与靶基因调控网络，利用上述结果进行联合分析。表6 表明：有23个差异表达靶基因(log2(FC)≥1或≤−1，P<0.05) 隶属于5个miRNA家族和一个新的miRNA (P<0.05)被鉴定，且miRNA-靶基因呈现负相关调控关系，表明这些靶基因可能是通过miRNAs 转录抑制作用而被降解。

表6 不同氮形态处理条件下灰杨根尖差异表达miRNA-靶基因Table 6 Significantly differentially expressed miRNAs-target pairs in root tips of P.×canescens under different nitrogen forms treatments

通过RT-qPCR 实验验证了其中7个miRNA-靶基因的表达模式，结果显示：差异表达miRNAs及其靶基因RT-qPCR 分析与高通量测序产生的表达模式相似，且表达模式呈负相关关系 (图5)。

图5 不同氮形态处理下miRNAs-靶基因对差异表达分析Fig.5 Validation of significantly differentially expressed miRNAs and their targets under different nitrogen forms treatments by sRNA-seq and RT-qPCR

3 讨论

研究表明，基因的表达受miRNAs 调节，不同氮形态可能导致miRNAs 及其调控的靶基因差异表达，从而影响植物根尖生长和发育[24]。本研究利用高通量测序平台研究不同氮形态处理下，杨树根尖差异miRNAs 和靶基因的表达模式。

在不同氮形态处理10 d 后，收获杨树根尖40 mm 区段材料，经过小RNA 高通量测序分析，获得96个显著差异表达miRNAs，其中，vvi-MIR399dp3_1ss13GA 上调表达倍数最高。前人研究表明，miR399 参与调控水稻多种营养饥饿反应[25]。在本研究中，miR399 靶基因为泛素化连接酶，且与miR399 呈负相关关系。有研究证明，泛素分子主要通过泛素活化酶、泛素结合酶和泛素化连接酶将靶蛋白泛素化，泛素化的蛋白最后被26S 蛋白酶体识别和降解[26]，该过程对植物营养缺乏等非生物胁迫发挥着重要作用[27]。因此，可以推测miR399 及其靶基因泛素化连接酶在不同氮形态胁迫响应中发挥着重要的调控作用。在本试验条件下，部分差异表达miRNAs 与杨树根尖生长发育相关，这与在水稻中的研究结果相似[15]。如miR171 家族成员，其靶基因之一是DELLA protein RHT-1。DELLA蛋白家族是GA 信号途径中的负调控因子，可以抑制GA 途径的基因表达从而抑制植物生长[28]。另有研究表明，DELLA-GRF4介导植物生长与氮代谢的协同调控机制，DELLA的积累不仅导致水稻生长矮化，而且降低了氮素利用效率[23]。在本研究中，相比较于硝态氮处理，铵态氮处理诱导miR171 靶基因DELLA表达，从而导致杨树根尖发育受到抑制，根长变短。该研究说明杨树根尖部分差异表达miRNAs 及其靶基因有通过对不同氮形态处理的响应而改变其根尖形态结构的能力。

对差异靶基因KEGG 通路分析表明，一些靶基因富集参与氮代谢途径，如NADH-dependent glutamate synthase family protein和nitrate transmembrane transporter(NRT1.1) 。NADHdependent glutamate synthase family protein酶作为铵态氮吸收同化过程中的关键酶对植物幼苗期根系的初级铵同化很重要[29]。另有研究表明，在水稻NADH-gogat1突变体中，相比较于1 mmol·L−1硝态氮处理，1 mmol·L−1铵态氮处理抑制其根系的生长[29]，该研究结果与本研究结果相一致。相比较于硝态氮处理，铵态氮处理抑制靶基因NADH-dependent glutamate synthase family protein表达，从而导致杨树根尖发育受到抑制，这表明NADH-dependentglutamate synthase family protein 可能在植物根尖响应不同氮形态时发挥重要的作用。NRT1.1作为硝态氮转运体对硝态氮的反应范围广泛，从初始硝酸盐反应到长期发育变化，且对侧根发育的影响发挥着重大作用，即NRT1.1在低浓度硝态氮条件下抑制侧根生长，而在高浓度硝态氮条件下促进植物侧根生长[30]，该结论与本文研究结果一致。在本研究中，相比较于硝态氮处理，铵态氮处理下靶基因NRT1.1上调表达，抑制了杨树侧根的伸长生长。另外，对灰杨根尖氮代谢以及生长发育相关差异表达靶基因进行MapMan 分析表明，NAC基因家族出现频率最多。在拟南芥中，NAC4基因作为protein auxin signaling F-BOX 3(AFB3)的下游基因，参与硝态氮响应，从而影响根系形态结构[31]。在小麦中，TaNAC2-5A可直接结合硝酸盐转运体和谷氨酰胺合成酶基因的启动子区域，过表达TaNAC2-5A可促进小麦根系生长和硝酸盐流入[32]。在本研究中，相比较于铵态氮处理，硝态氮处理中NAC1基因高表达，从而促进根生长，这一结论与前人的研究结果相一致。

在联合分析中，共发现24 对miRNA-靶基因呈负相关关系，其中，miR396-GRF是联合分析中出现最多的一组关系对。GRF4作为miR396 的靶基因能够驱动水稻根对铵态氮的吸收[23]，同时GRF4能够驱动硝态氮转运体的转录水平，例如NRT1.1B和NRT2.3a、GRF4还能够驱动硝态氮同化酶基因nitrate reductase 1(NIA1)、NIA3和nitrite reductase 1(NiR1)的合成从而去调节氮代谢，且在水稻中，GRF4突变抑制植物生长[23]。在本研究中，相比较于硝态氮处理，在铵态氮处理中miR396下调表达，其靶基因GRF1/2/5/9均上调表达，从而促进杨树根生长，这一结论与前人的研究结果相一致。该结果说明miR396-GRF在不同氮形态胁迫响应中发挥着重要的调控作用。另外，表达差异最显著的靶基因为新miRNAPC-3p-42422_177 的3个靶基因，但其功能未知，有待进一步研究证实。

4 结论

本研究从2种不同氮形态处理杨树根尖miRNAs文库中鉴定出523个已知miRNAs 和42个新miRNAs，其中，有96个miRNAs 显著差异表达。与硝态氮处理相比，铵态氮处理下有44个miRNAs 上调表达，52个下调表达，其中，miR396-GRF模块引起笔者的关注。研究表明，miR396-GRF模块可能通过响应不同氮形态调控杨树根系形态构型。该研究不仅可以为我国转基因杨树研究储备基因资源信息，还可以加快氮吸收能力强的杨树优质速生良种的培育。