毛竹GRF基因家族全基因組鑒定與表達(dá)分析

阮詩(shī)雨，張智俊，陳家璐，馬瑞芳，朱豐曉，劉笑雨

（1. 浙江農(nóng)林大學(xué) 省部共建亞熱帶森林培育國(guó)家重點(diǎn)實(shí)驗(yàn)室，浙江 杭州 311300；2. 浙江農(nóng)林大學(xué) 林業(yè)與生物技術(shù)學(xué)院，浙江 杭州 311300）

GRF(general regulatory factor)蛋白質(zhì)最先由MOORE等[1]在牛腦中發(fā)現(xiàn)，并根據(jù)淀粉凝膠電泳上的遷移特性命名。GRF蛋白質(zhì)是一類高度保守的同源或異源的二聚體蛋白質(zhì)，具有多種功能，廣泛存在于真核生物中，如酵母Pichia guilliermondii、擬南芥Arabidopsis thaliana、水稻Oryza sativa、花生Arachis hypogaea等。已有研究[2]表明：GRF蛋白質(zhì)家族通過(guò)與磷酸化的靶蛋白質(zhì)相互作用參與植物信號(hào)傳導(dǎo)、細(xì)胞定位、轉(zhuǎn)錄調(diào)控和應(yīng)激反應(yīng)等多種重要生命活動(dòng)過(guò)程，在植物代謝調(diào)控和生物合成反應(yīng)中發(fā)揮著重要作用，如擬南芥GRF蛋白質(zhì)可以與感光系統(tǒng)中的蛋白質(zhì)相互作用調(diào)節(jié)根系生長(zhǎng)發(fā)育[3]；葡萄Vitis vinifera GRF蛋白質(zhì)參與冷熱應(yīng)激反應(yīng)[4]；木薯Manihot esculenta GRF蛋白質(zhì)主要分布在細(xì)胞質(zhì)中，作用于淀粉合成酶Ⅲ靶蛋白質(zhì)，對(duì)淀粉的合成起到負(fù)調(diào)控作用[2]；菊花Dendranthema morifolium GRF蛋白質(zhì)參與開(kāi)花和周期調(diào)控，鹽、冷等脅迫響應(yīng)過(guò)程[5]；動(dòng)物細(xì)胞中GRF蛋白質(zhì)還可通過(guò)調(diào)節(jié)細(xì)胞周期，影響細(xì)胞凋亡，參與多種信號(hào)通路等方式來(lái)調(diào)控腫瘤進(jìn)程[6]。GRF活化后可以使G2/M期阻滯從而起到負(fù)調(diào)控細(xì)胞周期，發(fā)揮抑制癌基因的作用[7]。在動(dòng)物中GRF蛋白質(zhì)的過(guò)表達(dá)可能轉(zhuǎn)化為一種致癌因子，促進(jìn)腫瘤的發(fā)生[8]，還可能與腫瘤細(xì)胞耐藥性有關(guān)[9]。毛竹Phyllostachys edulis用途廣泛，筍和葉具有食用、藥用價(jià)值；竹材多用于建筑制造、工藝品制作。毛竹林是一種重要的經(jīng)濟(jì)林，具有重要生態(tài)價(jià)值，其固碳作用機(jī)制在不同的生長(zhǎng)階段有所差異[10]。毛竹基因組草圖已公布，且大量轉(zhuǎn)錄組數(shù)據(jù)也可以從公共數(shù)據(jù)庫(kù)中獲取[11]。目前根據(jù)毛竹全基因組數(shù)據(jù)進(jìn)行基因家族分析已取得了一定的成果，如ZF-HD基因家族[12]、B3基因家族[13]、APX基因家族[14]等，也分析了毛竹快速生長(zhǎng)期的基因表達(dá)[15?16]。但對(duì)于毛竹GRF基因家族的全基因組數(shù)據(jù)分析尚未有相關(guān)報(bào)道。本研究通過(guò)毛竹公開(kāi)的相關(guān)測(cè)序結(jié)果，利用生物信息學(xué)的方法，從基因組及轉(zhuǎn)錄組數(shù)據(jù)入手，對(duì)毛竹GRF基因進(jìn)行全基因組的鑒定與表達(dá)分析，擬為進(jìn)一步明確GRF基因家族在毛竹重要生長(zhǎng)發(fā)育過(guò)程中的功能解析提供依據(jù)。

1 材料與方法

1.1 基因家族來(lái)源、鑒定及理化性質(zhì)分析

毛竹基因組序列、編碼序列(CDS)、蛋白質(zhì)序列和基因組GFF注釋文件均從以下站點(diǎn)ftp://parrot.genomics.cn/gigadb/pub/10.5524/100001_101000/100498/[12]下載。從 Pfam 數(shù)據(jù)庫(kù)[17]中下載隱馬可夫模型(HMM) PF00244.17的結(jié)構(gòu)域數(shù)據(jù)，并以此結(jié)構(gòu)域數(shù)據(jù)為種子模型，用HMMER[18]檢索本地毛竹蛋白質(zhì)數(shù)據(jù)庫(kù)。在Excel 2018中，將E-value設(shè)置為≤1E?20，對(duì)檢索結(jié)果排序整理，去除重復(fù)，獲得候選基因。進(jìn)一步從毛竹全基因組數(shù)據(jù)庫(kù)中提取得到GRF家族成員的基因、CDS、蛋白質(zhì)序列以及基因結(jié)構(gòu)和位置信息；利用在線工具ProtParam(https://web.expasy.org/protparam/)、ProtScale(https://web.expasy.org/protscale/)[19]以及SignalP 4.1[20]在線分析GRF家族各成員理化性質(zhì)等。

1.2 家族進(jìn)化分析

依據(jù)毛竹、擬南芥、水稻GRF家族成員蛋白質(zhì)序列，分別通過(guò)ClustalW多重比對(duì)，用MEGA 7.0軟件鄰位連接(neighbor-Joining, NJ)法構(gòu)建種內(nèi)和種間系統(tǒng)進(jìn)化樹(shù)，自檢值取1 000次抽樣[21]。

1.3 基因結(jié)構(gòu)、基序和保守結(jié)構(gòu)域預(yù)測(cè)

根據(jù)毛竹全基因組的GFF注釋文件基因位置信息，分析毛竹GRF家族的基因結(jié)構(gòu)并繪制基因結(jié)構(gòu)圖；利用在線網(wǎng)站 NCBI Conserve Domain(https://www.ncbi.nlm.nih.gov/cdd/)和 MEME(https://www.ncbi.nlm.nih.gov/cdd/)對(duì)GRF家族成員的保守結(jié)構(gòu)域(domain)和基序(motif)進(jìn)行預(yù)測(cè)[22]，并通過(guò)TBtools[23]將結(jié)果可視化。

1.4 啟動(dòng)子分析

提取毛竹GRF基因上游1 500 bp序列作為啟動(dòng)子序列信息，通過(guò)在線預(yù)測(cè)軟件PlantCare[24]預(yù)測(cè)毛竹GRF基因的順式作用元件，并整理預(yù)測(cè)結(jié)果，富集順式作用元件，利用TBtools上的Simple Biosequence viewer功能進(jìn)行可視化分析。

1.5 染色體定位及共線性分析

利用 MCScanX[25]獲取GRF家族種內(nèi)、種間共線性關(guān)系，并用 TBtools軟件 Amazing Super Circos[26]和Multipe Synteny Plot分別對(duì)種內(nèi)和種間的結(jié)果可視化。

1.6 基因表達(dá)分析

選取NCBI SRA數(shù)據(jù)庫(kù)中毛竹不同組織器官：根(登錄號(hào)為ERR105075、ERR105076)，花序(登錄號(hào)為 ERR105069、ERR105070、ERR105071)，葉 (登錄號(hào)為 ERR105067、ERR105068、ERR105075)，鞭(登錄號(hào)為ERR105073、ERR105074)和筍不同生長(zhǎng)高度：0.2 m(登錄號(hào)為SRR6131114、SRR131113、SRR6131115)，0.5 m(登錄號(hào)為 SRR131117、SRR6131118、SRR5710699)和 1.0 m(登錄號(hào)為 SRR5710701、SRR5710702、SRR5710697)的轉(zhuǎn)錄組數(shù)據(jù)，分別計(jì)算毛竹GRF基因的 TPM(transcripts per million reads)值表示基因的表達(dá)豐度。為方便統(tǒng)計(jì)，對(duì)每個(gè)表達(dá)數(shù)值取以2為底的對(duì)數(shù)(log2)，使用TBtools Amazing Heatmap繪制基因表達(dá)熱圖，用對(duì)數(shù)轉(zhuǎn)換預(yù)處理數(shù)據(jù)，再用正態(tài)標(biāo)準(zhǔn)化的方法處理數(shù)據(jù)。

1.7 蛋白質(zhì)三級(jí)結(jié)構(gòu)同源模建

利用SWISSMODEL(https://www.swissmodel.expasy.org/)在線軟件[27]預(yù)測(cè)GRF蛋白質(zhì)的3D結(jié)構(gòu)。模建結(jié)果使用SAVES v5.0(https://servicesn.mbi.ucla.edu/SAVES/)[19]進(jìn)行評(píng)估。

2 結(jié)果與分析

2.1 基因家族成員鑒定及理化特性分析

根據(jù)植物GRF隱馬可夫模型(PF00244.17)搜索毛竹相關(guān)基因組數(shù)據(jù)，獲得相關(guān)GRF家族成員，然后通過(guò)E-value(≤1E?20)篩選、保守結(jié)構(gòu)域、基序特征分析，去除相同轉(zhuǎn)錄本重復(fù)，最終篩選得到13個(gè)GRF家族成員(表1)。將獲得13個(gè)GRF家族成員按照其在scaffold的分布先后順序命名為PeGRF01～PeGRF13。進(jìn)一步對(duì)PeGRF作蛋白質(zhì)特性分析，13個(gè)GRF蛋白質(zhì)中長(zhǎng)度最短的為PeGRF10(256個(gè)氨基酸)，最長(zhǎng)的為PeGRF09(293個(gè)氨基酸)，平均長(zhǎng)度266.8個(gè)氨基酸；各GRF蛋白質(zhì)等電點(diǎn)最小的為4.70(PeGRF02)，最大的為5.29(PeGRF01)，平均等電點(diǎn)為4.82；各GRF蛋白質(zhì)分子量最小的為PeGRF04(28.65 kD)，最大的為 PeGRF09(32.41 kD)，平均分子量為 29.79 kD。

表1 毛竹 GRF基因及其蛋白質(zhì)理化特性Table 1 Characteristics of PeGRF family genes and their deduced proteins

2.2 GRF基因家族分類與進(jìn)化樹(shù)構(gòu)建

利用MEGA 7.0對(duì)13個(gè)毛竹GRF、14個(gè)擬南芥GRF和8個(gè)水稻GRF的氨基酸序列比對(duì)后，采用NJ法進(jìn)行系統(tǒng)聚類分析(圖1)，絕大部分毛竹基因家族成員和水稻處于同一分支，表明毛竹與水稻的進(jìn)化關(guān)系較近。

圖1 毛竹 (Pe)、擬南芥 (At)和水稻 (Os)GRF家族系統(tǒng)進(jìn)化樹(shù)分析Figure 1 Phylogentic analysis of GRF gene family from Phyllostachys edulis (Pe), Arabidopsis thaliana (At) and Oryza sativa (Os)

2.3 GRF家族基因結(jié)構(gòu)、基序及保守結(jié)構(gòu)域

對(duì)毛竹GRF基因結(jié)構(gòu)分析發(fā)現(xiàn)：內(nèi)含子數(shù)量存在差異，非ε組成員都包含4個(gè)外顯子和3個(gè)內(nèi)含子，它們?cè)谖恢蒙细叨缺Ｊ?。ε組成員都具有不同于非ε組的內(nèi)含子-外顯子結(jié)構(gòu)，具有2個(gè)額外的N-末端內(nèi)含子[21]。利用NCBI-CDD對(duì)毛竹GRF基因進(jìn)行保守結(jié)構(gòu)域分析，PeGRF蛋白質(zhì)均包含14/3/3結(jié)構(gòu)域，毛竹GRF基因家族14/3/3結(jié)構(gòu)域存在一定的保守性，但該結(jié)構(gòu)域的分布位置有一定分化。利用MEME在線工具對(duì)該基因家族的保守基序預(yù)測(cè)，基數(shù)設(shè)置為10，結(jié)果顯示(圖2)：Motif1～6在每個(gè)家族成員中均出現(xiàn)，屬于高度保守結(jié)構(gòu)，其余基序在家族成員中出現(xiàn)的頻率及所在位置均存在一定的差異。

圖2 GRF家族基序分布特征Figure 2 Motif distribution of GRF family gene from Ph. edulis

2.4 啟動(dòng)子特征

如圖3所示：篩選出的部分典型的順式調(diào)控元件，除核心啟動(dòng)子TATA-box(5個(gè))和CAAT-box(16個(gè))外，還有與激素相關(guān)的順式調(diào)控元件，包括與赤霉素相關(guān)的GARE-motif(5個(gè))、P-box(3個(gè))，與生長(zhǎng)素有關(guān)的AuxRR-core(3個(gè))、TGA-element(6個(gè))，與脫落酸有關(guān)的ABRE(42個(gè))，與水楊酸有關(guān)的TCA-element(5個(gè))；與外部條件有關(guān)的順式調(diào)控元件，包括參與低溫響應(yīng)的LTR(2個(gè))和光響應(yīng)的G-box(48個(gè))。推測(cè)毛竹GRF蛋白質(zhì)家族可能參與激素和非生物脅迫響應(yīng)，家族基因表達(dá)模式可能有所不同。

圖3 PeGRF 基因家族啟動(dòng)子的上游順式作用元件Figure 3 Upstream cis-acting elements of promotor from PeGRF gene family

2.5 染色體分布及共線性分析

利用毛竹基因組GFF注釋文件提取PeGRF在scaffold上的分布特征，結(jié)果顯示：毛竹GRF基因在scaffold上分布不均勻，不同的scaffold基因分布密度不同，scaffold7、14、16、18和21僅包含1個(gè)PeGRF，scaffold3、13、15和22上分別包含2個(gè)。

利用TBtools工具，將毛竹GRF基因種內(nèi)和種間的共線性關(guān)系進(jìn)行了可視化分析。從圖4A中可以看出：除PeGRF02、PeGRF03和PeGRF07不存在種內(nèi)共線性關(guān)系外，其余家族基因成員間均有顯著的共線性關(guān)系，說(shuō)明GRF基因家族存在基因復(fù)制現(xiàn)象，推測(cè)在進(jìn)化過(guò)程中GFR基因可能通過(guò)復(fù)制進(jìn)行家族成員數(shù)量的擴(kuò)張。但PeGRF不存在串聯(lián)重復(fù)基因。物種間的共線性關(guān)系是反映不同物種來(lái)源于同一個(gè)祖先的現(xiàn)象。從圖4B可以看出：毛竹與水稻的共線性關(guān)系要明顯多于擬南芥，這可能與水稻和毛竹同屬于禾本科Gramineae，進(jìn)化關(guān)系較近有關(guān)。

圖4 毛竹 PeGRF 家族染色體分布 (A)及共線性分析 (B)Figure 4 Chromosomal distribution of PeGRF genes in Ph. edulis (A) and their collinear relationships (B)

2.6 GRF家族基因表達(dá)模式

本研究基于毛竹RNA-Seq轉(zhuǎn)錄組數(shù)據(jù)，對(duì)毛竹不同組織(葉、花序、鞭及根)以及不同生長(zhǎng)高度(0.2、0.5、1.0 m)的毛竹筍中的GRF表達(dá)量繪制熱圖。由圖5可以看出：除PeGRF10，PeGRF09在不同組織和生長(zhǎng)高度保持較低的表達(dá)量外，其他成員均有較高的表達(dá)量。在毛竹不同組織中，根和花序的表達(dá)量相對(duì)于葉和鞭要稍高；非ε組的GRF基因均有較高的表達(dá)。在竹筍的不同生長(zhǎng)階段，非ε組的GRF基因保持較高的表達(dá)水平；ε組不同的基因表達(dá)量有增有減，如PeGRF05在竹筍生長(zhǎng)各個(gè)階段均有較高的表達(dá)量，且隨生長(zhǎng)進(jìn)程表達(dá)量不斷增高；PeGRF06表達(dá)量隨生長(zhǎng)進(jìn)程呈下降趨勢(shì)。推測(cè)不同家族成員在參與組織器官發(fā)育的過(guò)程中發(fā)揮不同的作用，但其中的內(nèi)在分子機(jī)制還值得進(jìn)一步研究。

圖5 毛竹 GRF 基因家族表達(dá)水平熱圖分析Figure 5 Heatmaps of expression level of PeGRF family genes in Ph. edulis

2.7 GRF家族蛋白質(zhì)空間三級(jí)結(jié)構(gòu)

由圖6所示：毛竹GRF蛋白質(zhì)由2個(gè)單體連接而成，每個(gè)單體由反向平行的9個(gè)α螺旋組成，每個(gè)單體都存在與配體(FSC3、FEC4)相互作用的結(jié)合位點(diǎn)，2個(gè)FSC配體均與殼梭孢素有關(guān)，單體間構(gòu)成同源或異源二聚體，總體呈“W”型[28?29]。

圖6 毛竹 GRF家族蛋白質(zhì) SWISSMODEL同源模建的三維空間結(jié)構(gòu)Figure 6 Predicted 3D protein structure of the GRF family from Ph.edulis by SWISSMODEL

3 討論

物種基因組全序列的測(cè)定推動(dòng)了生物信息學(xué)的迅速發(fā)展，在海量數(shù)據(jù)的基礎(chǔ)上，利用生物信息學(xué)手段，對(duì)物種基因家族進(jìn)行高效的統(tǒng)計(jì)分類和分析，預(yù)測(cè)基因家族的結(jié)構(gòu)、功能及作用機(jī)制，將極大地推動(dòng)相關(guān)功能基因的挖掘和農(nóng)藝性狀遺傳的改良進(jìn)程[30]。隨著2018年第2版毛竹基因組數(shù)據(jù)的公布以及大量毛竹轉(zhuǎn)錄組數(shù)據(jù)的共享，毛竹GRF基因家族的生物信息學(xué)分析成為可能[11]。本研究通過(guò)全基因組數(shù)據(jù)分析發(fā)現(xiàn)：毛竹GRF家族成員共13個(gè)，數(shù)量多于水稻，可能的原因是毛竹染色體經(jīng)過(guò)加倍，基因組數(shù)據(jù)遠(yuǎn)大于水稻；另外，共線性分析進(jìn)一步證實(shí)：正是通過(guò)基因復(fù)制擴(kuò)增，毛竹GRF在數(shù)量上有優(yōu)勢(shì)。毛竹GRF基因家族各成員間的理化性質(zhì)存在一定的差異，但均含有14/3/3蛋白質(zhì)結(jié)構(gòu)域，其中有6種基序在每個(gè)成員中均出現(xiàn)。根據(jù)基因結(jié)構(gòu)將PeGRF分為ε組和非ε組，其中ε組可能保留了祖先的蛋白質(zhì)功能，這與PIOTROWSKI等[31]和WANG等[32]的研究結(jié)果相似。

大量研究表明GRF蛋白質(zhì)參與激素信號(hào)的轉(zhuǎn)導(dǎo)。如在擬南芥的研究中發(fā)現(xiàn)：GRF參與油菜素類激素(BR)調(diào)控細(xì)胞核發(fā)育的途徑[33]；在煙草Nicotiana tabacum中，GRF參與赤霉素(GA)生物合成調(diào)控[34]；在水稻中，GRF表達(dá)同脫落酸(ABA)密切相關(guān)[35]。本研究發(fā)現(xiàn)：毛竹GRF順式作用元件存在許多激素相關(guān)元件。由此可以推測(cè)毛竹GRF蛋白質(zhì)可能介導(dǎo)激素信號(hào)的轉(zhuǎn)導(dǎo)過(guò)程。但毛竹GRF同其他激素的相互關(guān)系還需進(jìn)一步驗(yàn)證。

GRF蛋白質(zhì)參與了植物的生長(zhǎng)發(fā)育，特別是在花器官的發(fā)育中具有重要作用。PERTL等[36]證實(shí)隨著百合Lilium brownii var. viridulum花粉管的生長(zhǎng)，GRF蛋白質(zhì)的表達(dá)量也明顯增加。李兵娟[37]也證實(shí)雷竹Phyllostachys violascens GRF基因參與開(kāi)花調(diào)控機(jī)制。本研究通過(guò)轉(zhuǎn)錄組數(shù)據(jù)分析發(fā)現(xiàn)：GRF蛋白質(zhì)在花序組織中高表達(dá)，且表達(dá)量明顯高于竹葉和竹鞭，這表明毛竹GRF基因可能參與花序的發(fā)育和調(diào)控。除此之外，在研究毛竹GRF順式作用元件時(shí)還發(fā)現(xiàn)其啟動(dòng)子區(qū)域存在許多光響應(yīng)元件，結(jié)合光周期對(duì)植物開(kāi)花的作用機(jī)制以及在模式植物水稻上的研究[38]，GRF基因可能是通過(guò)光響應(yīng)元件接受外界環(huán)境信號(hào)從而觸發(fā)其高表達(dá)，最終影響毛竹花的發(fā)育。由于受毛竹花發(fā)育相關(guān)材料的限制，該假設(shè)將在后續(xù)實(shí)驗(yàn)驗(yàn)證。

毛竹GRF蛋白質(zhì)是以一個(gè)螺旋結(jié)構(gòu)為主的同源二聚體，二聚體界面內(nèi)包著多個(gè)疏水殘基和多個(gè)極性殘基，外周則由鹽橋連接，三級(jí)結(jié)構(gòu)呈“W”型，每個(gè)單體分別含有2個(gè)凹槽，可能用于結(jié)合配體靶蛋白質(zhì)。毛竹GRF蛋白質(zhì)序列在進(jìn)化譜系中高度保守，并且與配體結(jié)合的氨基酸殘基極端保守，這同SEHNKE等[28]發(fā)現(xiàn)的結(jié)果相似。另外，雖然毛竹GRF蛋白質(zhì)的N端和C端同源性較低，但可能通過(guò)堿性簇維持空間構(gòu)象的穩(wěn)定[28]。PAUL等[39]在研究擬南芥GRF蛋白質(zhì)時(shí)發(fā)現(xiàn)，GRF蛋白質(zhì)還可以通過(guò)結(jié)合磷酸化的蛋白質(zhì)，參與重力反應(yīng)等生理過(guò)程。GRF蛋白質(zhì)在進(jìn)化上高度保守，毛竹PeGRF可能也具有相似的分子作用機(jī)制。但毛竹GRF蛋白質(zhì)生物學(xué)功能與上述空間結(jié)構(gòu)之間的關(guān)系還需進(jìn)一步的探索。