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

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

（1. 浙江農林大學 省部共建亞熱帶森林培育國家重點實驗室，浙江 杭州 311300；2. 浙江農林大學 林業與生物技術學院，浙江 杭州 311300）

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

1 材料與方法

1.1 基因家族來源、鑒定及理化性質分析

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

1.2 家族進化分析

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

1.3 基因結構、基序和保守結構域預測

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

1.4 啟動子分析

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

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

利用 MCScanX[25]獲取GRF家族種內、種間共線性關系，并用 TBtools軟件 Amazing Super Circos[26]和Multipe Synteny Plot分別對種內和種間的結果可視化。

1.6 基因表達分析

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

1.7 蛋白質三級結構同源模建

利用SWISSMODEL(https://www.swissmodel.expasy.org/)在線軟件[27]預測GRF蛋白質的3D結構。模建結果使用SAVES v5.0(https://servicesn.mbi.ucla.edu/SAVES/)[19]進行評估。

2 結果與分析

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

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

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

2.2 GRF基因家族分類與進化樹構建

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

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

2.3 GRF家族基因結構、基序及保守結構域

對毛竹GRF基因結構分析發現：內含子數量存在差異，非ε組成員都包含4個外顯子和3個內含子，它們在位置上高度保守。ε組成員都具有不同于非ε組的內含子-外顯子結構，具有2個額外的N-末端內含子[21]。利用NCBI-CDD對毛竹GRF基因進行保守結構域分析，PeGRF蛋白質均包含14/3/3結構域，毛竹GRF基因家族14/3/3結構域存在一定的保守性，但該結構域的分布位置有一定分化。利用MEME在線工具對該基因家族的保守基序預測，基數設置為10，結果顯示(圖2)：Motif1～6在每個家族成員中均出現，屬于高度保守結構，其余基序在家族成員中出現的頻率及所在位置均存在一定的差異。

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

2.4 啟動子特征

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

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

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

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

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

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

2.6 GRF家族基因表達模式

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

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

2.7 GRF家族蛋白質空間三級結構

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

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

3 討論

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

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

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

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