花發(fā)育相關基因分子進化與花發(fā)育調控網絡拓撲中心性的相關性研究

張頔 高曉陽 張軒 劉長寧

摘要 探索了花發(fā)育基因調控網絡的3個拓撲中心性（連接度、中間性和接近度）與其中基因分子進化速率的相關性。結果發(fā)現(xiàn)，隨著網絡中心性的增加，基因的序列將更加趨向于保守，即基因的進化速率與中心性參數呈負相關。這一結果與擬南芥蛋白質相互作用網絡中所觀察到的模式一致，也許是因為多效性制約了進化。

關鍵詞 花發(fā)育基因;基因調控網絡;分子進化;網絡拓撲中心性

中圖分類號 Q75文獻標識碼 A文章編號 0517-6611（2021）08-0001-04

doi：10.3969/j.issn.0517-6611.2021.08.001

開放科學（資源服務）標識碼（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

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

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

筆者收集和整理了擬南芥的花發(fā)育相關核心基因，并計算了其蛋白質編碼序列在4個十字花科植物（Arabidopsis lyrata、Brassica oleracea、Brassica rapa和Capsella rubella）中的進化速率，進一步評估了基因的進化速率與其所在花發(fā)育調控網絡中的拓撲中心性（連接度、中間性和接近度）之間的相關性。

1 材料與方法

1.1 花發(fā)育基因調控網絡

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

1.2 直系同源基因的識別

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

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

1.3 網絡中心性的計算

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

1.4 基因進化參數（dN/dS、dN和dS）的計算

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

1.5 參數間相關性的統(tǒng)計分析

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

2 結果與分析

2.1 網絡中心性的計算結果

使用的是Pajoro等[7]的研究中所總結的花發(fā)育基因調控網絡，網絡中包括38個編碼基因和201個互作關系。分別計算了網絡中每個節(jié)點反應網絡中心性的3個參數（連接度、中間性和接近度），計算結果見表1。

2.2 基因進化參數的計算結果

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

2.3 進化參數與網絡中心性之間的相關性分析

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

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

3 結論與討論

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

參考文獻

[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.