張頔 高曉陽 張軒 劉長寧

摘要 探索了花發(fā)育基因調(diào)控網(wǎng)絡(luò)的3個(gè)拓?fù)渲行男裕ㄟB接度、中間性和接近度）與其中基因分子進(jìn)化速率的相關(guān)性。結(jié)果發(fā)現(xiàn)，隨著網(wǎng)絡(luò)中心性的增加，基因的序列將更加趨向于保守，即基因的進(jìn)化速率與中心性參數(shù)呈負(fù)相關(guān)。這一結(jié)果與擬南芥蛋白質(zhì)相互作用網(wǎng)絡(luò)中所觀察到的模式一致，也許是因?yàn)槎嘈灾萍s了進(jìn)化。

關(guān)鍵詞 花發(fā)育基因;基因調(diào)控網(wǎng)絡(luò);分子進(jìn)化;網(wǎng)絡(luò)拓?fù)渲行男?/p>

中圖分類號 Q75文獻(xiàn)標(biāo)識碼 A文章編號 0517-6611（2021）08-0001-04

doi：10.3969/j.issn.0517-6611.2021.08.001

開放科學(xué)（資源服務(wù)）標(biāo)識碼（OSID）：

Correlation between Molecular Evolution of Flower Development Related Genes and the Topological Centralities of the Flower Development Regulatory Network

ZHANG Di1，2， GAO Xiao-yang1， ZHANG Xuan1，2 et al

（1.Xishuangbanna Tropical Botanical Garden， Chinese Academy of Sciences， Menglun， Yunnan 666303;2.University of Chinese Academy of Sciences， Beijing 100049）

Abstract In this study， we investigated the correlation between three topological centralities （connectivity， betweenness and closeness）of regulatory network of flower development， and the molecular evolutionary rate of the related genes. It was found that with the increase of network centralities， the gene sequences would tend to be more conservative， i.e.， the evolutionary rate of genes was negatively correlated with the network centrality parameters. These results were consistent with the patterns observed in the Arabidopsis thaliana protein-protein interaction network， may be because pleiotropy constrains evolution.

Key words Flower development genes;Gene regulatory network;Molecular evolution;Network topological centrality

基因總是處在一個(gè)特定的生物功能網(wǎng)絡(luò)環(huán)境中發(fā)揮作用。這些功能網(wǎng)絡(luò)是由分子及其相互間的作用構(gòu)成，如蛋白質(zhì)相互作用網(wǎng)絡(luò)是由蛋白質(zhì)及其之間的相互作用構(gòu)成，而基因調(diào)控網(wǎng)絡(luò)是由轉(zhuǎn)錄因子、被調(diào)控的靶基因及它們之間的轉(zhuǎn)錄調(diào)控關(guān)系組成。目前大規(guī)模的基因組、蛋白質(zhì)組和相互作用組等組學(xué)數(shù)據(jù)的涌現(xiàn)，為了解這些真實(shí)生物中運(yùn)行的網(wǎng)絡(luò)提供了新的機(jī)遇，同時(shí)也為分子進(jìn)化研究提供了一個(gè)新的視角。從傳統(tǒng)地關(guān)注單個(gè)基因的分子進(jìn)化，轉(zhuǎn)移到在整個(gè)網(wǎng)絡(luò)中的分子進(jìn)化研究時(shí)，人們通常會關(guān)注基因在網(wǎng)絡(luò)中的拓?fù)鋵傩耘c基因進(jìn)化速率的關(guān)聯(lián)，之前已有許多研究發(fā)現(xiàn)生物網(wǎng)絡(luò)中反應(yīng)基因中心性的拓?fù)鋵傩裕ㄟB接度、中間性和接近度）與進(jìn)化速率有相關(guān)性。如在酵母和果蠅的蛋白質(zhì)-蛋白質(zhì)相互作用網(wǎng)絡(luò)中的連接度[1-2]、中間性以及接近度[3]，人類、擬南芥、水稻、番茄、葡萄和玉米的基因共表達(dá)網(wǎng)絡(luò)中的連接度[4-5]和酵母的轉(zhuǎn)錄調(diào)控網(wǎng)絡(luò)中的中間性相繼被報(bào)道[6]，但對花發(fā)育基因調(diào)控網(wǎng)絡(luò)的研究相對較少。

擬南芥的花發(fā)育過程是由一個(gè)復(fù)雜精細(xì)的基因調(diào)控網(wǎng)絡(luò)控制[7]。首先，在該網(wǎng)絡(luò)的上游，開花時(shí)間受到多條成花途徑的調(diào)控，這些途徑匯聚在信號整合因子上。然后這些整合基因會激活分生組織身份基因，最后再激活花器官特征基因，進(jìn)而調(diào)節(jié)不同花器官身份（如花萼、花瓣、雄蕊、心皮和胚珠）的分化過程。研究發(fā)現(xiàn)花發(fā)育調(diào)控網(wǎng)絡(luò)幾個(gè)階段的基因進(jìn)化速率之間有所差異[8]，但在擬南芥的花器官細(xì)胞身份調(diào)控網(wǎng)絡(luò)中，發(fā)現(xiàn)進(jìn)化速率與所研究的拓?fù)鋵傩灾g無顯著相關(guān)性[9]?；ㄆ鞴偌?xì)胞身份調(diào)控網(wǎng)絡(luò)能否代表整體花發(fā)育網(wǎng)絡(luò)的性質(zhì)，目前尚不清楚。

筆者收集和整理了擬南芥的花發(fā)育相關(guān)核心基因，并計(jì)算了其蛋白質(zhì)編碼序列在4個(gè)十字花科植物（Arabidopsis lyrata、Brassica oleracea、Brassica rapa和Capsella rubella）中的進(jìn)化速率，進(jìn)一步評估了基因的進(jìn)化速率與其所在花發(fā)育調(diào)控網(wǎng)絡(luò)中的拓?fù)渲行男裕ㄟB接度、中間性和接近度）之間的相關(guān)性。

1 材料與方法

1.1 花發(fā)育基因調(diào)控網(wǎng)絡(luò)

花發(fā)育基因調(diào)控網(wǎng)絡(luò)中的數(shù)據(jù)主要來自于Pajoro等[7]，對參與花發(fā)育的基因及其相互作用關(guān)系進(jìn)行了綜述，其中的調(diào)控關(guān)系主要是通過染色質(zhì)免疫共沉淀（ChIP）試驗(yàn)來確定。除2個(gè)miRNA后，基因調(diào)控網(wǎng)絡(luò)包含38個(gè)蛋白質(zhì)編碼基因和201個(gè)轉(zhuǎn)錄調(diào)控關(guān)系（圖1）。借助擬南芥信息資源（TAIR）[10]中的基因功能數(shù)據(jù)信息，調(diào)控網(wǎng)絡(luò)中的基因被分成3個(gè)主要類別。

1.2 直系同源基因的識別

十字花科的5個(gè)已測序物種（Arabidopsis thaliana、Arabidopsis lyrata、Brassica oleracea、Brassica rapa和Capsella rubella）的蛋白質(zhì)和CDS（nucleotide coding sequences）序列的數(shù)據(jù)分別從Phytozome、NCBI和Ensembl Plants基因組數(shù)據(jù)庫上下載。

為了找出每個(gè)擬南芥基因在其他4個(gè)物種中的直系同源基因，在擬南芥的蛋白質(zhì)序列和另一個(gè)物種的蛋白質(zhì)序列之間分別進(jìn)行了雙向最優(yōu)局部比對搜索BLAST[11]（使用1e-15的E-value）。結(jié)果再次進(jìn)行過濾和篩查。直系同源基因過濾和篩查控制標(biāo)準(zhǔn)：比較少的缺失氨基酸，以及相似度高的基因。

1.3 網(wǎng)絡(luò)中心性的計(jì)算

對于每個(gè)基因，使用Python軟件包NetworkX[12]分別計(jì)算了連接度（degree）、中間性（betweenness）和接近度（closeness）3種拓?fù)渚W(wǎng)絡(luò)中心性。其中，連接度是指一個(gè)節(jié)點(diǎn)所連接的節(jié)點(diǎn)數(shù)，也被稱作“connectivity”;中間性是指所有最短路徑通過一個(gè)節(jié)點(diǎn)的比例;接近度是指一個(gè)節(jié)點(diǎn)與所有其他節(jié)點(diǎn)的平均距離的倒數(shù)。

1.4 基因進(jìn)化參數(shù)（dN/dS、dN和dS）的計(jì)算

首先使用MAFFT[13]（--auto參數(shù)）對每一組直系同源基因的蛋白質(zhì)序列進(jìn)行多重序列比對。隨后，根據(jù)蛋白質(zhì)的比對結(jié)果，使用PAL2NAL[14]對CDS序列進(jìn)行了基于密碼子的多重序列比對。使用Gblocks0.91b軟件[15]將比對結(jié)果中質(zhì)量較差的部分進(jìn)行了裁剪，使用的參數(shù)：-t=c -b4=5 -b5=h。此外，應(yīng)用PhyML 2.4軟件[16]中的最大似然方法進(jìn)行了系統(tǒng)發(fā)育分析，其中配置參數(shù)用SMS[17]進(jìn)行優(yōu)化選擇?；虻倪M(jìn)化速率（dN/dS）、非同義替換率（dN）和同義替換率（dS）是基于上述CDS密碼子比對和進(jìn)化樹，使用PAML 4[18]中的CODEML軟件進(jìn)行計(jì)算。

1.5 參數(shù)間相關(guān)性的統(tǒng)計(jì)分析

運(yùn)用R語言3.4.4（http：∥www.r-project.org/）的環(huán)境進(jìn)行有關(guān)統(tǒng)計(jì)分析。網(wǎng)絡(luò)中心性（連接度、中間性和接近度）與進(jìn)化相關(guān)參數(shù)（dN/dS、dN和dS）之間的相關(guān)性，采用斯皮爾曼等級相關(guān)系數(shù)（spearmans rank correlation coefficient）來進(jìn)行衡量（使用了R中的Hmisc包）。如果相應(yīng)零假設(shè)檢驗(yàn)的P< 0.05，則認(rèn)為被檢測的網(wǎng)絡(luò)屬性與進(jìn)化相關(guān)參數(shù)之間有顯著相關(guān)性。在整個(gè)計(jì)算流程中，自行編寫Perl和Python腳本進(jìn)行數(shù)據(jù)格式調(diào)整。

2 結(jié)果與分析

2.1 網(wǎng)絡(luò)中心性的計(jì)算結(jié)果

使用的是Pajoro等[7]的研究中所總結(jié)的花發(fā)育基因調(diào)控網(wǎng)絡(luò)，網(wǎng)絡(luò)中包括38個(gè)編碼基因和201個(gè)互作關(guān)系。分別計(jì)算了網(wǎng)絡(luò)中每個(gè)節(jié)點(diǎn)反應(yīng)網(wǎng)絡(luò)中心性的3個(gè)參數(shù)（連接度、中間性和接近度），計(jì)算結(jié)果見表1。

2.2 基因進(jìn)化參數(shù)的計(jì)算結(jié)果

通過蛋白質(zhì)雙向最優(yōu)比對和質(zhì)量過濾后，得到了花發(fā)育相關(guān)基因在其他4種十字花科植物中的直系同源基因。結(jié)果顯示，在研究的擬南芥38個(gè)基因中，有35個(gè)基因在4種植物中都找到了雙向最優(yōu)比對的基因，但其中FLM、PI和SNZ基因在某些植物中沒有找到同源基因。通過對沒有缺失值的35組直系同源基因使用PAML中的最大似然算法估算出了dN/dS、dN和dS。從dN對dS的比值（dN/dS）推斷自然選擇的影響。通常，適應(yīng)性變化可以通過分子水平計(jì)算非同義替代速率dN與同義替代速率dS的比值進(jìn)行分析。如果沒有選擇作用，或沒有很強(qiáng)的有害突變，同義與非同義替代的速率相同，則dN/dS=1;如果存在負(fù)選擇，則dN/dS<1;如果存在正選擇，則dN/dS>1。因此，dN/dS不僅可以用來檢測選擇作用，還可以用來確定選擇方向。所研究基因的dN/dS值都小于0.4（表2），平均值為0.178 8。表明這些基因總體上都是在純化選擇下進(jìn)化的，只是在進(jìn)化過程中受到了不同的選擇約束。

2.3 進(jìn)化參數(shù)與網(wǎng)絡(luò)中心性之間的相關(guān)性分析

分別計(jì)算了基因編碼區(qū)進(jìn)化參數(shù)與網(wǎng)絡(luò)中心性之間的相關(guān)性，結(jié)果見表3。經(jīng)過統(tǒng)計(jì)檢驗(yàn)，發(fā)現(xiàn)基因序列的非同義替換率（dN）與同義替換率的（dS）的比值（dN/dS）與網(wǎng)絡(luò)中心性（連接度、中間性和接近度）呈顯著負(fù)相關(guān)，這可能指示著處于網(wǎng)絡(luò)中央的基因受到了更多的功能限制，而傾向于減少非同義替換的純化選擇。dN的相關(guān)性也反映了類似的負(fù)相關(guān)趨勢。其中與中間性的負(fù)相關(guān)性最顯著，與連接度的負(fù)相關(guān)性次之，而與接近度的負(fù)相關(guān)性不顯著。另外在dS與網(wǎng)絡(luò)中心性之間沒有發(fā)現(xiàn)相關(guān)性。

該研究結(jié)果與之前在蛋白質(zhì)-蛋白質(zhì)相互作用網(wǎng)絡(luò)[3，19]和共表達(dá)網(wǎng)絡(luò)[4-5]中所得出的研究結(jié)果一致。似乎在這些分子網(wǎng)絡(luò)中，連接度越高的基因可能會因?yàn)榫哂卸嘈?，而在進(jìn)化上更加保守[20];這表明由于越中心的基因序列改變會對生物體產(chǎn)生更有害的影響。那么為何在擬南芥的花器官細(xì)胞身份調(diào)控網(wǎng)絡(luò)中沒有發(fā)現(xiàn)進(jìn)化速率與拓?fù)鋵傩灾g有顯著的相關(guān)性[9]，原因之一可能在于，根據(jù)Liu等[8]所述花發(fā)育基因調(diào)控網(wǎng)絡(luò)所分成的幾個(gè)功能子網(wǎng)，該花器官細(xì)胞身份調(diào)控網(wǎng)絡(luò)大約代表整個(gè)網(wǎng)絡(luò)中的一個(gè)子網(wǎng);所以雖然對于全局花發(fā)育基因調(diào)控網(wǎng)絡(luò)而言，進(jìn)化速率受到了基因中心性的影響，但同時(shí)又受到了局部功能約束的影響。如Szedlak等[21]研究發(fā)現(xiàn)，人類基因調(diào)控網(wǎng)絡(luò)中，基因的進(jìn)化特性與節(jié)點(diǎn)的中心性度量相關(guān)，同時(shí)在基因聚類簇內(nèi)部的進(jìn)化速率卻相對均一[21]。

3 結(jié)論與討論

該研究探討了在花發(fā)育基因調(diào)控網(wǎng)絡(luò)中每個(gè)基因的網(wǎng)絡(luò)拓?fù)渲行男园ㄟB接度、中間性和接近度對其編碼區(qū)序列進(jìn)化速率的影響。結(jié)果發(fā)現(xiàn)，花發(fā)育基因調(diào)控網(wǎng)絡(luò)總體上在純化選擇下進(jìn)化，但隨著網(wǎng)絡(luò)中心性的增加，基因的序列將更加趨向于保守。這種趨勢對于連接度、中間性和接近度而言方向都是相同的，只是在程度和顯著性上略有差別。該研究為在網(wǎng)絡(luò)背景下理解花發(fā)育相關(guān)基因的分子進(jìn)化提供了新的數(shù)據(jù)。由于目前關(guān)于擬南芥花發(fā)育基因調(diào)控網(wǎng)絡(luò)數(shù)據(jù)的數(shù)量和質(zhì)量仍在不斷發(fā)展中[1]，因此該研究所得出的相關(guān)性結(jié)論也受限于所選用的調(diào)控網(wǎng)絡(luò)數(shù)據(jù)來源。另外，在實(shí)際的生物系統(tǒng)中，除網(wǎng)絡(luò)拓?fù)渲行男酝?，有可能還有其他生物學(xué)參數(shù)也影響了基因序列的進(jìn)化，如表達(dá)水平和功能類別。因此，在今后的相關(guān)性研究中，隨著可用數(shù)據(jù)信息增加，可以研究更多的基因和生物學(xué)參數(shù)，以便進(jìn)一步地了解花發(fā)育調(diào)控網(wǎng)絡(luò)的保守性和可進(jìn)化性。

參考文獻(xiàn)

[1]LEMOS B，BETTENCOURT B R，MEIKLEJOHN C D，et al.Evolution of proteins and gene expression levels are coupled in Drosophila and are independently associated with mRNA abundance，protein length，and number of protein-protein interactions [J].Molecular biology and evolution，2005，22（5）：1345-1354.

[2]FRASER H B，HIRSH A E，STEINMETZ L M，et al.Evolutionary rate in the protein interaction network [J].Science，2002，296（5568）：750-752.

[3]HAHN M W，KERN A D.Comparative genomics of centrality and essentiality in three eukaryotic protein-interaction networks [J].Molecular biology and evolution，2005，22（4）：803-806.

[4]MASALIA R R，BEWICK A J，BURKE J M.Connectivity in gene coexpression networks negatively correlates with rates of molecular evolution in flowering plants [J].PLoS One，2017，12（7）：1-10.

[5]JORDAN I K，MARIO-RAMREZ L，WOLF Y I，et al.Conservation and coevolution in the scale-free human gene coexpression network [J].Molecular biology and evolution，2004，21（11）：2058-2070.

[6]JOVELIN R，PHILLIPS P C.Evolutionary rates and centrality in the yeast gene regulatory network [J].Genome biology，2009，10（4）：1-10.

[7]PAJORO A，BIEWERS S，DOUGALI E，et al.The（r）evolution of gene regulatory networks controlling Arabidopsis plant reproduction：A two-decade history [J].Journal of experimental botany，2014，65（17）：4731-4745.

[8]LIU Y，GUO C C，XU G X，et al.Evolutionary pattern of the regulatory network for flower development：Insights gained from a comparison of two Arabidopsis species [J].Journal of systematics and evolution，2011，49（6）：528-538.

[9]DAVILA-VELDERRAIN J，SERVIN-MARQUEZ A，ALVAREZ-BUYLLA E R.Molecular evolution constraints in the floral organ specification gene regulatory network module across 18 angiosperm genomes[J].Molecular biology and evolution，2014，31（3）：560-573.

[10]BERARDINI T Z，REISER L，LI D，et al.The Arabidopsis information resource：Making and mining the “gold standard” annotated reference plant genome [J].Genesis，2015，53（8）：474-485.

[11]CAMACHO C，COULOURIS G，AVAGYAN V，et al.BLAST+：Architecture and applications [J].BMC Bioinformatics，2009，10（1）：1-9.

[12]HAGBERG A A，SCHULT D A，SWART P J.Exploring network structure，dynamics，and function using networkx [C]//VAROQUAUX G，VAUGHT T，MILLMAN J，et al.Proceedings of the 7th Python in Science Conference.Pasadena，CA USA：[s.n.]，2008.

[13]KATOH K，STANDLEY D M.MAFFT multiple sequence alignment software version 7：Improvements in performance and usability [J].Molecular biology and evolution，2013，30（4）：772-780.

[14]SUYAMA M，TORRENTS D，BORK P.PAL2NAL：Robust conversion of protein sequence alignments into the corresponding codon alignments [J].Nucleic acids research，2006，34：W609-W612.

[15]CASTRESANA J.Selection of conserved blocks from multiple alignments for their use in phylogenetic analysis [J].Molecular biology and evolution，2000，17（4）：540-552.

[16]GUINDON S，DUFAYARD J F，LEFORT V，et al.New algorithms and methods to estimate maximum-likelihood phylogenies：Assessing the performance of PhyML 3.0 [J].Systematic biology，2010，59（3）：307-321.

[17]LEFORT V，LONGUEVILLE J E，GASCUEL O.SMS：Smart model selection in PhyML [J].Molecular biology and evolution，2017，34（9）：2422-2424.

[18]YANG Z.PAML 4：Phylogenetic analysis by maximum likelihood [J].Molecular biology and evolution，2007，24（8）：1586-1591.

[19]ALVAREZ-PONCE D，F(xiàn)ARES M A.Evolutionary rate and duplicability in the Arabidopsis thaliana protein-protein interaction network [J].Genome biology and evolution，2012，4（12）：1263-1274.

[20]HE X，ZHANG J.Toward a molecular understanding of pleiotropy [J].Genetics，2006，173（4）：1885-1891.

[21]SZEDLAK A，SMITH N，LIU L，et al.Evolutionary and topological properties of genes and community structures in human gene regulatory networks [J].PLoS Computational Biology，2016，12（6）：1-16.