毛竹CesA基因鑒定及其表達調控研究

王新悅 朱成磊 楊克彬 李廣柱 李 真 袁婷婷 高志民

(國際竹藤中心竹藤資源基因科學與基因產業(yè)化研究所，國家林業(yè)和草原局/北京市共建竹藤科學與技術重點實驗室，北京 100102)

纖維素是植物細胞壁的主要成分，是十分重要的可再生資源[1]。纖維素合成酶(cellulosesynthase, CesA)對纖維素合成起關鍵性作用，通常以CesA復合體(cellulose synthase complex, CSC)的形式發(fā)揮功能，該復合體由6個三聚體形成，而三聚體由多個CesA組成[2]，CSC在植物的纖維素合成與沉積過程中發(fā)揮重要作用[3]。研究表明，CesA包含多個保守區(qū)域，其中D、DxD、D和QxxRW是該家族的特征結構域[4-5]。不同物種中CSC發(fā)揮功能需要多個CesA[6]，自第一個CesA基因在棉花(Gossypiumspp.)中被確定后[7]，已在擬南芥(Arabidopsisthaliana)[8]、水稻(Oryzasativa)[4]和二穗短柄草(Brachypodiumdistachyon)[9]中分別鑒定出10、11和10個CesA基因。研究發(fā)現(xiàn)，擬南芥中AtCesA1/3/6與初生細胞壁合成相關[10]，AtCesA4/7/8與次生細胞壁相關[11]；水稻中OsCesA1/3/6與初生壁相關，OsCesA4/7/9與次生壁相關[4]，而且發(fā)現(xiàn)赤霉素可以通過信號轉導與NAC29/31 結合調控OsMYB61轉錄因子，OsMYB61與OsCesA4/7/9上游啟動子結合，進而聯(lián)級調控次生壁纖維素的合成[12]。

竹子具有生長快、材性好等特點，在緩解我國木材供給壓力方面發(fā)揮著重要作用。纖維素是竹材的主要成分，CesA基因的表達水平直接影響到纖維素的含量和品質，從而影響竹材材性[13]。因此，研究CesA基因的分子特征及其表達調控對竹子的開發(fā)利用具有重要價值。初步研究表明，毛竹(Phyllostachysedulis)各組織纖維素的含量與CesA基因表達量有一定的關聯(lián)[14]，筍中纖維素含量隨基因表達量上升而增加[15]，并預測PeCesA6/9/13對初生壁和次生壁的形成均起到重要作用，而PeCesA1/2/3/4/5主要參與次生壁的合成[16]。另有研究發(fā)現(xiàn)巨龍竹(Dendrocalamussinicus)中2個CesA基因與初生壁纖維素生物合成有關[17]。雖然在毛竹基因組草圖中鑒定了19個PeCesAs[18]，并對其中16個進行了定量分析[19]，但其功能尚未被深入研究，尤其在表達調控方面更是空白。染色體水平毛竹基因組的發(fā)布[20]為更全面研究其中的CesA提供了組學數(shù)據(jù)。本研究在對毛竹CesA基因家族成員全基因組鑒定的基礎上，通過基因表達模式篩選出與次生壁纖維素合成相關的CesA基因及MYB轉錄因子，并用實時熒光定量PCR(quantitative real-time PCR，qRT-PCR)和酵母單雜交技術進行驗證，以期為揭示竹子纖維素生物合成與調控分子機制提供理論參考。

1 材料與方法

1.1 試驗材料

以江西省林科院實驗林場的毛竹為對象，采集不同高度筍(0.5、1.0、2.0、4.0、6.0和8.0 m)的代表性節(jié)間(第15節(jié))，每種樣品平均隨機取3株，液氮速凍后存于-80℃冰箱，用于研究CesA與MYB在不同高度筍中的表達模式。

1.2 試驗方法

1.2.1 毛竹CesA基因家族成員鑒定及理化性質分析 通過PLAZA數(shù)據(jù)庫(https://bioinformatics.psb.ugent.be/plaza/)下載水稻、二穗短柄草和擬南芥的CesA氨基酸序列，利用TBtools基于毛竹基因組數(shù)據(jù)[21]查找CesA同源序列，在Pfam(http://pfam.xfam.org/)網站利用CesA蛋白查找保守結構域[22]，再運用Smart BLAST(https://blast.ncbi.nlm.nih.gov/smartblast/?LINK_LOC=BlastHomeLink#)進一步對候選基因進行逐條比對分析，保留保守結構域完整的基因，根據(jù)上述基因在Scaffold上的位置順序進行命名。采用在線工具Softberry(http://linux1.softberry.com/all.htm)分析CesA的基本理化性質。

1.2.2 毛竹CesA基因的生物信息學分析 通過TBtools中BLASTP程序對毛竹氨基酸序列進行雙向比對，閾值E<10-5，獲得共線性基因對[23]，通過TBtools內置的MCScanX程序對毛竹進行物種內的CesA共線性結果進行可視化展示[24]。并計算基因對之間的同義(Ks)和非同義(Ka)核苷酸替換率[25]。用TBtools軟件分析毛竹CesA的基因結構；用在線軟件MEME(http://meme.nbcr.-net/meme/intro.html)分析毛竹CesA的蛋白保守基序，利用Plant-mPLoc(http://www.csbio.sjtu.edu.cn/bioinf/plant-multi/#)在線軟件預測毛竹CesA的亞細胞定位。利用MEGA7.0內置的MUSCLE程序對毛竹、水稻、擬南芥和二穗短柄草的CesA氨基酸序列進行多序列比對，并采用鄰接法(neighbor-joining, NJ)構建系統(tǒng)進化樹(使用系統(tǒng)默認值，其中校驗參數(shù)Bootstrap值設置為重復1 000 次)[26]，模型為Poisson model[27]，序列中缺失位點選擇完全刪除，根據(jù)擬南芥[10-11]與水稻[4]中報道次生壁與初生壁相關的基因對進化樹的分支進行分類。

1.2.3 毛竹CesA基因表達模式分析 利用NCBI下載的毛竹26個組織/發(fā)育階段(鞭、筍、根、葉片、葉鞘和筍芽)的轉錄組數(shù)據(jù)(登錄號：SRX2408703～SRX2408728)[20]，篩選其中CesA基因的FPKM值，并取對數(shù)Log2(FPKM+1)作為基因的表達量，利用TBtools軟件中最長距離法與層次聚類法繪制表達譜熱圖，并利用基迪奧生物網站(https://www.omicshare.com/tools/) spearman系數(shù)繪制相關性分析圖。利用BambooNET數(shù)據(jù)庫(http://bioinformatics.cau.edu.cn/bamboo/)中Mutual rank與Pearson correlation coefficient(PCC)系數(shù)查找與毛竹CesA基因具有共表達關系的基因，于基迪奧生物網站繪制共表達網絡圖。

1.2.4 毛竹CesA基因在不同高度筍中的表達模式分析 根據(jù)毛竹CesA基因序列設計的特異性定量引物，由北京博邁德基因技術有限公司合成(表1)。采用TRIzol法(Promega, 美國, Z3101)提取毛竹各樣品的總RNA，用DNase Ⅰ (TaKaRa, 日本, 2270A)除去基因組DNA后，使用反轉錄試劑盒(Promega, 美國, A5003)合成cDNA。使用qTOWER熒光定量PCR儀(Analytikjena, 德國, qTOWER 2.2)，體系參照Roche Light Cycler?480 SYBR Green 1 Master kit試劑盒(Roche, 中國, 4707516001)，程序：95℃預變性6 min；95℃變性10 s，62℃退火10 s，45個循環(huán)[28]。以PeTIP41為內參基因[29]，3次生物學重復，采用2-ΔΔCT法[30]分析基因的相對表達量，利用IBM SPSS 26.0對3個技術重復的均值和標準差進行繪圖。

1.2.5 酵母單雜交分析 根據(jù)PeMYB35(PH02Gene26889)的蛋白質編碼區(qū)(coding sequence, CDS)序列和PeCesA1啟動子的PeCesA1pro序列，設計添加載體構建所需酶切位點的引物，進行PCR擴增。反應體系(20 μL)如下：PCR Buffer 2.0 μL，5′寡核苷酸 9.0 μL，3′寡核苷酸 9.0 μL。反應條件：95℃ 30 s，72℃ 2 min，37℃ 2 min，25℃ 2 min，1個循環(huán)。分別使用SmaI酶、EcoRI酶和SacI酶進行酶切，將PeMYB35的CDS序列和PeCesA1pro分別連接到表達載體pGADT7-Rec2 和pHis2的多克隆位點，將測序正確的質粒分別命名為pGADT7-Rec2-MYB35和pHis2-PeCesA1pro。另外，根據(jù)PeCesA1pro中3個GAMYB元件(CAACCTAC、CAACCCAA和CAACCAAC)分別設計引物進行擴增；同時利用EcoRI和SacI酶對pHis2載體進行酶切，對PCR產物與膠回收產物進行連接。連接體系：10×T4 Ligase 1.0 μL，T4 DNA連接酶 1.0 μL，pHis2 4.0 μL，元件PCR產物 4.0 μL，16℃過夜連接。將測序正確的質粒分別命名為pHis2-Y1、pHis2-Y2和pHis2-Y3。

將構建好的pHis2-PeCesA1pro、pHis2-Y1、pHis2-Y2 和pHis2-Y3質粒分別與pGADT7-Rec2-MYB35共轉化到酵母感受態(tài)細胞(Y187)中，同時設置陽性(p53His2+pGADT7-Rec2-53)和陰性(p53His2+pGADT7-Rec2-MYB35)對照，將共轉化質粒的酵母菌液分別涂布在DDO (-Leu/-Trp)篩選培養(yǎng)基上培養(yǎng)，篩選后在含有3-AT (20 mmol·L-1)的TDO (-His/-Leu/-Trp) 平板培養(yǎng)2～3 d[31]。

2 結果與分析

2.1 毛竹CesA基因家族成員的鑒定與進化分析

利用擬南芥和水稻的同源基因序列，在毛竹基因組數(shù)據(jù)庫中進行BLAST分析，共鑒定獲得21個CesA基因，分布在11個基本骨架(Scaffold)上，其中Scaffold14上有4個基因，Scaffold4和Scaffold10上有3個基因，其他Scaffold上有1～2個基因，根據(jù)其在Scaffold上順序分別命名為PeCesA1～PeCesA21(圖1)。PeCesAs編碼蛋白氨基酸長度為904～1 093 aa，預測分子量為101.10～123.63 kDa之間，理論等電點在5.95～8.60之間，大多數(shù)呈中性或堿性蛋白；亞細胞定位預測結果表明，21個PeCesAs均定位于質膜(表2)，說明PeCesAs可能在質膜上行使功能。

圖1 PeCesAs的Scaffold定位和共線性分析Fig.1 The scaffold location and collinearity analysis of PeCesAs

共線性分析發(fā)現(xiàn)，有19個PeCesAs之間存在共線性關系，且所有基因對的非同義對同義取代比(Ka/Ks)均小于1.0，說明PeCesAs基因在進化過程中可能經歷了較強的純化選擇[32-33]。利用毛竹、水稻、擬南芥和二穗短柄草的52個CesA的氨基酸序列構建的系統(tǒng)發(fā)育樹表明，所有CesAs被聚類在7個分支，在各分支中多數(shù)PeCesAs與二穗短柄草及水稻的CesA聚類在一起，與擬南芥的CesA較遠(圖2)。另外，研究發(fā)現(xiàn)水稻中OsCesA4/7/9以及擬南芥中AtCesA4/7/8主要參與次生細胞壁纖維素的合成，PeCesA1/5/10/11/17與OsCesA4/7/9和AtCesA4/7/8聚類在較近的分支，因此推測上述編碼蛋白在毛竹次生細胞壁纖維素合成中可能具有相似的功能。

注：Pe、Os、At和Bd分別代表毛竹、水稻、擬南芥和二穗短柄草。Note: Pe, Os, At and Bd indicate Phyllostachys edulis, Oryza sativa, Arabidopsis thaliana and Brachypodium distachyon, respectively.圖2 基于不同植物CesA家族成員構建的系統(tǒng)發(fā)育樹Fig.2 Phylogenetic tree constructed on the base of CesA family members in different plant species

2.2 PeCesAs基因結構及其編碼蛋白的保守基序

基因結構分析表明，21個PeCesAs具有8～14個內含子(圖3)，9～15個外顯子，其中具有14個外顯子的基因有11個，13個外顯子的基因有5個，15與9個外顯子的基因分別有2個，12個外顯子的基因只有1個(PeCesA1)，PeCesA18比PeCesA14多1個外顯子與1個內含子。各基因之間的差異除了外顯子的數(shù)量外，還表現(xiàn)在內含子的位置差異，由此推測這些差異可能會導致基因功能不同。

圖3 PeCesAs基因結構分析Fig.3 Gene structure analysis of PeCesAs

蛋白保守基序分析表明，在PeCesAs中共鑒定出20個保守基序，其中4個纖維素合成酶家族特征性保守結構域(DD、DCD、ED和QVLRW)分別分布在Motif 6、Motif 10、Motif 13和Motif 14中[34]。Motif 1、Motif 2～Motif 15為所有成員共有的保守基序；Motif 20為PeCesA1/4/6/7/14/15/18/19/21共同缺失的保守基序；此外，PeCesA7/14還缺少Motif 16～Motif 20；PeCesA16缺少Motif 19 (圖4)。

圖4 PeCesAs保守基序分析Fig.4 Conserved motif analysis of PeCesAs

2.3 PeCesAs組織特異性表達分析

基于毛竹26個組織的轉錄組數(shù)據(jù)分析表明，21個PeCesAs在毛竹不同組織以及同一組織不同生長時期的表達模式均存在著一定差異(圖5)。所有PeCesAs在各個組織中均檢測到表達，每個組織中各成員的表達量各不相同，說明其可能不同程度地參與毛竹不同組織的纖維素合成。表達量相似的基因被聚類在同一分支，如PeCesA3/6/7/8/9/19在毛竹根和筍的中部和基部表達量均較高，在葉片和筍的尖部表達量極低；而PeCesA12/13在所有組織中表達量均極低；PeCesA2/14/16/18/20在大部分組織中均表達，其中5個PeCesAs (PeCesA1/5/10/11/17)在筍的中部和基部的表達量隨著筍高度增加而上升，在3.0 m筍中的相對表達量顯著高于其他高度的筍。因此，推測這些PeCesAs在毛竹筍的纖維素合成過程中發(fā)揮著更為重要的作用。

圖5 PeCesAs在不同組織中的表達分析Fig.5 Expression analysis of PeCesAs in different tissues

2.4 PeCesAs與PeMYBs相關分析

通過共表達網絡預測發(fā)現(xiàn)，PeMYB35 (PH02Gene26889)和PeMYB37 (PH02Gene12586)與5個PeCesAs (PeCesA1/5/10/11/17)存在較強的共表達關系(圖6-A)。利用毛竹不同高度筍轉錄組數(shù)據(jù)分析PeMYBs與PeCesAs表達量的相關性，發(fā)現(xiàn)2個PeMYBs與5個PeCesAs均具有極強相關性(PCC>0.8) (圖6-B)，其中PeMYB35較PeMYB37與PeCesAs的相關性更強。對PeCesAs的上游啟動子分析發(fā)現(xiàn)，PeCesA1/5/10/11/17的啟動子中均包含GAMYB (CAACCNNN)元件(圖6-C)。在水稻中已證明MYB轉錄因子能夠與GAMYB結合，調控CesA的表達[12]。因此，推測PeMYB35和PeMYB37可能通過與PeCesA1/5/10/11/17啟動子中的GAMYB結合來調控毛竹纖維素的合成。

注：A：PeCesAs與PeMYB35共表達網絡預測；B：PeCesAs與PeMYB35在不同高度筍中表達量相關性分析；C：PeCesAs啟動子中的GAMYB元件分析。Note: A: Co-expression network prediction of PeCesAs and PeMYB35. B: Correlation analysis of PeCesAs and PeMYB35 expression in different height shoots. C: GAMYB in the promoters of PeCesAs.圖6 PeCesAs與PeMYB35關系圖Fig.6 The relationship diagram of PeCesAs and PeMYB35

2.5 PeMYB35與PeCesAs調控關系驗證

為驗證PeMYB35與5個PeCesAs的調控關系，采用qRT-PCR方法檢測PeMYB35和PeCesA1/5/10/11/17在不同高度筍中的表達模式。圖7結果顯示，隨著筍高度的增加，PeMYB35和5個PeCesAs基本呈現(xiàn)相似表達趨勢，即不同程度地先下降后上升再下降，該趨勢與圖5相似，在圖7-A中均在6.0 m高筍中的表達量達到最高值，之后有所下降。另外，酵母單雜交結果表明，陽性對照(p53His2+pGADT7-Rec2-53)和pHis2-PeCesA1pro+pGADT7-Rec2-MYB35、pHis2-Y1+pGADT7-Rec2-MYB35、pHis2-Y2+pGADT7-Rec2-MYB35和pHis2-Y3+pGADT7-Rec2-MYB35能夠在TDO/3-AT培養(yǎng)基上正常生長，而陰性對照(p53His2+pGADT7-Rec2-MYB35)不能正常生長(圖7-B)。由此說明，PeMYB35能夠與PeCesA1啟動子中3個不同的GAMYB相結合，進而調控基因的表達。

注：A：筍中PeMYB35和PeCesAs的表達定量分析；*和**分別表示其他高度筍中基因相對表達量與0.5 m筍中的差異顯著性水平為P<0.05和P<0.01。B：PeMYB35與PeCesA1pro共轉化酵母。Note: A: Expression analysis of PeMYB35 and PeCesAs in shoots with different heights. *, ** indicate significant differences at 0.05 and 0.01 level respectively between the relative expression in higher shoots and that in 0.5 m shoots. B: Co-transformation of yeast with PeMYB35 and PeCesA1pro.圖7 PeMYB35與PeCesAs調控關系驗證Fig.7 Verification of the regulatory relationship between PeMYB35 and PeCesAs

3 討論

纖維素含量嚴重影響著竹子的生長與竹材材性，且影響采后竹筍的品質[35]。CesA的表達是影響纖維素生物合成的關鍵因素。隨著越來越多物種中CesA基因家族的結構特征、基因功能驗證及轉錄調控被研究報道，竹子中CesA的研究倍受重視[36]。對毛竹中CesA基因的結構特征、表達模式，以及轉錄調控進行分析，對于揭示CesA基因在竹材形成過程中的作用具有重要的參考價值。本研究從毛竹中鑒定出21個CesA家族成員，與Peng等[18]的鑒定結果相似，其具有的保守結構域是纖維素合成的催化位點[37]。CesA家族基因的數(shù)量是擬南芥、水稻等二倍體物種的2倍左右，這可能與毛竹為異源四倍體有關[38]。本研究系統(tǒng)進化分析發(fā)現(xiàn)，PeCesAs聚類在7個分支上，與共線性分析結果相一致，呈現(xiàn)共線性的蛋白被聚類在同一個分支上，多數(shù)PeCesAs與水稻和二穗短柄草的CesA聚類到較近的分支，這與毛竹、水稻和二穗短柄草親緣關系較近的前人研究結果相一致[28]。

纖維素是細胞壁的重要組成部分，21個PeCesAs在毛竹26個不同組織中均檢測到表達，在筍和根中表達量較高，該結果與前人研究結果一致[18]，但本研究在芽和鞭中也檢測到了表達，這可能是取樣與轉錄組數(shù)據(jù)的差異造成的。此外CesA基因的表達具有傾向性，如白花泡桐中PfCesA4在莖中的表達量最高[39]，CiCesA3則在中間錦雞兒的葉中表達量最高[40]，毛白楊的PtoCesA4/7/8/17/18在木質部中高豐度表達[41]。不同的PeCesAs在不同組織中的表達豐度差異明顯，表明各基因成員的功能存在一定的差異。除了PeCesA12和PeCesA13在所有組織中表達較低外，其他多數(shù)PeCesAs均在次生壁木質化程度較高的鞭、根和筍(1.5、3.0和6.0 m)中表達豐度較高，表明這些PeCesAs的表達與毛竹細胞次生壁形成相關，與前人研究結果一致[4]。RNA-Seq數(shù)據(jù)分析結合qRT-PCR 發(fā)現(xiàn)PeCesA1/5/10/11/17主要在1.5、3.0和6.0 m筍中表達，且分別是OsCesA4/7/9[4]和AtCesA4/7/8[11]的同源基因，推測該基因可能在毛竹筍纖維素合成過程中發(fā)揮著重要作用，其具體生物學功能尚需開展深入研究。

纖維素的生物合成調控是一個復雜的網絡，涉及多種轉錄因子，其中MYB轉錄因子在調控次生壁纖維素合成中發(fā)揮著重要作用。在擬南芥中，AtMYB46通過與AtCesA4/7/8的啟動子中的M46RE元件結合，激活AtCesA4/7/8的表達，從而使纖維素含量升高[42]。在水稻中OsMYB58/63能與OsCesA7啟動子中AC-II 和SMRE3元件結合，上調基因的表達[43]，而NAC29/31 則通過調控MYB61，MYB61再與OsCesA4/7/9 中上游啟動子的GAMYB元件結合調控CesA的表達，從而實現(xiàn)對纖維素生物合成的調控，而SLR1和NAC29/31相互作用則會阻斷這條調控途徑[12]。在毛竹中PeMYB35和PeMYB37分別與擬南芥和水稻中的AtMYB46與OsMYB61為同源基因[44]。本研究發(fā)現(xiàn)并證明PeMYB35和PeMYB37均與5個PeCesAs (PeCesA1/5/10/11/17)之間存在共表達關系，在筍中的表達趨勢相似，表達量具有較強的相關性；酵母單雜交結果證實轉錄因子PeMYB35能夠與PeCesA1的上游啟動子中的MYB元件結合發(fā)生相互作用。由此推測，PeMYB35和PeMYB37可能通過調控PeCesAs (PeCesA1/5/10/11/17)的表達來參與毛竹纖維素的生物合成，但具體的調控分子機制以及對纖維素含量的影響尚需要通過進一步研究驗證。

4 結論

本研究從毛竹中鑒定了21個CesA基因家族成員(PeCesA1～PeCesA21)，轉錄組表達譜分析揭示出不同PeCesAs在毛竹不同組織與筍的不同發(fā)育時期的表達存在明顯差異，表明上述基因功能具有一定的差異性。不同高度筍的轉錄組表達譜和qRT-PCR結果均證實PeMYB35與PeCesA1/5/10/11/17呈現(xiàn)相似的表達趨勢，相關性極強(PCC>0.8)，且發(fā)現(xiàn)PeCesA1/5/10/11/17的上游啟動子均含有MYB結合位點GAMYB，表明MYB可能通過與PeCesAs啟動子的結合來調控其表達，PeMYB35與PeCesA1pro共轉化酵母的結果證明了這一推測。