白樺BpNAC012基因調(diào)控白樺木質(zhì)部發(fā)育表達譜分析

郭依萍 劉佳欣 于 穎 王 超 楊傳平

（林木遺傳育種國家重點實驗室，東北林業(yè)大學(xué)，哈爾濱 150040）

NAC 基因是植物特異性轉(zhuǎn)錄因子基因家族之一，其命名是由矮牽牛基因NAM［1］以及擬南芥基因ATAF 和CUC2 的首字母縮略詞構(gòu)成［2］。NAC 轉(zhuǎn)錄因子具有多種功能，如參與植物次生長，在細(xì)胞分裂和植株衰老中發(fā)揮作用，參與激素調(diào)控和信號轉(zhuǎn)導(dǎo)［3～4］。植物的次生生長包括維管組織形成、次生細(xì)胞壁形成、木質(zhì)化、細(xì)胞程序化死亡以及心材形成等過程［5～6］。近年研究發(fā)現(xiàn)多個NAC 基因?qū)?xì)胞次生壁的形成起著正調(diào)控作用。例如擬南芥SND1（Secondary wall-associated NAC domain 1），NST1/2/3（NAC secondary wall thicken promoting factor 1/2/3）和VND6/7（Vascular-related NAC domain 6/7）。在擬南芥次生壁正常形成的過程中，SND2/3、MYB103/85/52/54/69/42/43/20 和KNAT7等11 個SND1 轉(zhuǎn)錄因子是必需的。抑制SND2/3、MYB103/85/52/54 和KNAT7 的表達，能顯著減少纖維細(xì)胞次生壁增厚；而SND2/3 和MYB103 的過量表達，則能促進纖維中次生壁的增厚。SND1 在莖的維管束間纖維和木質(zhì)部纖維中特異表達，過表達SND1促進非厚壁細(xì)胞中次生壁的沉積；而抑制SND1 的表達，纖維中缺失次生壁［7～8］。NAC 轉(zhuǎn)錄因子對次生壁合成的調(diào)控是通過對下游靶基因的調(diào)控完成的。例如，VND7 能激活下游轉(zhuǎn)錄因子基因及一些參與次生壁形成、細(xì)胞壁化學(xué)修飾和細(xì)胞凋亡等非轉(zhuǎn)錄因子基因的表達［9］。

前期研究結(jié)果顯示，BpNAC012 基因調(diào)控白樺（Betula platyphylla Suk.）次生細(xì)胞壁的合成［10］。本研究主要以次生生長相關(guān)的白樺BpNAC012 基因為研究對象，對其在白樺木質(zhì)部中調(diào)控的下游基因進行研究。通過對過表達、抑制表達轉(zhuǎn)基因株系及對照株系莖部轉(zhuǎn)錄組進行分析，鑒定差異表達基因，初步確定BpNAC012基因調(diào)控白樺木質(zhì)部發(fā)育的基因表達途徑，為白樺的遺傳改良奠定基礎(chǔ)。

1 材料與方法

1.1 實驗材料

分別取3 個月苗齡的轉(zhuǎn)基因白樺和野生型白樺，去除葉片和根，作為轉(zhuǎn)錄組測序材料，每個處理3 次生物學(xué)重復(fù)，每次生物學(xué)重復(fù)5～7 株（抑制表達株系莖小，7 株可滿足建庫要求）。WT為野生型株系，OE為基因過表達株系，RS為抑制表達株系。

1.2 建庫測序流程

在RNA 樣品檢測合格后，利用帶有Oligo（dT）的磁珠富集總RNA 中的mRNA。隨后將mRNA 打斷成短片段，以mRNA 為模板，利用random hexamers 合成cDNA 第一鏈，利用PCR 合成cDNA 第二鏈，AMPure XP beads 純化雙鏈cDNA。PCR 擴增獲得測序所用文庫。為保證文庫質(zhì)量符合測序要求，使用Qubit2.0 對文庫進行初步定量，隨后使用Agilent 2100 對文庫的插入長度進行檢測，插入長度符合要求后，使用Q-PCR 方法準(zhǔn)確定量文庫的有效濃度后（有效濃度＞2 nmol·L-1），進行Illumina HiSeq測序。

1.3 生物信息學(xué)分析

利用Corset程序?qū)D(zhuǎn)錄本進行層次聚類，以聚類后序列為參考，進行后續(xù)分析。

1.3.1 測序質(zhì)量評估

通過測序錯誤率分布檢查來檢測測序數(shù)據(jù)的質(zhì)量，用Qphred=-10log10（e）來衡量，其中測序錯誤率用e表示。將測序得到的原始測序序列（Sequenced Reads）或者raw reads去除帶接頭（adapter）的reads；去除N（N表示無法確定堿基信息）比例大于10%的reads；去除低質(zhì)量reads（質(zhì)量值Qphred ＜=20的堿基數(shù)占整個reads的50%以上的reads）。

1.3.2 轉(zhuǎn)錄本拼接

對于無參考基因組的項目，獲得clean reads后，采用Trinity 進行轉(zhuǎn)錄本拼接。利用比對上轉(zhuǎn)錄本的reads數(shù)和表達模式對轉(zhuǎn)錄本進行Corset層次聚類（https：//code.google.com/p/corset-project/）。對轉(zhuǎn)錄本及聚類序列長度分別進行拼接轉(zhuǎn)錄本長度分布統(tǒng)計。

1.3.3 基因功能注釋及CDS預(yù)測

利用生物信息學(xué)數(shù)據(jù)庫進行基因功能注釋，包括：Nr，Nt，Pfam，KOG/COG，Swiss-prot，KEGG，GO。從比對結(jié)果中提取出轉(zhuǎn)錄本的ORF 編碼框信息，并按照標(biāo)準(zhǔn)密碼子表將編碼區(qū)序列翻譯成氨基酸序列（按照5′-＞3'的順序）；對于沒有比對結(jié)果或者比對未預(yù)測出結(jié)果的序列，則采用estscan（3.0.3）軟件預(yù)測其ORF。

1.3.4 基因表達水平分析

采用RSEM 軟件將每個樣品的clean reads 往參考序列（ref）上做mapping。進行參考序列比對。比對結(jié)果進行統(tǒng)計，得到樣品比對的readcount 數(shù)目，并對其進行FPKM（expected number of Fragments Per Kilobase of transcript sequence per Millions base pairs sequenced）轉(zhuǎn)換，來進行基因的表達水平分析。

1.3.5 差異表達分析

利用基因表達水平分析中得到的readcount數(shù)據(jù)進行基因差異表達分析。采用DESeq2 軟件，篩選閾值為padj＜0.05 且|log2FoldChange|＞1，對原有假設(shè)檢驗得到的p-value進行校正。

1.3.6 差異基因富集分析

利用GOseq 進行GO 富集分析。以KEGG（Kyoto Encyclopedia of Genes and Genomes）Pathway 為單位，應(yīng)用超幾何檢驗，找出差異基因相對于所有有注釋的基因顯著富集的pathway，進行Pathway 顯著性富集分析。繪制KEGG 富集散點圖，KEGG 富集程度通過Rich factor、qvalue 和富集到此通路上的基因個數(shù)來衡量。其中Rich factor指差異表達的基因中位于該pathway 條目的基因數(shù)目與所有有注釋基因中位于該pathway 條目的基因總數(shù)的比值。qvalue 是做過多重假設(shè)檢驗校正之后的Pvalue，qvalue 的取值范圍為［0，1］，越接近于零，表示富集越顯著。

2 結(jié)果與分析

2.1 RNA質(zhì)量分析

RNA 質(zhì)量決定著文庫的質(zhì)量，RNA 提取質(zhì)量分析結(jié)果顯示如圖1，樣品質(zhì)量滿足建庫測序要求，且總量滿足2次以上建庫需要。

2.2 表達譜測序結(jié)果評估

2.2.1 測序質(zhì)量評估

測序過程存在錯誤率，測序錯誤率分布檢查可以反映測序數(shù)據(jù)的質(zhì)量。一般情況下，單個堿基位置的測序錯誤率應(yīng)該低于1%。從表1可以看出，本文庫測序錯誤率均低于0.1%，說明測序質(zhì)量符合要求。

從表2 可以看出，6 個文庫的G/C 含量符合正常規(guī)律。接頭序列、未知序列、低質(zhì)量序列等不合格序列的百分比低于2%，cleanread 應(yīng)占原始數(shù)據(jù)的90%以上，這樣的數(shù)據(jù)為合格數(shù)據(jù)。從表1 和附圖1 可以看出，該文庫cleanread 占原始數(shù)據(jù)的比例均符合后續(xù)分析的要求。

2.2.2 轉(zhuǎn)錄本拼接

對轉(zhuǎn)錄本及聚類序列長度分別進行統(tǒng)計，結(jié)果見表2、附表1 和附圖2。從圖表中可以看出，本次分析轉(zhuǎn)錄組的基因長度分布符合進一步分析的要求。

表1 數(shù)據(jù)產(chǎn)出質(zhì)量情況一覽表Table 1 List of data output quality

表2 拼接長度頻數(shù)分布情況Table 2 List of frequency distribution of splice length

2.3 基因功能注釋

為獲得基因功能全面信息，對本研究文庫獲得的數(shù)據(jù)進行7 大數(shù)據(jù)庫基因功能注釋，包括：Nr，Nt，Pfam，KOG/COG，Swiss-prot，KEGG，GO 數(shù)據(jù)庫。本項目數(shù)據(jù)在7 大數(shù)據(jù)庫中的注釋成功率情況見表3。從注釋成功率可以看出，獲得功能注釋的基因比率符合后續(xù)基因功能分析的要求。

表3 基因注釋成功率統(tǒng)計Table 3 Gene annotation success rate statistics

2.4 基因表達水平分析

2.4.1 參考序列比對

參考序列比對統(tǒng)計結(jié)果見附表2。在RNAseq 技術(shù)中，F(xiàn)PKM 分析是目前最為常用的基因表達水平估算方法。所以我們將readcount數(shù)進行了FPKM轉(zhuǎn)換用于下一步差異表達分析。

2.4.2 基因差異表達分析及差異基因篩選

本研究3個處理間共獲得差異表達基因1 476條。從圖2 中可以看出，與對照相比，過表達株系OE 中上調(diào)的基因有627 條，下調(diào)的基因有229 條，抑制表達株系中上調(diào)的基因有299條，下調(diào)表達的基因有207 條。說明在轉(zhuǎn)基因株系中，BpNAC012基因的表達影響一系列的基因表達變化，進而調(diào)控木質(zhì)部發(fā)育性狀的變化。

2.5 差異基因GO富集分析

GO 功能顯著性富集分析給出與基因組背景相比，在差異表達基因中顯著富集的GO 功能條目，從而給出差異表達基因與哪些生物學(xué)功能顯著相關(guān)。結(jié)果如表4 所示。差異表達基因主要富集在轉(zhuǎn)錄因子相關(guān)的GO 分類中。說明轉(zhuǎn)錄因子在本實驗樣本的基因表達差異及性狀差異中起到重要作用。統(tǒng)計被顯著富集的各個GOterm 中的基因數(shù)，以柱狀圖的形式展示，如圖3 所示。結(jié)合差異表達基因數(shù)量分析，從圖中可以看出，過表達BpNAC012 能夠調(diào)控更多的基因表達變化。而抑制表達BpNAC012 更多的是影響蛋白修飾和轉(zhuǎn)運類基因的表達。

表4 樣品中差異基因富集的Gene Ontology分類Table 4 Gene Ontology classification of differential gene enrichment in samples

2.6 差異基因KEGG富集分析

差異基因KEGG 富集散點圖（見圖4）圖形化展示了KEGG 富集分析結(jié)果。從圖4中可以看出，涉及受體信號通路、營養(yǎng)代謝、氨基酸合成及苯丙烷生物合成相關(guān)代謝通路基因比較富集。其中苯丙烷生物合成，是次生細(xì)胞壁合成的關(guān)鍵代謝通路，因此說明BpNAC012 的表達，調(diào)控了次生細(xì)胞壁合成的代謝通路。

2.7 表達差異最大的前十位基因分析

對過表達株系相對于野生型表達量上調(diào)最大的10個基因和抑制表達株系相對于野生型表達量下調(diào)最大的10 個基因進行分析（見表5～6）?？勺C明這些基因強烈受BpNAC012 的調(diào)控，并在過表達或抑制表達株系發(fā)育的生理程序中起著重要作用。

2.8 BpNAC012 調(diào)控木質(zhì)部發(fā)育及次生壁合成相關(guān)基因的表達

在BpNAC012過表達株系中，鑒定到一些和纖維素合成相關(guān)的基因（見表7），和一些細(xì)胞壁相關(guān)酶類。同時，鑒定了一些與木質(zhì)素合成及沉積相關(guān)的基因，和木質(zhì)部發(fā)育相關(guān)的基因。這些基因可能受BpNAC012 的誘導(dǎo)，參與過表達轉(zhuǎn)基因白樺木質(zhì)部發(fā)育及次生細(xì)胞壁的合成。

在抑制表達株系中，鑒定到一些木質(zhì)部發(fā)育及細(xì)胞壁合成相關(guān)基因下調(diào)表達（見表8），如一些與纖維素合成相關(guān)的基因下調(diào)表達。

2.9 BpNAC012 調(diào)控其它轉(zhuǎn)錄因子基因的差異表達

在過表達和抑制表達株系中，分別鑒定到了與野生型對照差異表達的轉(zhuǎn)錄因子（見表9～11）。其中，與野生型對照相比，在過表達株系中鑒定了97 個上調(diào)表達的轉(zhuǎn)錄因子，15 個下調(diào)表達的轉(zhuǎn)錄因子，在抑制表達株系中，鑒定了42個上調(diào)表達的轉(zhuǎn)錄因子，15個下調(diào)表達的轉(zhuǎn)錄因子。

表6 抑制表達株系中差異表達倍數(shù)最高的前10位基因Table 6 The top 10 genes with the highest differential expression folds in the suppressed expression lines

3 討論

BpNAC012 對白樺次生細(xì)胞壁合成的調(diào)控是通過下游基因表達變化實現(xiàn)的，與對照相比，過表達株系OE 中上調(diào)的基因有627 條，下調(diào)的基因有229 條，抑制表達株系中上調(diào)的基因有299 條，下調(diào)表達的基因有207 條。兩個轉(zhuǎn)錄組中共有140條基因差異表達。說明在轉(zhuǎn)基因株系中，Bp?NAC012 基因的表達影響了一系列的基因表達變化，進而調(diào)控次生壁的合成，這些差異表達基因可能參與白樺次生細(xì)胞壁及木質(zhì)部形成的調(diào)控。

在轉(zhuǎn)基因株系中，表達變化最明顯的基因，可能對于BpNAC012 轉(zhuǎn)基因造成的變化起著重要作用，對過表達株系相對于野生型表達量上調(diào)最大的10個基因和抑制表達株系相對于野生型表達量下調(diào)最大的10個基因進行分析。UTP—glucose-1-phosphate uridylyltransferase，EC number 2.7.7.9，同義名稱UDP glucose pyrophosphorylase，UDPG pyrophosphorylase，UGPase［11］，是 一 種參 與 糖 代 謝 的酶。它以葡萄糖-1-磷酸和UTP 為底物合成UDP

葡萄糖；UTP-葡萄糖-1-磷酸尿苷酰基轉(zhuǎn)移酶是在所有三個領(lǐng)域（細(xì)菌、真核生物和古細(xì)菌）中發(fā)現(xiàn)的酶，因為它是糖原生成和細(xì)胞壁合成的關(guān)鍵參與者［12］。UDP-葡萄糖焦磷酸化酶基因過表達可提高黃麻（Corchorus capsularis L.）的纖維素含量［13］。說明UGPase 在營養(yǎng)生長、細(xì)胞壁形成及纖維素合成中起重要作用。在BpNAC012抑制表達株系中，該基因顯著下調(diào)表達，可能與抑制表達株系纖維素含量下降有關(guān)。值得注意的是No lysine kinase 1 isoform 1，同時出現(xiàn)在過表達株系上調(diào)倍數(shù)前十中，和抑制表達株系下調(diào)倍數(shù)前十中，說明該基因與BpNAC012 的表達密切協(xié)同，受BpNAC012 的強調(diào)控。No lysine kinase 1（WNK）具有絲氨酸/蘇氨酸激酶的催化賴氨酸的獨特特征，因而得名。WNK1 在植物體內(nèi)離子轉(zhuǎn)運的調(diào)節(jié)、激素對滲透和離子脅迫的調(diào)節(jié)、晝夜節(jié)律和開花時間、蛋白質(zhì)轉(zhuǎn)運中、離子通道調(diào)控等過程中起重要作用［14］。因此說明BpNAC012 表達的改變，在某種程度上可能影響了轉(zhuǎn)基因株系的信號傳導(dǎo)途徑。

表7 過表達株系木質(zhì)部發(fā)育相關(guān)上調(diào)表達基因Table 7 Genes related to up-regulated expression of xylem development in overexpressing lines

表8 抑制表達株系木質(zhì)部發(fā)育相關(guān)下調(diào)表達基因Table 8 Genes related to down-regulated expression of xylem development-inhibited expression lines

表9 過表達株系中差異表達轉(zhuǎn)錄因子Table 9 Differentially expressed transcription factors in overexpressing lines

表10 抑制表達株系中差異表達轉(zhuǎn)錄因子Table 10 Inhibition of differentially expressed transcription factors in expressing lines

過量表達BpNAC012可促進木質(zhì)部發(fā)育，次生壁增厚。因此本研究進行差異表達基因分析時，關(guān)注了和這些性狀相關(guān)的基因。CesA是細(xì)胞壁生物合成過程中催化纖維素合成的關(guān)鍵酶。該蛋白參與纖維素合成代謝途徑，參與次生細(xì)胞壁的形成，是木質(zhì)部細(xì)胞壁增厚所需［15］。本研究中，Bp-NAC012 過表達株系中鑒定的2 條CesA 相關(guān)基因的表達量比對照高，而在抑制表達株系中鑒定的1條CesA 基因的表達量比野生型對照低。表明Bp-NAC012 通過增強過表達株系中CesA 基因的表達促進纖維素生物合成。細(xì)胞壁結(jié)構(gòu)的調(diào)節(jié)在植物的生長和分化中起著關(guān)鍵的作用。Xyloglucan endotransglucosylase/hydrolase family protein（XTHs）催化木葡聚糖鏈的分解和分子連接，在細(xì)胞壁的松弛和重新排列中起作用。本研究在過表達株系中鑒定了2 條上調(diào)表達的Xyloglucan endotransglucosylase/hydrolase family protein，綜合分析結(jié)果說明，該基因受BpNAC012 的誘導(dǎo)，在白樺木質(zhì)部發(fā)育中起重要的作用。endo-1，4-beta-glucanase 家族基因在植物生長和發(fā)育的多方面具有功能，能促進植物生長和側(cè)根發(fā)育，增強細(xì)胞擴展。楊樹的endo-1，4-beta-glucanase 是木質(zhì)部組織特異的，在纖維素生物合成中起作用，能夠促進植物的生長［16］。本研究中，在抑制表達株系中，鑒定了1 個下調(diào)表達的endo-1，4-beta-glucanase 基因，綜合分析，說明endo-1，4-beta-glucanase 隨著BpNAC012表達的下調(diào)而下調(diào)，受其調(diào)控，和抑制表達轉(zhuǎn)基因株系的木質(zhì)部發(fā)育及纖維素合成受阻相關(guān)。

本研究在過表達和抑制表達株系中，分別鑒定到與野生型對照差異表達的轉(zhuǎn)錄因子AP2-EREBP、Tify、WRKY、Orphans、bHLH、NAC、GRAS、C2H2、SWI/SNF、SET、MYB、HSF、AUX/IAA、ARF、Sigma70-like、mTERF、G2-like、C3H、HB、FAR1、CPP、ABI3VP1、TRAF、GNAT、C2C2-Dof和LOB 等。在過表達株系中，差異表達數(shù)量最多的轉(zhuǎn)錄因子是AP2-EREBP 類轉(zhuǎn)錄因子，而且全部是上調(diào)表達；其次是Tify 類轉(zhuǎn)錄因子，WRKY 類轉(zhuǎn)錄因子；此外bHLH、NAC、AUX/IAA 和MYB 轉(zhuǎn)錄子表現(xiàn)了不同的表達模式。在抑制表達株系中，上調(diào)表達的轉(zhuǎn)錄因子數(shù)量明顯少于過表達株系，其中差異表達數(shù)量最多的，依然是AP2-EREBP 類轉(zhuǎn)錄因子，并且有1 個下調(diào)；上調(diào)表達的WRKY、Tify和Orphars類轉(zhuǎn)錄因子數(shù)量也明顯下降；此外，在過表達株系中都上調(diào)的轉(zhuǎn)錄因子，在抑制表達株系中出現(xiàn)了下調(diào)表達，如C3H 和C2H2。在過表達株系中，GRAS、SWI.SNF 和HSF 類轉(zhuǎn)錄因子都有大于3條的上調(diào)表達，而在抑制表達株系中沒有檢測到；在抑制表達株系中，TRAF、C2C2-dof 和GNAT類轉(zhuǎn)錄因子都下調(diào)表達，而在過表達表達株系中沒有檢測到。以上結(jié)果說明，在BpNAC012調(diào)控的轉(zhuǎn)基因株系中，存在多個層級的轉(zhuǎn)錄調(diào)控，BpNAC012能夠調(diào)控其他轉(zhuǎn)錄因子的表達；在過表達和抑制表達株系中，存在不同的、復(fù)雜的調(diào)控模式，一些特異轉(zhuǎn)錄因子，在各自的性狀變化中起著相應(yīng)的調(diào)控作用。