轉錄組學對轉基因楊樹表型變化的代謝調控機制解析

和成成，吳照群，許 楠，張秀省，杜 鑫，楊海靈

(北京林業大學 生物科學與技術學院，北京 100083)

細胞分裂是所有細胞生物的基本特征之一，在動、植物和微生物生長發育過程中發揮著極其重要的作用。對細胞周期的研究一直是科學研究的熱點，它對動、植物的生長發育以及人類疾病研究都至關重要。與動物不同，大部分植物屬于胚后發育生物，這說明植物器官的形成與細胞分裂、分化以及細胞生長都聯系緊密。細胞周期的正常進行需要許多周期調控因子的共同參與。其中，Cyclin作為一種隨細胞周期的進行而不斷合成和降解的蛋白，在維持細胞周期正常進行的過程中發揮十分重要的作用。一般認為，D型細胞周期蛋白(CYCD)可以與CDKA形成復合物，控制著細胞周期G1/S的過渡[1]。在擬南芥中，共鑒定出10個D型Cyclin[2]；而楊樹中鑒定出24個D型Cyclin[3]，其中D2/4亞類在擬南芥和楊樹中分別有3個和2個成員[3-4]。這表明木本植物中的D2/4類Cyclin在進化過程中更加保守。

在草本植物的研究中，已經發現CYCD對植物生長發育至關重要。葉片受傷或邊緣被去除會促進CYCD;1的表達，通過與CDKA相互作用來增強CDKA激酶活性，進而造成葉片邊緣細胞重新進入細胞周期[5]。ICK2/KRP2可與CYCD2;1相互作用，ICK2/KRP2通過保持CYCD2;1失活來抑制根分枝，調節中柱鞘細胞對生長素波動的響應[6]。擬南芥KRP2抑制CDK活性，過表達KRP2導致葉片呈現鋸齒狀，葉片表皮細胞和葉肉細胞增大，同時抑制葉片核內周期[7]。擬南芥CYCD3參與調控維管發育，cycd3;1-3缺失突變體植物的莖和下胚軸明顯變細，形成層活性降低，次生生長受阻[8]。過表達AtCYCD3;1使得葉片無法形成正常的柵欄組織和海綿組織，同時葉表皮細胞由大量不完全分化的小的多邊形細胞組成[9]。CYCD3的缺失會損害芽分生組織功能并導致對細胞分裂素的應答降低[10]。AtCYCD2;1的轉基因品系提高了CDK激酶活性，縮短了G1期持續時間并減小了在S期和分裂階段的細胞大小。轉基因品系的葉細胞較小且數量更多，表皮細胞和柵欄組織細胞大小分別為對照的59%和69%，而葉片則是正常大小[11]。在木本植物中，過表達PtaCYCD1;2會造成楊樹產出波紋狀表型，葉片和莖明顯卷曲，且維管組織所占的比例增大[12]。過表達PtoCYCD3;3造成楊樹葉片變大且褶皺，莖增粗，同時出現分枝[3]。無論是在草本還是木本植物中，CYCD的表達量變化(缺失和過表達)會造成植株出現一系列的表型變化，如葉片形態、株高和莖粗等。

課題組前期研究發現，OE-PtoCYCD2;1造成楊樹生長減緩，株高降低，葉片向遠軸面卷曲且莖粗變細[13]。本研究主要以OE-PtoCYCD2;1楊樹為實驗材料，通過轉錄組學測序來研究PtoCYCD2;1對楊樹生長發育的影響，為研究木本植物D型細胞周期蛋白功能提供理論基礎。

1 材料和方法

1.1 實驗材料

選取1年生的經過qRT-PCR鑒定過的野生型(WT)和過表達(OE-PtoCYCD2;1) ‘741楊’(OE)進行扦插擴繁[13]。扦插后的植株生長于溫室中，栽培條件為20～25 ℃，光照強度3 500 Lx，16 h光照/8 h黑暗。選取生長3個月的幼苗幼葉進行轉錄組測序，WT和OE-PtoCYCD2;1均設置3個生物學重復，共6個樣本。將幼葉取下后迅速置于液氮中速凍，-80 ℃保存。

1.2 方 法

1.2.1 總RNA提取和文庫構建轉錄組測序委托北京諾禾致源科技股份有限公司完成。采用標準提取方法從組織中提取RNA，隨后對RNA樣品進行嚴格質控，質控方式主要是通過Agilent 2100 bioanalyzer精確檢測RNA完整性。

建庫起始RNA為總RNA，通過Oligo(dT)磁珠富集帶有polyA 尾的mRNA，隨后在Fragmentation Buffer中用二價陽離子將得到的mRNA 隨機打斷。以片段化的mRNA 為模版，隨機寡核苷酸為引物，在M-MuLV 逆轉錄酶體系中合成cDNA 第一條鏈，隨后用RNaseH 降解RNA 鏈，并在DNA 聚合酶 I 體系下，以dNTPs 為原料合成cDNA 第二條鏈。純化后的雙鏈cDNA經過末端修復、加A 尾并連接測序接頭，用AMPure XP beads 篩選370～420 bp 左右的cDNA，進行PCR 擴增并再次使用AMPure XP beads 純化PCR 產物，最終獲得文庫。將樣本庫檢合格后的樣本進行Illumina測序。

1.2.2 測序結果的統計和功能注釋為了保證數據分析的質量及可靠性，需要對原始數據進行過濾。主要包括去除帶接頭(adapter) 的reads、去除含N(N 表示無法確定堿基信息)的reads、去除低質量reads(Qphred ≤20 的堿基數占整個read長度的50%以上的reads)。同時，對clean data 進行Q20、Q30和GC 含量計算。后續所有分析均是基于clean data進行的高質量分析。

選擇NCBI上的毛果楊(Populustrichocarpa)基因組作為參考。使用HISAT2 v2.0.5 構建參考基因組的索引，并使用HISAT2 v2.0.5 將配對末端clean reads 與參照基因組比對。采用StringTie(1.3.3b)進行新基因預測。采用subread軟件中的featureCounts工具用于計算比對上的reads映射到每個基因的讀數。同時根據基因的長度計算每個基因的FPKM。

1.2.3 數據的統計與分析使用DESeq2 軟件(1.20.0)進行差異表達基因分析。作為有生物學重復的樣本，分析過程以|log2FC| ≥ 0 且 FDR(Q值) <0.05 作為差異基因篩選條件。通過clusterProfiler 3.4.4軟件實現差異表達基因的GO、KEGG、PFAM和轉錄因子富集分析，其中轉錄因子富集分析使用PlantTFDB數據庫對轉錄因子進行預測分析。

1.2.4 qRT-PCR驗證RNA提取使用天根生化科技(北京)有限公司的RNAprep Pure植物總RNA提取試劑盒；cDNA合成使用天根生化科技(北京)有限公司的FastKing cDNA第一鏈合成試劑盒(去基因組)；qRT-PCR使用天根生化科技(北京)有限公司的SuperReal PreMix Plus (SYBR Green)試劑盒進行實驗。用Primer 5進行引物設計，引物序列見表1。使用2-ΔΔCT計算所選差異表達基因的相對表達水平。

1.2.5 生理指標測定植物葉片葉綠素含量測定采用分光光度法，浸提法提取葉綠素[14]，提取液溶液為丙酮∶無水乙醇∶水(4.5∶4.5∶1)；植物葉片可溶性糖含量檢測使用索萊寶生物科技有限公司的植物可溶性糖含量檢測試劑盒；植物葉片和莖段木質素含量檢測使用索萊寶生物科技有限公司的木質素含量檢測試劑盒；以上3種生理指標測定時均分別取3株WT和OE-PtoCYCD2;1植株作為生物學重復，并且測定3次作為技術重復。

表1 qRT-PCR引物序列

2 結果與分析

2.1 測序數據質量分析

所有樣品在經過llumina平臺測序后生成原始的FASTQ數據，對原始測序數據進行質控統計。由表2可知，所有樣品的clean reads均達到6G以上，Q20的百分比均在98%以上，Q30的百分比在94%以上，GC含量均在44%左右，這些數據表明測序數據良好，可以進行下一步分析研究。

2.2 差異表達基因分析

通過計算FPKM，6個樣本中共統計得到42 816個基因。將OE-PtoCYCD2;1的所有基因表達量與WT相比，以|log2FC| ≥ 0且FDR(Q值)<0.05為篩選標準進行差異表達基因篩選，共篩選得到1 269個差異表達基因，其中上調基因700個，下調表達基因569個(圖1)。

2.2.1 GO分析對統計到的差異表達基因進行GO功能富集分析，主要分為生物過程(BP)、分子功能(MF)和細胞組成(CC) 3大類，氣泡圖中展示了在三大類中富集最顯著的30個小類(圖2，A)。其中細胞壁生物發生、植物型細胞壁生物發生、細胞碳水化合物代謝過程、幾丁質的響應和植物型細胞壁組織或生物發生是差異基因富集最顯著的生物學過程；葡萄糖基轉移酶活性、UDP-葡萄糖基轉移酶活性、UDP-糖基轉移酶活性、天冬氨酸型內肽酶活性和天冬氨酸型肽酶活性是差異基因富集最顯著的分子功能；Cul4A-RING E3泛素連接酶復合物、Cul4-RING E3泛素連接酶復合物、光系統Ⅱ析氧復合體、胞外區和類囊體膜是差異基因富集最顯著的細胞組分。

2.2.2 KEGG分析對篩選到的差異基因進行KEGG富集分析，氣泡圖中展示了富集最顯著的20個代謝通路(圖2，B)，其中碳代謝、半胱氨酸和甲硫氨酸代謝、甘氨酸、絲氨酸和蘇氨酸代謝、光合生物中的碳固定、植物晝夜節律、氨基酸的生物合成、苯丙烷生物合成、乙醛酸和二羧酸代謝、磷酸肌醇代謝和淀粉和蔗糖代謝是富集最顯著的代謝通路。

2.2.3 Pfam分析對篩選到的差異基因進行Pfam富集分析，氣泡圖中展示了富集最顯著的10種結構域(圖3，A)。其中AP2結構域、N 端木聚糖酶抑制劑、C 端木聚糖酶抑制劑、真核天冬氨酰蛋白酶和糖基轉移酶家族富集最顯著。

通過使用PlantTFDB對轉錄因子進行預測分析，共預測出2 422個轉錄因子，其中包括145個差異表達基因，這些差異表達基因分別屬于31個轉錄因子基因家族。圖3，B中展示了數目最多的10種轉錄因子，其中數目最多的為ERF轉錄因子且均上調表達，這與結構域富集結果相同。植物AP2/ERF類轉錄因子參與植物生長發育過程[15]，這可能暗示OE-PtoCYCD2;1主要通過調控ERF轉錄因子來調控植物的生長發育。

表2 樣本測序數據質量匯總

2.2.4 轉錄組數據驗證為驗證轉錄組數據可靠性，選取9個差異表達基因進行qRT-PCR實驗驗證(圖4)。結果顯示，9個差異表達基因的qRT-PCR結果與RNA-seq測定的差異基因表達水平變化趨勢基本一致，表明了本研究所用RNA-seq結果具有極高的可靠性。

2.3 PtoCYCD2;1對楊樹的影響分析

2.3.1 木質素合成OE-PtoCYCD2;1與WT相比，明顯顯示出植株株高降低，葉片卷曲[13]。同時GO富集分析顯示，在植物細胞壁合成相關的生物學過程中的相關基因表達明顯發生改變。進一步分析發現，有8個差異表達基因富集到木質素合成通路中，且均下調表達(圖5)。其中有5個漆酶(Laccase)基因均下調表達，這可能表明PtoCYCD2;1主要通過影響LAC表達來影響植物木質素合成。這表明OE-PtoCYCD2;1會影響楊樹木質素合成途徑中相關基因表達，這可能是造成OE-PtoCYCD2;1植株株高降低，莖變細的原因。

2.3.2 卡爾文循環KEGG富集分析顯示，29個差異表達基因富集到碳代謝通路中，上調表達基因2個，下調表達基因27個。卡爾文循環作為植物中光合碳同化的基本途徑，在植物生長發育中至關重要。對碳代謝通路中富集的差異表達基因進一步分析發現，其中8個與卡爾文循環有關的被富集到，且均下調表達(圖6)。該通路中有3個果糖-1,6-二磷酸醛縮酶(fructose-1,6-bisphosphate aldolase,FBA)基因均下調表達，可能表明PtoCYCD2;1主要通過影響FBA表達來影響植物木質素合成。總之，OE-PtoCYCD2;1會影響楊樹卡爾文循環途徑中相關基因表達，對植物生長發育產生影響。這可能是造成OE-PtoCYCD2;1植株整體生物量降低的原因。

2.4 PtoCYCD2;1對楊樹生理指標的影響

通過觀察WT和OE-PtoCYCD2;1的生長過程，發現與WT相比，OE-PtoCYCD2;1的葉片顏色更綠，由此猜測可能由于葉片的葉綠素含量差異造成的。通過檢測WT和OE-PtoCYCD2;1幼葉和成熟葉的葉綠素含量，結果發現，幼葉和成熟葉的葉綠素a(Chla)分別增加52.98%和57.22%；幼葉和成熟葉的葉綠素b(Chlb)分別增加71.44%和138.16%，較Chla增加的比例更大，尤其是成熟葉；幼葉和成熟葉的總葉綠素含量分別增加57.36%和78.22%(表3)。結果說明：OE-PtoCYCD2;1會提高單位質量下葉片(無論是幼葉還是成熟葉)的葉綠素含量，尤其是Chlb的增加更加明顯。

植物體內可溶性糖含量通常作為植物碳代謝變化的指標之一。通過檢測WT和OE-PtoCYCD2;1的幼葉和成熟葉的可溶性糖含量，結果發現，與WT相比，OE-PtoCYCD2;1的幼葉可溶性糖含量變化不大，但成熟葉的可溶性糖含量下降了12.72%(表4)。結果表明，OE-PtoCYCD2;1會降低單位質量下葉片的可溶性糖含量，且對成熟葉的影響更大。

轉錄組結果表明，OE-PtoCYCD2;1植株中與木質素合成相關的基因下調表達，因此我們檢測了WT和OE-PtoCYCD2;1的葉片和莖段中木質素含量。結果發現，與WT相比，OE-PtoCYCD2;1幼葉中木質素含量下降了4.48%，而成熟葉中木質素含量下降了8.03%；OE-PtoCYCD2;1幼莖中木質素含量下降了20.03%，成熟莖中木質素含量下降最多，達到了31.63%(表4)。該結果表明，OE-PtoCYCD2;1會降低單位質量下葉片和莖部的木質素含量，在成熟莖中下降更加明顯。

表3 2種材料葉片葉綠素含量

表4 2種材料葉片可溶性糖和木質素含量

3 討 論

植物生長發育離不開細胞周期的正常進行。與動物相比，植物細胞具有全能性，單個植物細胞可以生長成為一個新的個體，這離不開細胞分裂和分化。細胞周期蛋白(Cyclin)、細胞周期依賴性激酶(CDK)和細胞周期依賴性激酶抑制劑(ICK)形成三元復合物[1]，控制著細胞周期的正常進行，其中D型細胞周期蛋白(CYCD)參與G1/S的轉換。本研究在轉錄組中鑒定出所有的CDK和ICK，但發現只有CDKG;4上調表達。在草本植物擬南芥中，CDKG亞家族中有兩個成員：CDKG1和CDKG2[16]。CDKG1蛋白在高溫下維持花粉育性和產量有重要作用。缺乏 CDKG1的植物表現出正常的營養生長，在高溫環境下生長不育，但在較低溫度下可育。相比之下，缺乏 CDKG2 的植物在高溫和低溫下都可以完全繁殖[17]。CDKG1 和 CDKG2 均已顯示介導下游基因的 AS。在營養組織中CDKG1 調節剪接因子 U2AF65A 的 AS，而在花粉中它調節 CalS5 mRNA 的 AS[18-19]。在楊樹中，CDKG基因家族擴張至5個成員[3]，它們之間的功能尚未有報道。OE-PtoCYCD2;1植株中PtoCDKG;4表達上調，它們之間的聯系和造成這種變化的原因仍需進一步研究。

CLR5為AP2/ERF轉錄因子家族蛋白的一個成員，它受生長素調控表達并可以通過正向調節 A 型 RR、OsRR1 來抑制細胞分裂素信號傳導，從而促進冠根啟動[20]。PtaERF1在韌皮部組織中表達，可能暗示PtaERF1在維管組織發育和功能中具有潛在作用[21]。BOLITA基因(屬于 AP2/ERF 家族)的轉基因植物的葉子小于非轉基因植物，葉片變小是由于細胞大小和細胞數量均減少所造成[22]。本研究發現有26個AP2/ERF轉錄因子均上調表達，表明PtoCYCD2;1可能主要通過影響AP2/ERF轉錄因子的表達來調控植物的生長發育。

木質素作為植物生長發育中的一種次生代謝物質，在植物的生長發育和抗性方面具有重要的生物學功能[23]。木質素是植物細胞壁的主要成分之一，是一種分子量大、組成和結構復雜的天然酚類聚合物。木質素代謝參與植物生長發育，干擾木質素生物合成會導致植物生長發育受抑制和畸形發育[24]。2個擬南芥CCR1敲除突變體的花序莖表現出木質素水平降低，且兩者都具有植株矮小和延遲衰老的表型[25]。本研究通過轉錄組分析發現8個木質素合成通路基因下調表達，且葉片和莖中單位質量下的木質素含量均下降，又因為OE-PtoCYCD2;1植株的葉片更小，莖段更細，表明OE植株的總體質量會有所下降[13]，導致植株總體木質素總量更少，這可能是造成OE植株株高降低、莖粗變細的原因。

擬南芥CYCD2的表達受到蔗糖的誘導[6,26]，這可能表明CYCD2的表達在植物碳代謝上存在一定的相關性。本研究發現，共有29個差異基因富集到碳代謝通路上，2個差異基因上調表達，27個差異基因下調表達。其中，一些重要的碳代謝相關基因，如3-磷酸甘油醛脫氫酶(GADPH)、磷酸甘油異構酶(TIM)、丙酮酸激酶(PK)、蘋果酸脫氫酶(MDH)等均下調表達。其中卡爾文循環通路有8個差異表達基因也下調表達。以植物體內的可溶性糖作為檢測植株碳代謝變化的指標，發現葉片的可溶性糖含量降低，又由于OE植株的葉片面積變小，整體植株的質量有所下降[13]，推測植株總體可溶性糖含量降低，與轉錄組分析結果相呼應。表明PtoCYCD2;1可能調控植物的碳代謝水平，降低植物對CO2的利用效率，造成過表達植株的總生物量下降。

綜上，根據OE-PtoCYCD2;1的轉錄組數據，結合已發表的轉基因表型和文獻分析，認為PtoCYCD2;1的過表達能上調多個AP2/ERF轉錄因子的表達，并且下調碳代謝通路和木質素通路中一些關鍵基因的表達，減少植株中蔗糖等碳源和木質素的含量，從而造成過表達植株表型改變，植株總體生物量下降。本研究對木本植物中細胞周期蛋白表達量改變引起植株表型變化方面的功能探討提供了很好的研究思路和方向。