轉(zhuǎn)錄組學對轉(zhuǎn)基因楊樹表型變化的代謝調(diào)控機制解析

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

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

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

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

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

1 材料和方法

1.1 實驗材料

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

1.2 方 法

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

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

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

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

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

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)；植物葉片可溶性糖含量檢測使用索萊寶生物科技有限公司的植物可溶性糖含量檢測試劑盒；植物葉片和莖段木質(zhì)素含量檢測使用索萊寶生物科技有限公司的木質(zhì)素含量檢測試劑盒；以上3種生理指標測定時均分別取3株WT和OE-PtoCYCD2;1植株作為生物學重復，并且測定3次作為技術(shù)重復。

表1 qRT-PCR引物序列

2 結(jié)果與分析

2.1 測序數(shù)據(jù)質(zhì)量分析

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

2.2 差異表達基因分析

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

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

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

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

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

表2 樣本測序數(shù)據(jù)質(zhì)量匯總

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

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

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

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

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

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

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

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

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

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

3 討 論

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

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

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

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

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