陳 沫，何沙娥，陳少雄，歐陽林男，張 程，張維耀

(1.國家林業(yè)和草原局桉樹研究開發(fā)中心，廣東 湛江 524022；2.南京林業(yè)大學，江蘇 南京 210037)

合理的林分密度能夠顯著促進林木的徑向生長[1-2]，從而提高木材產(chǎn)量，解析種植密度促進林木徑向生長的內在機制意義重大。林木的徑向生長（即次生木質部生成）來源于維管形成層的活動。維管形成層細胞向外分化產(chǎn)生次生韌皮部，向內分化產(chǎn)生次生木質部，其中，次生木質部細胞不斷積累，使植株徑向增粗[3]。由維管形成層分化產(chǎn)生成熟的次生木質部細胞需要經(jīng)歷木質部母細胞分化和分裂、木質部細胞增大、次生壁加厚、細胞程序化凋亡和心材形成多個階段[4]，各個階段在分子水平上均受到一大批基因的調控。因此，研究木質部在不同種植密度下的基因表達模式，挖掘和鑒定響應種植密度的木質部生成關鍵基因并解析其作用機制是研究種植密度影響林木徑向生長內在機制的重要內容。

以往以模式植物擬南芥（Arabidopsis thaliana(Linn.) Heynh.）和楊屬（Populus）為研究對象，分離鑒定了一些參與木質部生成各個階段的關鍵調控基因，在木質部生成方面取得了重要進展。諸多研究表明，調控因子如受體激酶和轉錄因子等在木質部生成過程中發(fā)揮重要的調控作用，如PTL（PETAL LOSS，花瓣缺失）轉錄因子作為抑制子參與調控擬南芥根系發(fā)育過程中的形成層細胞分裂和木質部細胞分化活動[5]；類受體激酶PXY 和其家族成員PXL1、PXL2 參與調控擬南芥維管組織的細胞分裂活動[6-7]；類受體激酶FEI 和SHOU 參與次生壁物質合成，fei突變體纖維素合成減少，而SHOU 能夠恢復fei突變體表型[8-9]。次生壁合成也需要一些轉錄因子參與，Chen 等[10]在毛果楊（Populus trichocarpaTorr.&Gray）中繪制了一個4 層級的細胞壁合成轉錄調控網(wǎng)絡，其中包含17 個轉錄因子和27 個細胞壁成分合成基因，它們協(xié)同調控木質部生成；最近在楊屬、松屬（Pinus）等樹種中又相繼證實了多個NAC 和MYB 家族成員參與木材細胞壁物質合成過程[11-13]。由此可見，次生木質部生成的分子機制十分復雜，仍然有許多未知亟待解析。

桉屬（Eucalyptus）與楊屬為世界公認的速生人工林樹種，但二者木材的生物量和質量存在明顯差異。這2 個樹種在木質部生成調控機制上可能存在差異，但目前關于桉屬植物木質部生成相關重要基因的解析及其調控機制的研究相對較少[14-15]，同時，關于種植密度影響徑向生長的分子基礎研究尚未見報道。因此，本研究以不同種植密度下尾巨桉（Eucalyptus urophylla×EucalyptusgrandisW.Hill ex Maiden）的主莖為試驗材料，聯(lián)合PacBio Iso-Seq 和RNA-Seq 轉錄組技術，通過對尾巨桉木質部差異表達轉錄本（DETs）的表達模式展開分析，挖掘和鑒定桉樹中響應種植密度的木質部生成關鍵基因，以期為科學闡明種植密度影響林木徑向生長的分子機制奠定理論基礎。

1 材料與方法

1.1 試驗材料

試驗地位于廣東省湛江市南方國家級林木種苗示范基地內（111°38′ E，21°30′ N）。試驗材料為長勢均一、主莖通直、無病蟲害的尾巨桉DH32-29 無性系組織培育苗。2019 年7 月14 日，將幼苗移植至高0.5 m 直徑0.3 m 的種植袋中，并按照10 株·m?2高種植密度（株行距0.3 m×0.3 m）和5 株·m?2低種植密度（株行距0.6 m×0.6 m）移放至長×寬×高為5.0 m×4.5 m×0.5 m 的樣地內，每個處理設置3 個重復，共計6 個樣地。在幼苗生長過程中，水分和營養(yǎng)充分供應且保持一致。苗木種植后定期測定生長性狀，于直徑及其增長率表現(xiàn)出顯著差異的時期進行取樣（2019 年10 月15 日）（表1）。每個樣地選擇1 株平均木，取距離地面10 cm 處主莖的木質部、韌皮部，共計12 個樣本。液氮冷凍后，帶回試驗室置于?80℃的超低溫冰箱保存?zhèn)溆?。樣品命名為H-X1～3、L-X1～3、HP1～3、L-P1～3，其中，H 為高種植密度，L 為低種植密度，X 表示木質部，P 表示韌皮部，數(shù)字1～3 代表生物學重復。

表1 種植密度對尾巨桉幼苗生長性狀的影響Table 1 Effect of planting density on growth traits for seedlings of E.urophylla×E.grandis

1.2 轉錄組測序及結果分析

使用植物RNA 提取試劑盒DP130508（天根生物公司）分別提取12 個樣本的RNA，等量混合后使用SMARTer? PCR cDNA Synthesis Kit 試劑盒合成mRNA 的全長cDNA，通過PacBio 平臺進行全長轉錄組測序（百邁客公司）；使用CD-HIT軟件[16]去除轉錄本中的冗余序列；使用BLAST 軟件(version 2.2.26)[17]將得到非冗余轉錄本序列和NR、eggNOG、Pfam、GO、Swissprot、KOG、COG、KEGG 數(shù)據(jù)庫比對進行轉錄本功能注釋；使用iTAK 軟件[18]進行調控因子預測。

對12 個樣本的RNA 分別進行cDNA 合成和文庫構建，委托百邁客公司進行Illumina 二代測序。使用RSEM[19]軟件，通過比對到三代轉錄本上的位置信息，對轉錄本的表達水平進行定量，采用FPKM（Fragments Per Kb of transcript per Million fragments mapped）作為衡量指標。

1.3 DETs 篩選及分析

使用DESeq2[20]對樣本H-X1、H-X2、H-X3和L-X1、L-X2、L-X3 進行分析，獲得2 個種植密度間DETs。篩選標準為差異倍數(shù)（Fold Change）≥2 且錯誤發(fā)現(xiàn)率（FDR）<0.01；使用百邁客云平臺熱圖繪制工具繪制DETs 熱圖，每個表達量的數(shù)值均取以10 為底數(shù)的對數(shù)（Log10）[21]；使用STRING 數(shù)據(jù)庫[22]對DETs 構建蛋白互作網(wǎng)絡。

1.4 DETs 實時熒光定量PCR 驗證

以轉錄組試驗中2 個種植密度下尾巨桉幼苗的木質部RNA 為模板，使用TSK302S RT6 cDNA Synthesis Kit Ver 2（擎科生物公司）合成cDNA 模板，從443 個DETs 中隨機選取9 個DETs，進行qRT-PCR 驗證轉錄組測序結果的準確性；選擇9 個DETs 進行組織表達分析，材料為尾巨桉幼苗的葉片、韌皮部、第2 節(jié)間木質部（木質部2）和第7 節(jié)間木質部（木質部7）4 個部位，RNA 的提取采用植物RNA 提取試劑盒V1.5(百菲特生物公司)。利用Primer 5.0設計qRT-PCR引物，以arp4(Actin-related protein 4)編碼基因作為內參，引物見表2。使用2×T5 Fast qPCR Mix(SYBR Green I)熒光定量試劑盒（擎科生物公司）進行熒光定量檢測。每個樣品設置4 個技術重復，采用2?ΔΔCt計算相對表達量。使用SPSS 25.0 處理數(shù)據(jù)，使用Excel 2016 軟件繪圖，以P<0.01 表示差異顯著。

表2 qRT-PCR 引物序列Table 2 Primers for qRT-PCR

2 結果與分析

2.1 轉錄組測序統(tǒng)計

對高、低2 個種植密度下尾巨桉主莖木質部和韌皮部進行Pacbio 全長轉錄組測序，共獲得247 536條環(huán)形一致序列（CCS），其中，209 654 條為全長非嵌合序列；對全長非嵌合序列進行聚類分析，獲得81 317 條一致性序列；進一步對一致性序列進行校正，最終獲得高質量轉錄本（準確度＞99%）79 965 條，低質量轉錄本1 281 條，高質量轉錄本達到98%以上；經(jīng)CD-HIT 去冗余后，獲得非冗余轉錄本序列45 490條（表3）。

表3 PacBio Iso-Seq 測序數(shù)據(jù)評估統(tǒng)計表Table 3 Statistical table of PacBio Iso-Seq data evaluation

對高、低兩個種植密度下尾巨桉主莖木質部進行RNA-Seq 測序，共獲得130 292 914 條150 bp雙端Reads。6 個樣本中，測序獲得最多的Clean Reads 為23 322 440 條，最少為19 998 315 條；獲得Clean Data 堿基數(shù)最多為6 988 132 460 bp，最少為5 990 895 036 bp；GC 含量為50.30%～51.05%，%≥Q30（Clean Data 質量值≥30 的堿基所占的百分比）為93.88%～94.41%（表4），GC 含量和%≥Q30均值分別為50.60%、94.19%。

表4 RNA-Seq 測序數(shù)據(jù)評估統(tǒng)計表Table 4 Statistical table of RNA-Seq data evaluation

2.2 不同種植密度下尾巨桉主莖轉錄組分析

將三代轉錄組測序獲得的45 490 個非冗余轉錄本與NR、eggNOG、Pfam、GO、Swissprot、KOG、COG、KEGG 數(shù)據(jù)庫比對并進行功能注釋，有44 888 條轉錄本至少比對到一個數(shù)據(jù)庫中，注釋率達98.68%（表5），其中，NR 數(shù)據(jù)庫注釋率達98.47%，eggNOG、Pfam、Swissprot、GO、KOG數(shù)據(jù)庫注釋到的轉錄本數(shù)量均超過50%，KEGG數(shù)據(jù)庫注釋到的轉錄本數(shù)量最低。

表5 轉錄本在各數(shù)據(jù)注釋情況Table 5 Statistics of BLAST annotation of different databases

2.3 不同種植密度木質部轉錄本表達模式分析

為了充分了解不同種植密度下尾巨桉木質部轉錄本的表達模式，以H-X 為對照組，L-X 為試驗組，篩選獲得差異表達轉錄本443 個，其中，215個DETs 在低種植密度條件下上調表達，228 個DETs 下調表達。注釋結果顯示，443 個DETs 中包含了60 個假定的調控因子，歸為37 個家族（圖1A、圖2），其中，28 個差異表達調控因子在低種植密度下尾巨桉的木質部中上調表達，編碼21 類調控因子。這些調控因子與多種生物學過程有關，GARP-ARR-B、Trihelix 和RLK-Pelle_LRRXI-1 參與調控維管形成層活動；C3H、MYB、NAC和RLK-Pelle_LRR-XIIIa 參與調控細胞壁成分合成；RLK-Pelle_LRR-III 參與次生生長過程；RLK-Pelle_RLCK-V 參與細胞周期調控；GARP-G2-like 和WNK_NRBP 參與節(jié)律調控。32 個DETs 下調表達，所對應調控因子的功能多與外界脅迫有關，如熱激轉錄因子HSF、參與鹽脅迫的bHLH 和STE_STE7、參與干旱脅迫的RLK-Pelle_WAK_LRK10L-1、參與病蟲害應答的RLK-Pelle_LRR-VI-2 和RLKPelle_LRR-Xa、對脫落酸（ABA）敏感的bHLH 和RLK-Pelle_CrRLK1L-1 等。除編碼調控因子的轉錄本外，一些轉錄本的表達模式也發(fā)生了顯著的變化，在低種植密度植株的木質部中細胞壁物質合成相關蛋白（蔗糖合成酶，纖維素合成酶，漆酶，糖基轉移酶，糖基水解酶）和木質素物質合成相關編碼基因多為上調表達（圖1B）。

圖1 不同種植密度下尾巨桉木質部中80 條差異表達轉錄本的表達水平Fig.1 Heatmap showing the expression levels of 80 planting density-regulated differentially expressedtranscripts in xylem for E.urophylla×E.grandis

圖2 不同種植密度下尾巨桉木質部中60 條差異表達調控因子的類型分布Fig.2 Classification of 60 differentially expressed regulatory factors in xylem for E.urophylla×E.grandis under different planting densities

2.4 差異表達轉錄本蛋白互作網(wǎng)絡分析

通過STRING 數(shù)據(jù)庫對443 個DETs 進行蛋白互作網(wǎng)絡分析，檢測到一簇相互作用的基因簇，共包含30 個轉錄本（圖3，附表1）。這些基因主要與參與翻譯后修飾、蛋白質轉換和伴侶蛋白，細胞周期控制、細胞分裂和染色體分割，轉錄等功能。MIX-PB_transcript_67597 作為網(wǎng)絡的核心基因，與MIX-PB_transcript_30100 和MIX-PB_transcript_52911占據(jù)了網(wǎng)絡的中心位置（圖3）。數(shù)據(jù)庫注釋表明：MIX-PB_transcript_67597 編碼類受體激酶家族RLK-Pelle_LRR-XI-1 的PXL2，注釋功能為轉錄。轉錄本MIX-PB_transcript_52911 編碼蛋白CUL1 參與目標蛋白的泛素化；MIX-PB_transcript_30100 編碼蛋白T15D22.7 參與硫醇氧化酶和硫氧還蛋白家族蛋白的合成，二者在數(shù)據(jù)庫的注釋功能均為細胞周期控制、細胞分裂和染色體分割。這3 個基因在低種植密度植株的木質部中均上調表達，可能在尾巨桉主莖徑向生長中發(fā)揮關鍵作用（附表1）。

附表 1 蛋白互作網(wǎng)絡中30 個參與種植密度響應的差異表達轉錄本注釋Attached list 1:The annotation of 30 differentially expressed transcripts responding to planting density in protein interaction network

圖3 不同種植密度下尾巨桉木質部中差異表達轉錄本蛋白互作網(wǎng)絡Fig.3 The interaction network of proteins of differentially expressed transcripts in xylem for E.urophylla×E.grandis under different planting densities

2.5 qRT-PCR 驗證

以轉錄組測序中2 個種植密度下的木質部RNA為模板，隨機選擇9 個DETs 進行qRT-PCR 驗證，qRT-PCR 結果與轉錄組測序的表達結果基本一致（圖4），證明轉錄組測序數(shù)據(jù)準確可靠。

圖4 9 個差異表達轉錄本的qRT-PCR 驗證Fig.4 Confirmation of the expression profiles of 9 differentially expressed transcripts by qRT-PCR

為了驗證木質部生成候選基因的功能，通過qRT-PCR 對9 個木質部生成候選基因進行分析，檢測了它們在不同組織的基因表達情況。這9 個DETs 分別是轉錄因子基因PTL、MYB46、C3H14、NAC86，細胞壁合成酶基因CesA、LAC17以及PXL2、T15D22.7、CUL1。qRT-PCR 結果（圖5）：這些編碼基因在葉片中幾乎不表達，在韌皮部、木質部均有不同程度的表達，其中，PXL2在木質部2 中的表達量最高，在韌皮部和木質部7 中表達水平也較高；CUL1和T15D22.7在韌皮部和木質部中高水平表達；PTL在韌皮部和木質部7 中表達量最高，顯著高于葉片和木質部2；NAC86在韌皮部中表達水平最高，但在木質部2 和木質部7 中也有較高表達水平；MYB46、C3H14、CesA和LAC17的表達模式相對一致，在木質部2 的表達量最高，木質部7 中略低，而韌皮部中表達量極低。這些結果進一步證明了測序結果的可靠性，同時也表明這些基因可能在木質部生成中發(fā)揮重要作用。

圖5 木質部候選基因的組織表達模式Fig.5 The expression of candidate genes of xylem in various tissues

3 討論

為了研究種植密度影響林木徑向生長的分子機制，本研究聯(lián)合PacBio Iso-Seq 和RNA-Seq 測序技術，在高、低2 個種植密度下尾巨桉DH32-29幼苗直徑及直徑生長量開始出現(xiàn)顯著差異時進行轉錄組測序，獲得45 490 條在主莖表達的非冗余轉錄本，共有44 888（98.68%）個轉錄本得到注釋，其中，在NR 數(shù)據(jù)庫中有98.47%（44 854/45 490）的轉錄本能夠預測到同源序列，注釋率較高，可能是因為巨桉（Eucalyptus grandisHill）全基因組測序已經(jīng)完成，剩余602（1.32%）個轉錄本未得到注釋，可能是尾巨桉中的新基因。由于缺乏尾巨桉全基因組序列，本研究結合二代和三代轉錄組測序結果進行分析有利于數(shù)據(jù)的分類和定量研究，提高了序列的準確性和可靠性。

木材的生長性狀是由多基因決定的高度復雜性狀。調控因子在木材形成過程中發(fā)揮了重要的調控作用，改變調控因子的表達就可促使多個功能基因共同發(fā)揮作用，從而達到性狀改良的效果[10,23]。因此，本研究重點關注木質部中響應種植密度的60 個差異表達調控因子。高種植密度植株木質部中優(yōu)勢表達的調控因子多數(shù)與脅迫相關，這與歐美楊（Populus×euramericana(Dode) Guinier）中參與脅迫響應的轉錄因子在高種植密度植株葉片中上調表達相一致[24]；而低種植密度中優(yōu)勢表達的調控因子主要與細胞分裂活動和細胞壁物質合成相關。這可能是因為在高種植密度條件下，植株生長過密，易受多種生物或非生物因子的脅迫，從而引起脅迫相關基因上調表達以響應劣勢的生長環(huán)境，而在低種植密度條件下有利于植株生長，促進了木質部生成相關基因的表達，因而使其徑向生長顯著大于高種植密度。

為了更好的挖掘和鑒定響應種植密度的木質部生成關鍵基因，本研究進一步對低種植密度植株木質部中優(yōu)勢表達的調控因子進行分析，發(fā)現(xiàn)了一些重要的調控因子編碼基因PXL2、MYB46、C3H14、NAC86和PTL，它們在低種植密度植株木質部中均上調表達（圖1A）。蛋白互作網(wǎng)絡和組織表達分析結果表明，類受體激酶PXL2與2 個關鍵節(jié)點基因CUL1和T15D22.7互作并具有一致的表達模式（圖3），這與擬南芥中的研究相吻合，PXL2與PXL1、PXY共同協(xié)作調控維管組織發(fā)育過程中細胞分裂活動[6]；此外，楊屬植物中PXL2的同源基因PXY通過調控維管組織的細胞分裂活動參與木質部生成[25]，由此表明，PXL2、CUL1和T15D22.7可能共同協(xié)作調控尾巨桉主莖維管組織的細胞分裂活動，這可能是造成不同種植密度下植株徑向生長差異的重要原因之一。

MYB 和鋅指蛋白C3H 是植物中普遍存在的轉錄因子，越來越多的研究表明，這些轉錄因子在植物次生壁的合成中發(fā)揮重要作用[26-31]。在擬南芥中，MYB46及其楊屬和桉屬植物中的同源基因PtrMYB21和EgMYB2是次生壁合成調控網(wǎng)絡的主控開關成員，能夠直接調控其它轉錄因子和次生壁合成相關基因的表達[26-30]；AtC3H14及其楊屬植物中同源基因PdC3H17/18是MYB46的直接靶基因，可正向調控次生壁的合成及木質部生成相關MYBs基因的表達[31]。本研究表明，MYB46、C3H14、CesA和參與木質素合成的漆酶基因LAC17等主要在木質部中表達，且均在直徑生長量顯著增大的植株上調表達，推測MYB46、C3H14 可能參與調控LACs和CesAs的表達，并且MYB46、C3H14 可能通過參與尾巨桉木質部生成過程中次生壁合成的調控對種植密度做出響應。這些結果表明木質部生成過程中轉錄調控的復雜性，進一步揭示這些調控因子在木質部生成中的作用具有重要意義。

值得關注的是，在擬南芥根系發(fā)育過程中，NAC86直接參與韌皮部內篩管分子成熟過程的調控[32]，PTL作為抑制子參與調控木質部細胞分化活動[5]，尚未見二者參與木本植物木質部生成調控的報道。在本研究中，NAC86、PTL在低種植密度植株的木質部中上調表達（圖1A），在韌皮部和木質部中均有表達（圖5），這表明擬南芥同源基因NAC86和PTL可能同時參與木本植物中木質部和韌皮部的生成，PTL可能在木本植物的維管組織發(fā)育過程中發(fā)揮正向調控作用，暗示這2 個基因可能具有與草本植物中不同的功能，值得進一步研究。

基于以上研究結果，筆者提出了關于種植密度影響尾巨桉徑向生長的分子調控網(wǎng)絡模型（圖6）。該模型中，種植密度通過影響調控維管組織細胞分裂的PXL2、調控次生壁合成的MYB46和C3H14、以及調控木質部細胞分化的PTL及其下游功能基因的表達模式進而促進了植株的徑向生長。然而，該模型中各部分成員尤其是PTL、NAC86的具體功能和作用機制仍有待深入研究。為此，接下來將進一步開展基因功能研究，明確模型中關鍵成員如PTL、NAC86和PXL2在木材形成中的功能，并開展調控機制研究，明確其上下游調控基因，同時將研究結果推廣到更多研究群體中，驗證研究結論普遍性，以使研究結果得以為生產(chǎn)實踐應用。

圖6 種植密度影響尾巨桉徑向生長的分子調控網(wǎng)絡模型Fig.6 The regulatory network model of radial growth proposed for the E.urophylla×E.grandis related to planting density

4 結論

本研究從全基因組轉錄水平和基因表達兩方面，鑒定了參與種植密度響應的顯著表達的基因，獲得了參與桉樹徑向生長的候選調控因子，提出了種植密度影響尾巨桉徑向生長的分子調控網(wǎng)絡模型，并提出NAC86和PTL在木本植物徑向生長與草本植物生長中的功能具有差異，該研究對于揭示桉樹人工林密度調控的分子機制和木質部生成關鍵基因的研究具有重要意義。