基于硝態氮或銨態氮條件下楊樹根尖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模塊可能通過響應不同氮形態調控楊樹根系形態構型。該研究不僅可以為我國轉基因楊樹研究儲備基因資源信息，還可以加快氮吸收能力強的楊樹優質速生良種的培育。