白樺BpNAC012基因調控白樺木質部發育表達譜分析

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

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

NAC 基因是植物特異性轉錄因子基因家族之一，其命名是由矮牽?；騈AM［1］以及擬南芥基因ATAF 和CUC2 的首字母縮略詞構成［2］。NAC 轉錄因子具有多種功能，如參與植物次生長，在細胞分裂和植株衰老中發揮作用，參與激素調控和信號轉導［3～4］。植物的次生生長包括維管組織形成、次生細胞壁形成、木質化、細胞程序化死亡以及心材形成等過程［5～6］。近年研究發現多個NAC 基因對細胞次生壁的形成起著正調控作用。例如擬南芥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 轉錄因子是必需的。抑制SND2/3、MYB103/85/52/54 和KNAT7 的表達，能顯著減少纖維細胞次生壁增厚；而SND2/3 和MYB103 的過量表達，則能促進纖維中次生壁的增厚。SND1 在莖的維管束間纖維和木質部纖維中特異表達，過表達SND1促進非厚壁細胞中次生壁的沉積；而抑制SND1 的表達，纖維中缺失次生壁［7～8］。NAC 轉錄因子對次生壁合成的調控是通過對下游靶基因的調控完成的。例如，VND7 能激活下游轉錄因子基因及一些參與次生壁形成、細胞壁化學修飾和細胞凋亡等非轉錄因子基因的表達［9］。

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

1 材料與方法

1.1 實驗材料

分別取3 個月苗齡的轉基因白樺和野生型白樺，去除葉片和根，作為轉錄組測序材料，每個處理3 次生物學重復，每次生物學重復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 擴增獲得測序所用文庫。為保證文庫質量符合測序要求，使用Qubit2.0 對文庫進行初步定量，隨后使用Agilent 2100 對文庫的插入長度進行檢測，插入長度符合要求后，使用Q-PCR 方法準確定量文庫的有效濃度后（有效濃度＞2 nmol·L-1），進行Illumina HiSeq測序。

1.3 生物信息學分析

利用Corset程序對轉錄本進行層次聚類，以聚類后序列為參考，進行后續分析。

1.3.1 測序質量評估

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

1.3.2 轉錄本拼接

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

1.3.3 基因功能注釋及CDS預測

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

1.3.4 基因表達水平分析

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

1.3.5 差異表達分析

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

1.3.6 差異基因富集分析

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

2 結果與分析

2.1 RNA質量分析

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

2.2 表達譜測序結果評估

2.2.1 測序質量評估

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

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

2.2.2 轉錄本拼接

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

表1 數據產出質量情況一覽表Table 1 List of data output quality

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

2.3 基因功能注釋

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

表3 基因注釋成功率統計Table 3 Gene annotation success rate statistics

2.4 基因表達水平分析

2.4.1 參考序列比對

參考序列比對統計結果見附表2。在RNAseq 技術中，FPKM 分析是目前最為常用的基因表達水平估算方法。所以我們將readcount數進行了FPKM轉換用于下一步差異表達分析。

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

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

2.5 差異基因GO富集分析

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

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

2.6 差異基因KEGG富集分析

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

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

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

2.8 BpNAC012 調控木質部發育及次生壁合成相關基因的表達

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

在抑制表達株系中，鑒定到一些木質部發育及細胞壁合成相關基因下調表達（見表8），如一些與纖維素合成相關的基因下調表達。

2.9 BpNAC012 調控其它轉錄因子基因的差異表達

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

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

3 討論

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

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

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

表7 過表達株系木質部發育相關上調表達基因Table 7 Genes related to up-regulated expression of xylem development in overexpressing lines

表8 抑制表達株系木質部發育相關下調表達基因Table 8 Genes related to down-regulated expression of xylem development-inhibited expression lines

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

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

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

本研究在過表達和抑制表達株系中，分別鑒定到與野生型對照差異表達的轉錄因子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 等。在過表達株系中，差異表達數量最多的轉錄因子是AP2-EREBP 類轉錄因子，而且全部是上調表達；其次是Tify 類轉錄因子，WRKY 類轉錄因子；此外bHLH、NAC、AUX/IAA 和MYB 轉錄子表現了不同的表達模式。在抑制表達株系中，上調表達的轉錄因子數量明顯少于過表達株系，其中差異表達數量最多的，依然是AP2-EREBP 類轉錄因子，并且有1 個下調；上調表達的WRKY、Tify和Orphars類轉錄因子數量也明顯下降；此外，在過表達株系中都上調的轉錄因子，在抑制表達株系中出現了下調表達，如C3H 和C2H2。在過表達株系中，GRAS、SWI.SNF 和HSF 類轉錄因子都有大于3條的上調表達，而在抑制表達株系中沒有檢測到；在抑制表達株系中，TRAF、C2C2-dof 和GNAT類轉錄因子都下調表達，而在過表達表達株系中沒有檢測到。以上結果說明，在BpNAC012調控的轉基因株系中，存在多個層級的轉錄調控，BpNAC012能夠調控其他轉錄因子的表達；在過表達和抑制表達株系中，存在不同的、復雜的調控模式，一些特異轉錄因子，在各自的性狀變化中起著相應的調控作用。