A contribution to the study of plant development evolution based on gene co-expression networks

ORIGINAL RESEARCH ARTICLE published: 05 August 2013 doi: 10.3389/fpls.2013.00291 A contribution to the study of plant development evolution based on gene co-expression networks Francisco J. Romero-Campero1, Eva Lucas-Reina2, Fatima E. Said2,JoséM.Romero 2and Federico Valverde2* 1Department of Computer Science and Artificial Intelligence, Universidad de Sevilla, Sevilla, Spain 2Molecular Plant Development and Metabolism, Instituto de Bioquímica Vegetal y Fotosíntesis, Consejo Superior de Investigaciones Científicas y Universidad de Sevilla, Sevilla, Spain Edited by: Madelaine E. Bartlett, Brigham Young University, USA Reviewed by: Stephan Wenkel, University of Tuebingen, Germany Enamul Huq, The University of Texas at Austin, USA Idan Efroni, New York University, USA *Correspondence: Federico Valverde, Molecular Plant Development and Metabolism Group,Instituto de Bioquímica Vegetal y Fotosíntesis, Consejo Superior de Investigaciones Científicasy Universidad de Sevilla, 49th, Americo Vespucio Avenue, 41092 Sevilla, Spain e-mail: federico.valverde@ ibvf.csic.es Phototrophic eukaryotes are among the most successful organisms on Earth due to their unparalleled efficiency at capturing light energy and fixing carbon dioxide to produce organic molecules. A conserved and efficient network of light-dependent regulatory modules could be at the bases of this success. This regulatory system conferred early advantages to phototrophic eukaryotes that allowed for specialization, complex developmental processes and modern plant characteristics. We have studied light-dependent gene regulatory modules from algae to plants employing integrative-omics approaches based on gene co-expression networks. Our study reveals some remarkably conserved ways in which eukaryotic phototrophs deal with day length and light signaling. Here we describe how a family of Arabidopsis transcription factors involved in photoperiod response has evolved from a single algal gene according to the innovation, amplification and divergence theory of gene evolution by duplication. These modifications of the gene co-expression networks from the ancient unicellular green algae Chlamydomonas reinhardtii to the modern brassica Arabidopsis thaliana may hint on the evolution and specialization of plants and other organisms. Keywords: photoperiod, evolution, gene co-expression networks, Chlamydomonas,Physcomitrella,Arabidopsis INTRODUCTION Day length, or photoperiod, regulates strategic developmental processes in plants. It constitutes a key external signal to feed the circadian rhythms with the needed information to maintain the set of the day (Imaizumi, 2010) and is also critical to discern seasons (Jackson, 2009). Thus, photoperiodic signals regulate crucial developmental responses such as dormancy; germination; senescence or the transition from vegetative to reproductive stages (Valverde, 2011). The floral transition is one of the most conserved evolutionary processes in angiosperms (Romero and Valverde, 2009; Serrano et al., 2009) due to its central role in producing new plant generations and in the transmission of acquired characteristics. Many of the key genes involved in the floral transition are conserved in an inter-species dependent manner as well as many of the external cues that trigger the reproductive response (Amasino, 2010). According to this, the response to the three most important external agents that control reproduction, namely temperature, nutrients and day length, are regulated by a set of inter-species gene regulatory networks with conserved functions (Ausín et al., 2005). The photoperiodic flowering pathway involves a series of genes that are influenced by light and circadian signals to orchestrate a response that ensures the best moment of the year to flower. Thus, long-day (LD) plants flower as the day lengthens and shortday (SD) plants flower when days start to shorten. In the model species Arabidopsis, a facultative LD plant, a detailed knowledge of the gene pathways involved in the photoperiod floral transition has been accumulated and solid mechanisms have been proposed to explain its flowering behavior (Fornara et al., 2010). In this pathway, the role of the gene CONSTANS (CO) is crucial because its fine control at several regulatory levels (Valverde et al., 2004) assure that the plant triggers the flowering response exactly at the precise moment of the year. It modulates the expression of the florigen FLOWERING LOCUS T (FT) gene in the leaf vasculature, whose protein, reaching the meristem, triggers the reproductive developmental program, eventually producing flowers (Corbesier et al., 2007). While many genes have been discovered in this regulatory module, dealing with protein stability (Jang et al., 2008; Lázaro et al., 2012), modulating the light response of the proteins (Yu et al., 2008)ortheirfunction(Kim et al., 2007), the integration of the pathway within other flowering routes or the output genes that trigger the flowering transition are less known. Due to the characteristic attached-to-the-soil behavior of plants, a complex network of regulatory processes is in the base of their physiological responses. These complex systems can be better analysed employing holistic and integrative approaches (Usadel et al., 2009; Tohge and Fernie, 2012). The last years have seen a bloom of massive data acquisition techniques to approach plant biology in a holistic and integrative way, particularly in transcriptional information (Metzker, 2010). Microarrays first and Next-Generation Sequencing (NGS) later have provided enormous amount of gene expression data for a multitude of plant species in different physiological/genotypic conditions and developmental stages (Schliesky www.frontiersin.org August 2013 | Volume 4 | Article 291 |1 Romero-Campero et al. Conserved light signaling in plants et al., 2012). These data can also be enriched by further experimental approaches and include this information into the gene co-expression networks generated. Therefore, when trying to analyze the developmental response of a plant to an external condition we can combine data from gene co-expression analyses generated using microarray and/or NGS approaches with physiological data such as the time to flower, the weight of the plant or the chlorophyll content. This helps to associate a particular plant behavior to a particular gene expression pattern. Thus, the levels of complexity poised by the enormous amounts of data generated at different regulatory levels, such as transcriptomics, proteomics and metabolomics can be approached with an integrative and holistic perspective employing network tools to reduce noise and find novel patterns of organization. When an evolutionary approach covering the history of gene networks within the phylogeny of organisms is employed, an interesting relation between function and diverse processes, such as developmental processes, starts to be unveiled. This could help to explain many of the intriguing intertwines between evolution and development observed in different organisms (Müller, 2007). Sequence similarity constitutes the classical approach to assign potential functions to genes and to study their evolutionary history (Lajoie et al., 2010). This methodology focuses on the comparison between individual genes and do not take into account that genes perform their function in coordination with many other genes. In this respect, sequence similarity has been shown to be incomplete when predicting phenotypic differences between species. For example, the genes involved in the human and chimpanzee brains share very high sequence similarities that do not correspond with the marked phenotypic differences between them (Oldham et al., 2006). Holistic and integrative approaches such as gene co-expression networks that take into account the orchestration among genes are emerging as powerful tools to predict gene function and to infer their evolutionary history. This methodology assumes that if a group of genes are co-expressed they should have similar functions and a common evolutionary history. In these networks, genes are represented as nodes and an edge is established between two nodes if the expression profiles of the corresponding genes exhibit a correlation value high enough to provide evidence of co-expression. The topological analysis of these networks such as the distribution of the number of neighbors of each gene, the number of co-expressed genes for a given gene, can provide information about their function and evolutionary history (Aoki et al., 2007; Usadel et al., 2009). In this paper, we have studied the evolution of the coexpression sub-networks or modules around the photoperiod central family of CO-Like genes, CONSTANS (CO) homologs, in three model species whose genome is available: the green unicellular alga Chlamydomonas reinhardtii (CrCO gene), the moss Physcomitrella patens (PpCOL genes) and the higher plant Arabidopsis thaliana (AtCOL genes). These species are landmarks of the evolutionary lineage of plants. Microarrays or RNA-seq database experiments have been used to construct gene coexpression networks. Phylogenetic analyses have been combined together with functional enrichment analyses in order to better understand the evolution of the CO-Like (COLs) family between the species and their functional specialization. This combination of gene expression data analysis, gene ontology (GO) term enrichment and phylogenetic studies constitutes a novel methodology to study gene function and evolution. The approach is not restricted to our case study and can be applied to the study of sub-networks or modules of other transcription factors. Through the analysis of the evolution of gene networks we can study gene duplication and diversification and how this has affected the networks. Additionally, our methodology could be a useful tool to identify homologous genes related to the same specific process and therefore predict true gene orthology. Finally, our analysis can also be used to explain why and to what extent different signaling pathways are linked in an organism and, therefore, constitute a valuable tool to understand plant plasticity. MATERIALS AND METHODS PLANT, ALGAL MATERIAL, AND GROWTH CONDITIONS Arabidopsis thaliana 35S:CO-GR (Simon et al., 1996) transgenic lines were grown in MS plates. Seeds were previously incubated 4daysat4 ◦C in the dark before sowing under 16 h light/8 h dark cycle (long-day, LD) with temperature ranging from 22◦C(day) to 18◦C(night)at75µE/m2light intensity. In the experiment employing dexamethasone (DEX) 1 µg/mL and cycloheximide (CHX) 1mM chemicals were added 10 days after sowing at ZT0 and leaf samples harvested 4 h after the drug treatment, considering Zeitgeber Time 0 (ZT0) the moment at which the lights are switched on. Chlamydomonas reinhardtii cell-wall deficient mutant CW15 (Davies and Plaskitt, 1971)andpNIA:CrCO (Serrano et al., 2009) transgenic line were grown in stirred conical cylindrical flasks containing Sueoka NO− 3(Sueoka et al., 1967)underLDconditions in control rooms at 22◦Cand50µE/m2. Algal cells were harvested at 4 days at ZT4, 4 h after the lights went on. Plants and algae were grown in a model SG-1400 phytotron (Radiber SA, Spain). RNA ISOLATION AND Q-PCR RNA was isolated from Arabidopsis seedlings (0.1 g leaf tissue) and Chlamydomonas (20 ml of an exponential phase culture) employing, in both cases, the TRIZOL (Invitrogen) protocol as described by the manufacturer. In short, the sample was mixed with 1ml of TRIZOL and 0.2 ml of chloroform, the mixture was then centrifuged at 16,000 g for 10 min at 4◦C. The supernatant was treated with 1 volume of 2-propanol, incubated 15 min at room temperature and centrifuged at 16,000 g for 10 min at 4◦C. 0.75 ml 3 M LiCl was added to the pellet and incubated for t>10 min at room temperature and centrifuged at 16,000 g for 10 min at 4◦C. The pellet was washed with 80% (v/v) ethanol and centrifuged at 16,000 g for 10 min at 4◦C. The final RNA sample was suspended in 30 µl of DEPC treated water and quantified employing a ND-1000 Spectrophotometer (Nanodrop). 1µg of TRIZOL isolated RNA was used to synthesize cDNA employing the Quantitec® Reverse kit (Qiagen) following the instruction recommended by the manufacturer. cDNA was diluted to a final concentration of 10 ng/µLandstoredat−20◦C until Q-PCR was performed. Primers to amplify the 3translated region of AtCO, AtSSS, AtFAD, AtGS, AtZEP, AtWRKY33, AtUBQ10, CrCO, CrSSS, CrFAD, CrGS, CrZEP,CrATG8,and Frontiers in Plant Science | Plant Evolution and Development August 2013 | Volume 4 | Article 291 |2 Romero-Campero et al. Conserved light signaling in plants CrTUB (Table 1) were designed employing an Oligo analyzer program (Integrated DNA technologies, http://eu.idtdna.com/ analyzer/Applications/OligoAnalyzer/). Q-PCR was performed in a Multicolor Real-Time PCR Detection System iQTM5 from Bio-Rad in a 10 µL reaction: primer concentration 0.2 µM, 10 ng cDNA and 5 µL SensiFAST TM SYBR & Fluorescein Kit (Bioline). Each sample was measured by triplicate. The Q-PCR program consisted in (1) 1 cycle (95◦C, 2 min); (2) 40 cycle (95◦C, 5s; 60 ◦CforArabidopsis primers and 65◦CtoChlamydomonas primers, 10 s and 72◦C,6s)(iii)1cycle(72 ◦C, 6 s). Fluorescence was measured at the end of each extension step and the melting curve was performed between 55 and 95◦C. The initial concentration of candidate and reference gene was calculated by means of LingRegPCR software version 11.0 (Ruijter et al., 2009). Normalized data was calculated dividing the average of four replicates of each sample of the candidate and reference genes. PHYLOGENETIC ANALYSIS Evolutionary relationships among the CO-Like genes of Chlamydomonas,Physcomitrella and Arabidopsis were analyzed using predicted amino-acid sequences from Phytozome v9.0 and aligned with the program MUSCLE (Edgar, 2004). This alignment was then used to generate a phylogenetic tree applying the Neighbor-Joining algorithm (Saitou and Nei, 1987)with the JTT+G 1.53 as substitution model (Jones et al., 1992). The bootstrap consensus tree inferred from 1000 replicates is taken to represent the evolutionary history of the analyzed proteins. The tree is drawn to scale, with branch lengths in the same units as those of the evolutionary distances used to infer the phylogenetic tree. Phylogenetic analyzes were conducted with MEGA5 (Tamura et al., 2011). Accession numbers of sequences usedinthealignmentareshowninTable 2. IDENTIFICATION OF CONSERVED MOTIFS The sequences of COL proteins from the different plants and algae were analyzed using the MEME program (http:// meme.nbcr.net/meme/) as described by Bailey et al. (2009). Only the three main motifs (2 B-boxes and CCT domain) were represented. The motifs were characterized using the Table 2 | COL protein family in Clamydomonas, Physcomitrella,and Arabidopsis employed in the phylogenetic studies of this work. Protein Gene number Protein Gene number CrCO g6302 CO At5g15840 PpCOL1 Pp1s371_27V6 AtCOL1 At5g15850 PpCOL2 Pp1s97_109V6 AtCOL2 At3g02380 PpCOL3 Pp1s364_5V6 AtCOL3 At2g24790 PpCOL4 Pp1s36_238V6 AtCOL4 At5g24930 PpCOL5 Pp1s26_5V6 AtCOL5 At5g57660 PpCOL6 Pp1s236_21V6 AtCOL6 At1g68520 PpCOL7 Pp1s195_82V6 AtCOL7 At1g73870 PpCOL8 Pp1s143_52V6 AtCOL8 At1g49130 PpCOL9 Pp1s108_97V6 AtCOL9 At3g07650 PpCOL10 Pp1s3_491V6 AtCOL10 At5g48250 AtCOL11 At4g15250 AtCOL12 At3g21880 AtCOL13 At2g47890 AtCOL14 At2g33500 AtCOL15 At1g28050 AtCOL16 At1g25440 The left column includes the short protein code number as found in the literature (PpCOL4-PpCOL10 named in this work) and the right column the gene code number from the corresponding genome databases (Phytozome and TAIR). Table 1 | Primers employed for Q-PCR experiments. Arabidopsis genes Sequence Amplified fragment size (bp) CO 5-CCAATGGACAGAGAAGCCAGG-35-GCATCGTGTTGAACCCTTGC- 3175 AtSSS 5-CTGGGGATCATCAGCTACACAATACG-35-CACGTGCGATTAGGAACAGCTC-381 AtFAD 5- CGTCGTTAAGTTCCTTCAAGCC -35- CATAGCTTCAATCGAACCGACAG -3157 AtGS 5- CCAGCTTCGAACATGGATCC -35-CCTAAGACATTGCTTGATAGAGAACAC-3167 AtZEP 5- CTCCGAAATCGACGAGGAAG -35- TGCAAGGAATAGCTGAAAGCAG -3166 AtUBQ10 5- GAAGTTCAATGTTTCGTTTCATGT -35- GGATTATACAAGGCCCCAAAA -3119 AtERF 5- CCAATGTTCAGCAGAATGCC -35- GGACGATGAGAAAGAATTAGGAG -385 AtWRKY33 5-TACCGGGCCTTTTGGTTA-35- CCACCACCAACAAAGTTTTG-381 Chlamydomonas genes Sequence bp CrCO 5- CTTCCCGCAAGGCGTATGC -35- GCCTCAATCTCCTCCTTCTTGGC -373 CrSSS 5- ACGTGTACCGCTCCATCAGC -35- GCAGCACTCTTGCACTATGCAG -3107 CrFAD 5- GACGAGAAGGTCAACTACAAGCC -35- GCTTGCTCAGCTCCGATTAGC -3150 CrGS 5- GCTACGGCTACCTGGAGGA -35- CATCGCTGCCCTTATTAGCTGG -3127 CrZEP 5- AAGAGGCAGGTTGGCTTAGTGC -35- GGTGTCTGTCAACGTGTGTAGC -3127 CrTUB 5- GTTGCATCGTTAGCGTGGACG -35- GCAGCAGCCAATGTTCAGACT -3170 CrERF 5- AGCCAGGCTCGCTGCAACTTCC -35- GGA AGT TGC AGC GAG CCT GGC T -3108 CrATG8 5- TCCCGATATCGACAAGAAG -35- TGCGGATGACGTACACAAAT -375 www.frontiersin.org August 2013 | Volume 4 | Article 291 |3 Romero-Campero et al. Conserved light signaling in plants Conserved Domain Search Service (http://www.ncbi.nlm.nih. gov/Structure/cdd/wrpsb.cgi) (Marchler-Bauer and Bryant, 2004). DATA COMPILATION AND PROCESSING This study is based on an ensemble of more than one TeraByte transcriptomic data from Chlamydomonas,Physcomitrella, and Arabidopsis obtained under related physiological and genotype conditions (Table 3). These data comprise RNA-seq experiments from the public database, the Sequence Read Archive (SRA, http://www.ncbi.nlm.nih.gov/Traces/sra/) (Wheeler et al., 2005) and microarray data from the Array Express Archive (ArrayExpress, http://www.ebi.ac.uk/arrayexpress/) (Rustici et al., 2013). For Chlamydomonas 50 RNA-seq data sets representing eight different genotypes under diverse physiological conditions (González-Ballester et al., 2010; Miller et al., 2010; Castruita et al., 2011; Kropat et al., 2011; Boyle et al., 2012; Fischer et al., 2012; Urzica et al., 2012)wereanalyzed.For Physcomitrella 13 RNA-seq data sets representing two different genotypes (Zemach et al., 2010; Xiao et al., 2011; Chen et al., 2012) were analyzed. Finally, in order to use comparable conditions Arabidopsis microarray data from experiments under similar conditions to those from Chlamydomonas and Physcomitrella (Gutiérrez et al., 2007; Long et al., 2010; Patterson et al., 2010; Iyer-Pascuzzi et al., 2011; Cheng et al., 2013)were analyzed. The experimental conditions analyzed in this study are diverse enough to capture the true co-expression among genes. These conditions include nutrient deficiency (nitrogen, iron, copper, and sulfur deprivation), oxidative stress, light stimuli, DNA damage and different developmental stages (Table 3). An estimation of gene expression in the transcriptomes generated in the different conditions under study was obtained following the methodology described in Trapnell et al. (2012). This approach uses a reference genome for each species. The Chlamydomonas reference genome corresponds to the release v5.3 which is assembled into 17 chromosomes and 37 additional unmapped scaffolds. The Chlamydomonas genomic information is based on the Augustus update u11.6. For Physcomitrella we used the reference genome release v1.6 which is assembled into 2106 scaffolds that have not been mapped yet to the 27 chromosomes that constitute its genome. The Physcomitrella genomic information is based on the Cosmoss update v1.6. Finally, the Arabidopsis genomic information is based on the TAIR v10 resource. This annotation and genome information was downloaded from Phytozome (http://www.phytozome.net), a webbased platform for green plant comparative genomics (Goodstein et al., 2012). The workflow followed in the RNA-seq data processing is sketched in Figure 1 and detailed in Figure S1. The first step consisted on the alignment of the short read sequences stored in fastq files to the corresponding reference genome using the software package Tophat (Trapnell et al., 2009)toproduceBAM files (binary alignment maps). In the second step these files together with the gene annotation information available were used to assemble the whole transcriptomes of Chlamydomonas and Physcomitrella under the conditions studied. The assembled Table 3 | Description of the experiments and datasets employed in this work. Experiment Accession Data base number 2137 WT cells were cultivated photoheterotrophically under Fe-replete (20 mM), Fe-deficient (1 mM) and Fe-limited (0.25 mM) conditions SRP010563 SRA 2137 WT cells were exposed to hydrogen peroxide for 0, 0.5, and 1 h SRP010084 SRA 4A+WT cells and sor1 mutant cells were grown in 12 h light/dark cycle for several days. Samples were taken at ZT6 SRP009273 SRA 2137 WT cells and crr1 mutant cells were cultivated in TAP or minimal medium under Cu-sufficient and Cu-defficient conditions SRP005483 SRA CW15 cells were cultivated in N-repleted, N-deprived conditions SRP003630 SRA D66WTcellsandars11 mutant cells were cultivated under continuous light in S-repleted and S-deprived conditions SRP002284 SRA Dark grown WT and double pubs/hy2 mutant protonema were irradiated with red light for 1 h SRP011279 SRA WT protonemal, caulonemal, and chloronemal tissues were collected at 3, 14, 24, 30 days of development SRP009201 SRA Samples were collected in control tissues and tissues treated with the DNA-DSB inducing agent bleomycin SRP004443 SRA Col-0 WT and pye mutant seedlings were grown in Fe-repleted and Fe-deprived media GSE21582 GEO Col-0 plants were grown hydroponically, transferred to a nitrogen free medium for 26 h and finally supplied with 1 mM nitrate or 1 mM ammonium GSE29589 GEO Col-0 plants were hydroponically in nutrient solutions with various concentrations of nitrate and sucrose E-MEXP-828 Array express Col-0 seedlings were cultivated in the dark to reduce endogenous H2O2then were treated with 5mM H2O2 GSE40574 GEO Five days old Col-0 seedlings were transferred to sulfur deficient media. Samples were collected at 0, 3,12, 24, 18, and 72 h GSE30098 GEO Samples from Col-0, CO overexpressor 35S:CO and co-2 mutant were collected in long-day conditions (16 h light/8 h dark) at ZT4 E-MTAB-1078 Array express The columns represent the species name, accession numbers and data bases acronyms for the different experimental datasets. A short descriptive paragraph of the experimental design (as stated by the authors in the databases) is also provided. Frontiers in Plant Science | Plant Evolution and Development August 2013 | Volume 4 | Article 291 |4 Romero-Campero et al. Conserved light signaling in plants FIGURE 1 | Schematic workflow including data collection and processing as well as gene co-expression network construction, analysis and validation. The figure simplifies the main steps followed in the workflow for the construction and analysis of the co-expression networks for Chlamydomonas,Physcomitrella,andArabidopsis.This workflow requires a reference genome and the corresponding Gene Ontology annotation for each species. The annotation was obtained from the Phytozome database. Gene expression estimation measured in FPKM was obtained by aligning short reads to the reference genome and assembling transcripts using Tophat and Cufflinks. Gene correlation and network construction were performed with different Bioconductor R packages. Cytoscape was used for network visualization. Functional patterns in the neighborhood of the genes of interest were identified with the web tools agriGO and GOrilla. Finally, Q-PCR was used for experimental validation. A detailed representation of this workflow can be seen in Figure S1. transcriptomes were stored in GTF files and resolved using the software package Cufflinks and its program Cuffmerge (Trapnell, 2010). Finally, gene expression levels from the different conditions integrated in our study were estimated using Cuffdiff, a program included in the Cufflinks package that takes as input the alignments files (BAM files) and the assembled whole transcriptomes (GTF files). In order to avoid biases due to the length of the different transcripts and number of reads generated in each experiment, Cufflinks estimates gene expression using as unit of measurement the Fragments Per Kb of exon per Million mapped reads (FPKM) (Mortazavi et al., 2008). For the rest of the analysis, exploration, manipulation and visualization of the data generated by Cufflinks we used the R package cummeRbund (Goff et al., 2011). The analysis of the affymetrix ath1-121501 microarray data obtained for Arabidopsis was performed following the methodology described in Hahne et al. (2008). The R package affy (Gautier et al., 2004) from the Bioconductor project (Gentleman et al., 2004) was used for quality control, background correction, normalization with the RMA algorithm (Irizarry et al., 2003), and estimation of gene expression levels in the different conditions under study in this work. RESULTS AND DISCUSSION CONSTRUCTION OF GENE CO-EXPRESSION NETWORKS FOR Chlamydomonas,Physcomitrella,ANDArabidopsis In this analysis, significant co-expression patterns among genes in the transcriptomes of Chlamydomonas,Physcomitrella,and Arabidopsis have been determined analysing the massive amount of gene expression data (see Materials and Methods section) that covers a wide variety of physiological conditions and genotypes (Figure 1). These co-expression patterns were represented using three different gene co-expression networks, one for each species considered in this study representing key steps in the evolution of photosynthetic organisms. First, in order to remove noise, only those genes that exhibited significant changes in at least one comparison between a condition and its corresponding control, were selected from the studied transcriptomes. This was performed according to the standard approach used for microarray data and its adaptation to the analysis of RNA-seq data (Bullard et al., 2010). The logarithm of the gene expression level measured in FPKM was computed and, using the delta method, the variance of the log odds estimated. Differentially expressed genes were then selected combining this information with a fold-change criterion of two in the expression level with respect to the corresponding control. Next, the expression profiles of the differentially expressed genes were extracted from the gene expression data generated according to the workflow presented in Materials and Methods section.TheabsolutevalueofthePearson correlation coefficient between gene expression profiles was used as a measurement of the level of co-expression between the corresponding genes. Finally, it was necessary to establish a correlation threshold above which it was assumed that two genes are significantly coexpressed. For the rational selection of a gene correlation threshold a criterion that seeks the generation of a scale-free network with a high density was used. This criterion was chosen since most biological networks characterized so far exhibit this scale-free property (Barabasi and Albert, 1999)andbecauseahigh-density network facilitates the identification of patterns in the neighborhood of genes (Aoki et al., 2007). For each species, the correlation value for which the R2of the linear regression for the logarithmic transform of the node degree distribution presented a maximum while keeping an average of ∼20 neighbors per gene was determined. For Chlamydomonas reinhardtii the correlation threshold was 0.90, producing a R2of 0.9274 and a p-value of 3.397e-08 for the scale-free property. For Physcomitrella patens the correlation threshold was 0.94 producing a R2of 0.9333 and a p-value of 8.242e-08 for the scale-free property. For Arabidopsis thaliana www.frontiersin.org August 2013 | Volume 4 | Article 291 |5 Romero-Campero et al. Conserved light signaling in plants the correlation threshold was 0.90, producing a R2of 0.7827 and ap-value of 2.964e-04 for the scale-free property. According to the above methodology three co-expression networks in which edges between genes represent significant co-expression patterns in the analysed conditions were generated. The gene co-expression network for Chlamydomonas consisted of 8443 genes and 138,575 significant co-expression relationships. For Physcomitrella the corresponding network was constituted by 9080 genes connected by 518,209 significant co-expression relationships. Finally, the gene co-expression network for Arabidopsis represented 6204 genes and their 665,034 co-expression relationships. These networks were imported into the software package Cytoscape (Smoot et al., 2011)fortheirvisualization using the organic layout (Figure 2). The corresponding Cytoscape files containing the specification of each network can be downloaded from the link http://ackermann.cs.us. es/web_network/chlamy_physco_arabidopsis.zip for further exploration and analysis. IN THE EVOLUTIONARY HISTORY OF THE CO-Like GENE FAMILY FROM Chlamydomonas TO Arabidopsis BOTH ESSENTIAL AND SECONDARY GENES CAN BE IDENTIFIED The location of genes in co-expression networks and the size of their nearby neighborhood can be used to determine their FIGURE 2 | Gene co-expression networks for Chlamydomonas, Physcomitrella,andArabidopsis.In the networks, nodes represent genes and edges between nodes show significant co-expression relation between them. CrCO,PpCOLs, and AtCOLs and their neighborhoods (genes three edges apart) are highlighted with different colors. In the Chlamydomonas network (bottom left) CrCO and its neighborhood (red color) occupy a central position suggesting its role as an essential regulator in the Chlamydomonas transcriptome. In the Physcomitrella network (top center) PpCOL1,PpCOL4 and their neighborhoods (red and blue colors, respectively) occupy a central position, while PpCOL5 (green color) is located in the periphery of the network. In the Arabidopsis network (bottom right) AtCOLs genes (different colors) are widely spread over the network indicating their role in multiple biological processes. Functional conservation, specialization and divergence during the evolutionary line of COLs can be observed in their neighborhood. Consistently, CrCO,PpCOLs, and AtCOLs genes are co-expressed with genes involved in light response (ZEP, Zeaxanthin epoxidase), starch metabolism (SSS, solubale starch synthase), nitrogen metabolism (GS, cytosolic glutamine synthase)and response to chemical stimulus (ERF, ethylene response factor)asitis suggested by their nearby localization in the corresponding networks. relevance in the entire transcriptome of the corresponding organism in the conditions under study (Aoki et al., 2007). The three gene co-expression networks analysed in this work are scale-free networks. In this type of networks most nodes are connected with few others. Nevertheless, there exists a small set of nodes that are connected to a large number. These nodes are called hubs (Kleinberg, 1999)andplayakeyroleinthe functioning and information propagation in the corresponding networks (Barabasi and Albert, 1999). The hubs in each network were determined using the R package igraph (Csardi and Nepusz, 2006). Only genes with a score in the top 5% were considered hubs. The same criterion was used for the three networks. Scale-free networks are robust against random perturbations since these are most likely to hit a node with only a few neighbors and therefore to disrupt only a small portion of the network. Nevertheless, scale-free networks are fragile against perturbations affecting the hub nodes. Therefore, the mutation of a hub node affects a significant number of other nodes connected to it producing a cascade phenomenon that reaches a large part of the entire network (Wang and Chen, 2003). Employing this criterion we can identify central and essential genes as well as peripheral and secondary ones in the evolutionary line from the single gene CrCO in Chlamydomonas to the 17 AtCOL homologues in Arabidopsis (Figure 2 and Figure S2). In the Chlamydomonas gene co-expression network, CrCO appears located near the core of the network as a hub gene with more than 50 neighbors (Figure 2). This suggests that CrCO is an essential gene in the Chlamydomonas transcriptome under the conditions studied in this work. A CrCO overexpression or mutation would predictably produce a disruption in the functioning of all these neighboring genes that would propagate quickly throughout the network affecting a large number of other genes and the biological processes they are involved in. Therefore, any major change in the expression of CrCO could result in an extensive change in the Chlamydomonas phenotype. This effect has been reported previously (Serrano et al., 2009)wheresilencingofCrCO using antisense RNA was shown to be detrimental for algal growth. Moreover, overexpression of CrCO was reported to produce massive changes in cellular morphology, chlorophyll content, growth rate and cell-cycle regulation. Using sequence similarity programs from the Phytozome web portal, ten genes were identified in the Physcomitrella genome exhibiting high similarity (>40%) with the Arabidopsis CO gene (At5g15840) or the Chlamydomonas CrCO gene (g6302). These genes that contained at least one B-box and CCT domain, are identified in the current annotation version v1.6 of the Physcomitrella genome draft (Table 2). Three of these genes have been studied previously (Shimizu et al., 2004; Zobell et al., 2005) and named PpCOL1,PpCOL2,andPpCOL3 respectively. Based on the phylogeny and following the same nomenclature, the rest of the genes were named PpCOL4 to PpCOL10. Only three of these genes PpCOL1,PpCOL4,andPpCOL5 are present in the Physcomitrella gene co-expression network generated here. The rest of the genes either were not expressed in the samples considered in this study or did not exhibit significant co-expression levels with any other gene in the network. Frontiers in Plant Science | Plant Evolution and Development August 2013 | Volume 4 | Article 291 |6 Romero-Campero et al. Conserved light signaling in plants This suggests that these genes play their main role in other physiological conditions not studied in this work. PpCOL1 and PpCOL4 are found at the center of the network and although they cannot be regarded as hubs they are connected to other genes that are co-expressed with more than 300 genes. Both genes are placed close to each other suggesting that they may be involved in similar or related biological processes. This would confer redundancy and robustness to the network, so that the alteration in the expression pattern of one of these two genes could be counteracted by the normal function of the other. In this respect, the apparent multiple gene duplication of CrCO that gave rise to the ten-gene family of PpCOLs has produced redundant hub genes offering robustness to the transcriptome of Physcomitrella.TheCrCO descendant gene in the Physcomitrella network, PpCOL5, appears in the periphery of the gene coexpression network acting as a hub with 122 neighboring genes. Therefore, PpCOL5 seems to be a key gene involved in secondary biological processes in this Physcomitrella transcriptome since the genes of this cluster are somehow disconnected from the rest of the network (Figure 2). This may suggest that, although a major modification of the expression of PpCOL5 would not be lethal, it could massively disrupt the biological processes in which it is involved. Finally, in the Arabidopsis gene co-expression network, seven AtCOL genes out of a 17-member gene family (Table 2) have been identified. These are CO,COL1,COL2,COL5,COL7,COL9,and COL13. The other AtCOL genes were not significantly expressed in the conditions considered in this study possibly because they are involved in biological processes not studied here. These AtCOLs appear spread all over the network, suggesting that they are involved in a large number of biological processes covering an important part of the transcriptome. COL1,COL2,andCOL7 are closely arranged at the core of the network (red cluster in Figure 2). In fact, COL1 and COL2 are directly connected. The three of them have more than 150 neighbors each and could independently be considered hubs playing a key role in the function and transmission of information in the core of the network. Since these genes are located near each other they are probably involved in related biological processes. A major disruption in the expression pattern of one of these three genes could be easily compensated by the correct functioning of the other two. In fact, single mutations in either COL1 or COL2 show no significant effect in the plant (Ledger et al., 2001). Therefore, it could be predicted that, in order to observe any appreciable change in the phenotype a multiple mutant should be generated. This seems to indicate that the new events of gene duplication that took place between Physcomitrella and Arabidopsis have increased the number of redundant hub genes and, as a consequence, the robustness against external perturbations of its transcriptome. The rest of the AtCOLs in our network are distributed toward the periphery and present fewer neighboring genes. It could be inferred then that these genes are not essential hubs and could be involved in secondary processes. For example, a mutation in CO limits its effects to the capacity of the plant to flower in response to photoperiod (Rédei, 1962) and its overexpression results mainly in an early flowering phenotype (Putterill et al., 1995). Secondary phenotypic effects not observed in the mutant could arise due to CO overexpression assuming the function of other AtCOLs.In a similar situation, COL5 and COL9 are located nearby CO and their overexpression have been reported to produce a mild early flowering and late flowering phenotypes, respectively (Cheng and Wang, 2005; Hassidim et al., 2009). For COL13 its situation in the network is analogous but the effect of its mutation or overexpression has not been studied. A similar flowering phenotype could be predicted for this gene when overexpressed. THE GENE NEIGHBORHOODS OF CrCO,PpCOLs,ANDAtCOLs REVEAL CONSERVATION OF CO-EXPRESSED BIOLOGICAL PROCESSES ACROSS THE PLANT EVOLUTIONARY LINEAGE In order to study the different biological processes associated to the genes that constitute the evolutionary line from the single copy CrCO in Chlamydomonas to the 17 AtCOLs in Arabidopsis (Figure S2), their neighborhood in the corresponding networks was characterized. For this reason, genes at a distance of at least three, that is, those genes that can be reached using paths with three edges from the genes of interest, were selected. Next, GO (The Gene Ontology Consortium, 2000) term enrichment over the annotation of these genes using the web-based software tools GOrilla (Eden et al., 2009)andAgriGO(Du et al., 2010)was performed,takingtheentiregenomeasbackground.Inthecase of Arabidopsis we were able to use directly the GO annotation information available at the resource TAIR v.10. Nevertheless, for Physcomitrella and Chlamydomonas such information is not available due to the lack of previous experimental studies. Instead, for these species, GO annotation inferred from the information in Phytozome and Pfam (Punta et al., 2012)abouttheprotein families to which each gene-encoded protein belongs was used. First, for the GO enrichment in each one of the three networks, the union of the neighborhood of the corresponding COL genes under study, was considered. More specifically, for Chlamydomonas,Physcomitrella,andArabidopsis the set of genes used to perform GO enrichment were the neighbors at a distance of three of CrCO,anyofthePpCOLs and any of the AtCOLs, respectively. At a first glance, a high overlapping between the significantly enriched GO terms in the previous three sets of neighboring genes could be observed. This indicates that the different biological processes in which CrCO,somePpCOLs and AtCOLs are involved are largely conserved across the species evolution. Additionally, new biological processes appeared in the enrichment of the corresponding sets of genes in Physcomitrella and Arabidopsis that both species share in common (Table 3). This provides evidences of newly acquired influence over different processes in Physcomitrella that are afterwards conserved and expanded in Arabidopsis. One of the most significantly enriched GO terms in the regions under study in the three networks from Chlamydomonas to Arabidopsis is Response to light stimulus. Genes located in the neighborhood of CrCO,PpCOLs,andAtCOLs significantly include genes involved in photoperception such as Cryptochromes (CRY)(Cre06.g295200,Pp1s488_10V6,At4g08920), in redox processes regulated by light such as Zeaxanthin epoxidase (ZEP) (Cre02.g095750,Pp1s321_9V6,At5g67030)andCytochrome P450 (CYP)(Cre07.g325000,Pp1s281_82V6,At5g05690)andinthe protection against high light intensity such as Light Harvesting www.frontiersin.org August 2013 | Volume 4 | Article 291 |7 Romero-Campero et al. Conserved light signaling in plants Complexes (LHC)(Cre17.g740950,Pp1s628_3V6,At4g10340)of the photosystem II. This response to light has been previously reported for CrCO (Serrano et al., 2009), PpCOL1 (Zobell et al., 2005)andCO (Valverde et al., 2004). Therefore, this interspecies network analysis could successfully predict gene function using the significantly enriched GO terms observed in the corresponding modules around the genes of interest. The next conserved GO terms are involved in metabolism such as Starch metabolic process; Lipid metabolic process and Nitrogen compound metabolic process. Consistently, in the regions of the three networks, genes involved in starch biosynthesis such as Starch Synthases (SSS)(Cre16.g665800, Pp1s234_74V6,At4g18240), in starch hydrolysis such as Betaamylases (BAM)(Cre01.g044100,Pp1s317_42V6,At2g32290), in fatty acid synthesis such as Long Chain Acyl-CoA synthases (LACS) (Pp1s113_124V6,At2g47240)and3-ketoacyl-CoA synthase (KCS) (Cre07.g320550, Pp1s268_29V6, At4g34250), in fatty acid modification such as Fatty Acid Desaturases (FAD)(Cre13.g590500, Pp1s98_209V6,At3g15850), Glutamine Synthetase (GS)(g13061, Pp1s19_281V6,At5g37600)andNitrate Transporter (NTR) (g18260,Pp1s283_88V6,At4g18480) can be found. Recently, the influence of CrCO and CO over the metabolism of starch and carbohydrates has been described (Serrano et al., 2009;Ortiz et al., unpublished results). The effect of these genes over the metabolism of lipids and nitrogen compounds suggested in this work remains to be explored experimentally. This reflects the possible application of this type of analysis, employing the evolutionary line of gene co-expression networks, for gene function prediction. This way, when a gene with an unknown function is consistently co-expressed across its evolutionary lineage with the same well-characterized genes, a putative function could be assigned to it. Finally, several GO terms defined as Response to chemical stimulus and Developmental process, were identified that significantly and specifically appeared enriched both in the Physcomitrella and Arabidopsis networks. Although these GO term did not appear significantly represented in the neighborhood of CrCO in the Chlamydomonas network we were still able to identify several algal genes that showed a similar annotation, including Ethylene Response Factors (ERF) (g15714,Pps265_50V6,At1g53910)andAuxin Transporters (AUX) (Cre16.g680200,Pp1s167_10V6,At5g43700). Other GO term such as Response to jasmonate only presented conserved genes in Physcomitrella and Arabidopsis networks such as Jasmonate Insensitive (JIN)(Pp1s11_350V6,At1g32640)andJasmonate Responsive (JR)(Pp1s200_12V6,At3g16470) genes. In fact, it has been suggested that Chlamydomonas can respond to some hormonal stimuli and that this effect has stretched and become more complex as plants developed intricate structures and functions, reaching the complexity of the current plant hormone responses (Riaño-Pachon et al., 2008). These could be identified as the precursors of genes associated to hormone signaling and experiments could be designed to demonstrate their function. Therefore, another application of this cross species network analysis could be to predict the evolution of physiological processes and their importance in the diversification of key regulatory responses (such as hormone responses) between different species. If this is the case, when applied to a species transcriptome, implying the study of regulatory modules of many key factors in a natural population context, this analysis could be a very useful tool to understand speciation. GENE SPECIALIZATION AND THE ESTABLISHMENT OF CONNECTIONS BETWEEN DIFFERENT BIOLOGICAL PROCESSES IS REFLECTED IN THE EVOLUTIONARY HISTORY FROM CrCO TO AtCOLs Previously it was shown that, in the multi-gene families of Physcomitrella and Arabidopsis, different genes appear at distant positions in the network showing distinct levels of importance in the organization and information transmission of the corresponding networks. In order to determine if this influences the biological functions of genes, GO term enrichment in the neighborhood of each individual gene of interest in this work, was performed. In the Chlamydomonas genome a single CrCO copy involved in core biological processes can be identified: response to light stimulus, starch metabolic process, lipid metabolic process and nitrogen compound metabolic process. Attending to our analysis CrCO could also be marginally involved in pre-developmental processes and response to chemical stimuli (Table 4 and Figure 3). In the Physcomitrella network three COL genes were identified: PpCOL1,PpCOL4,andPpCOL5. GO term enrichment in each of the neighborhood of these genes show that each one has diverged to specialize in a particular biological process. PpCOL1 exhibits a high level of co-expression with LHC, SSS and fatty acid denaturizes. This suggests a specialization of this gene in response to light and carbon metabolism, consolidating a connection that was already established in Chlamydomonas.Onthe other hand, in the neighborhood of PpCOL4 and PpCOL5, genes associated with nitrogen metabolism, aging, development and response to chemical stimuli can be found. This suggests that these processes, that were secondary to CrCO in Chlamydomonas, became related to new COL genes and possibly more relevant during the evolution of plants into species like Physcomitrella. Duplication and diversification from CrCO produced genes such as PpCOL4 and PpCOL5 that covered the new challenges related to maintaining a multicellular organism in order to coordinate, among others, its developmental processes. Of special interest among the genes highly co-expressed with these two genes is the Physcomythrella PEBP, Phosphatidyl-Etanolamine Binding Protein (Pp1s32_140V6). This gene shows a high sequence similarity with the FT gene, the mobile florigen hormone that constitutes the main target of CO in Arabidopsis. Therefore, another possible application of our approach would be to identify true functional orthologous genes by determining whether or not their neighboring genes of known functions in the corresponding networks are evolutionarily conserved. In the Arabidopsis network, a considerable increase in the complexity and diversification of the function of the seven AtCOL genes included in this study is apparent. A group of genes comprising COL1,COL2,COL5,andCOL7 are co-expressed with genes whose function is similar to the core processes associated with CrCO, such as response to light, starch, lipid and nitrogen Frontiers in Plant Science | Plant Evolution and Development August 2013 | Volume 4 | Article 291 |8 Romero-Campero et al. Conserved light signaling in plants Table 4 | GO terms significantly enriched in CrCO,PpCOLs, and AtCOLs neighborhood in the gene co-expression networks. Chlamydomonas reinhardtii Physcomitrella patens Arabidopsis thaliana Gene GO term*Representative genes** Gene GO term*Representative genes** Gene GO term*Representative genes** CrCO (g6302) Response to light stimulus GO:0009416 (2.73 ×10−5) Cre06.g295200 Cre07.g325000 Cre17.g740950 (24) PpCOL1 (Pp1s371_27V6) Response to light stimulus GO:0009416 (2.61 ×10−3) Pp1s281_82V6 Pp1s628_3V6 (13) COL1 (At5g15850), COL2 (At3g02380) COL5 (At5g57660), COL7 (At1g73870) Response to light stimulus GO:0009416 (3.26 ×10−6) At4g08920 At5g05690 At1g44575 (24) Response to cold GO:0009409 (6.73 ×10−6) At1g12860 At5g57560 At5g20630 (17) Starch metabolic process GO:0005982 (2.98 ×10−5) Cre16.g665800 Cre01.g044100 (7) Carbohydrate metabolic process GO:0005975 (1.2×10−3) Pp1s234_74V6 Pp1s317_42V6 (12) Starch metabolic process GO:0005982 (1.94 ×10−3) At4g18240 At5g24300 At2g32290 (7) Lipid metabolic process GO:0006629 (4.83 ×10−3) Cre13.g590500 g14829 Cre06.g256750 (20) Fatty acid biosynthetic process GO:0006633 (3.74 ×10−2) Pp1s268_29V6 Pp1s7_269V6 Pp1s98_209V6 (5) Lipid metabolic process GO:0006629 (1.45 ×10−2) At2g47240 At4g34250 At3g15850 (20) Nitrogen compound metabolic process GO:0006807 (5.63 ×10−3) g13061 Cre03.g207250 g18260 (72) PpCOL4 (Pp1s36_238V6) Nitrogen compound metabolic process GO:0006807 (4.7×10−3) Pp1s283_88V6 Pp1s128_117V6 (44) Nitrogen compound metabolic process GO:0006807 (3.34 ×10−5) At5g37600 At5g35630 At4g18480 At3g09650 (18) Aging GO:0007568 (1.82 ×10−2) Pp1s1_163V6 (6) Response to auxin stimulus GO:0009733 (7.46 ×10−4) At1g29440 At4g38850 At1g08810 (13) Developmental process GO:0032502 (4.52 ×10−4) Cre06.g271100 Cre06.g251700 Cre06.g299300 (53) PpCOL5 (Pp1s26_5V6) Developmental process GO:0032502 (1.19 ×10−4) Pp1s10_24V6 Pp1s232_84V6 Pp1s38_120V6 (16) Response to Gibberellin stimulus (0.029) At2g37640 At1g74670 At1g71030 (13) CO (At5g15840) Response to Jasmodic acid stimulus GO:0009753 (6.73 ×10−5) At1g32640 At3g16470 At1g19640 (9) Response to temperature stimulus GO:0009266 (1.86 ×10−3) At2g45660 At3g27660 (8) Response to carbohydrate stimulus GO:0009743 (2.15 ×10−2) At3g50060 At5g67300 (6) Response to chemical stimulus GO:0042221 (2.9×10−2) g15714 Cre16.g680200 Cre02.g095750 (11) Response to chemical stimulus GO:0042221 (6.81 ×10−3) Pp1s167_70V6 Pp1s265_50V6 Pp1s175_16V6 (23) Defense response GO:0006952 (6.39 ×10−4) At4g31800 At1g42560 (14) COL9 (At3g07650) Circadian rhythm GO:0007623 (4.56 ×10−5) At2g46830 At5g61380 (5) Response to cold GO:0009409 (3.65 ×10−3) At2g21660 At3g08730 (5) (Continued) www.frontiersin.org August 2013 | Volume 4 | Article 291 |9 Romero-Campero et al. Conserved light signaling in plants Gutiérrez,R.A.,Lejay,L.V.,Dean, A., Chiaromonte, F., Shasha, D. E., and Coruzzi, G. M. (2007). Qualitative network models and genome-wide expression data define carbon/nitrogen-responsive molecular machines in Arabidopsis. Genome Biol. 8, R7. Hahne, F., Huber, W., Gentleman, R., and Falcon, S. (2008). Bioconductor Case Studies. New York, NY: Springer. doi: 10.1007/978-0-387-77240-0 Hassidim,M.,Harir,Y.,Yakir,E.,Kron, I., and Green, R. M. (2009). Overexpression of CONSTANS-LIKE 5 can induce flowering in shortday grown Arabidopsis.Planta 230, 481–491. doi: 10.1007/s00425-009- 0958-7 Heinz, S., Benner, C., Spann, N., Bertolino, E., Lin, Y. C., Laslo, P., et al. (2010). Simple combinations of lineage-determining transcription factors prime cisregulatory elements required for macrophage and B cell identities. Mol. Cell 38, 576–589. doi: 10.1016/j.molcel.2010.05.004 Imaizumi, T. (2010). Arabidopsis circadian clock and photoperiodism: time to think about location. Curr. Opin. Plant Biol. 13, 83–89. doi: 10.1016/j.pbi.2009.09.007 Irizarry, R. A., Hobbs, B., Collin, F., Beazer-Barclay, Y. D., Antonellis, K. J., Scherf, U., et al. (2003). Exploration, normalization, and summaries of high density oligonucleotide array probe level data. Biostatistics 4, 249–264. doi: 10.1093/biostatistics/4.2.249 Iyer-Pascuzzi, A. S., Jackson, T., Cui, H., Petricka,J.J.,Busch,W.,Tsukagoshi, H., et al. (2011). Cell identity regulators link development and stress responses in the Arabidopsis root. Dev. Cell 21, 770–782. doi: 10.1016/j.devcel.2011.09.009 Jackson, S. D. (2009). Plant responses to photoperiod. New Phytol. 181, 517–531. doi: 10.1111/j.1469-8137.2008.02681.x Jang, S., Marchal, V., Panigrahi, K. C.,Wenkel,S.,Soppe,W.,Deng, X. W., et al. (2008). Arabidopsis COP1 shapes the temporal pattern of CO accumulation conferring a photoperiodic flowering response. EMBO J. 27, 1277–1288. doi: 10.1038/emboj.2008.68 Jones,D.T.,Taylor,W.R.,and Thornton, J. M. (1992). The rapid generation of mutation data matrices from protein sequences. Comput. Appl. Biosci. 8, 275–282. Khanna, R., Kronmiller, B., Maszle, D. R., Coupland, G., Holm, M., Mizuno, T., et al. (2009). The Arabidopsis B-box zinc finger family. Plant Cell 21, 3416–3420. doi: 10.1105/tpc.109.069088 Kim, W. Y., Fujiwara, S., Suh, S. S., Kim, J., Kim, Y., Han, L., et al. (2007). ZEITLUPE is a circadian photoreceptor stabilized by GIGANTEA in blue light. Nature 449, 356–360. doi: 10.1038/nature06132 Kimura, M., and Ohta, T. (1974). On some principles governing molecular evolution. Proc. Natl. Acad. Sci. U.S.A. 71, 2848–2852. doi: 10.1073/pnas.71.7.2848 Kleinberg, J. (1999). Authoritative sources in a hyperlinked environment. J. ACM 46, 604–632. doi: 10.1145/324133.324140 Kropat, J., Hong-Hermesdorf, A., Casero, D., Ent, P., Castruita, M., Pellegrini, M., et al. (2011). A revised mineral nutrient supplement increases biomass and growth rate in Chlamydomonas reinhardtii.Plant J. 66, 770–780. doi: 10.1111/j.1365-313X.2011.04537.x Lajoie, M., Bertrand, D., and El- Mabrouk, N. (2010). Inferring the evolutionary history of gene clusters from phylogenetic and gene order data. Mol. Biol. Evol. 27, 761–772. doi: 10.1093/molbev/msp271 Lázaro, A., Valverde, F., Piñeiro, M., and Jarillo, J. A. (2012). The arabidopsis E3 Ubiquitin ligase HOS1 negatively regulates CONSTANS abundance in the photoperiodic control of flowering. Plant Cell 24, 982–999. Ledger, S., Strayer, C., Ashton, F., Kay, S. A., and Putterill, J. (2001). Analysis of the function of two circadian-regulated CONSTANSLIKE genes. Plant J. 26, 15–22. doi: 10.1046/j.1365-313x.2001.01003.x Li,L.,Huang,Y.,Xia,X.,andSun,Z. (2006). Preferential duplication in thesparsepartofayeastprotein interaction network. Mol. Biol. Evol. 23, 2467–2473. doi: 10.1093/molbev/msl121 Long,T.A.,Tsukagoshi,H.,Busch,W., Lahner,B.,Salt,D.E.,andBenfey,P. N. (2010). The bHLH transcription factor POPEYE regulates response to iron deficiency in Arabidopsis roots. Plant Cell 22, 2219–2236. doi: 10.1105/tpc.110.074096 Manfield, I. W., Devlin, P., Jen, C. H., Westhead,D.R.,andGilmartin, P. M. (2007). Conservation, convergence, and divergence of light-responsive, circadianregulated, and tissue-specific expression patterns during evolution of the Arabidopsis GATA gene family. Plant Physiol. 143, 941–958. doi: 10.1104/pp.106.090761 Marchler-Bauer, A., and Bryant, S. H. (2004). CD-Search: protein domain annotations on the fly. Nucleic Acids Res. 1, W327–W331. doi: 10.1093/nar/gkh454 Metzker, M. L. (2010). Sequencing technologies – the next generation. Nat. Rev. Genet. 11, 31–46. doi: 10.1038/nrg2626 Miller, R., Wu, G., Deshpande, R.,Vieler,A.,Gartner,K., Li, X., et al. (2010). Changes in transcript abundance in Chlamydomonas reinhardtii following nitrogen deprivation predict diversion of metabolism. Plant Physiol. 154, 1737–1752. doi: 10.1104/pp.110.165159 Mortazavi, A., Williams, B., McCue, K., Schaer, L., and Wold, B. (2008). Mapping and quantifying mammalian transcriptomes by RNA- seq. Nat. Methods 5, 621–628. doi: 10.1038/nmeth.1226 Müller, G. B. (2007). Evo–devo: extending the evolutionary synthesis. Nat. Rev. Genet. 8, 943–949. doi: 10.1038/nrg2219 Nasvall,J.,Sun,L.,Roth,J.R.,and Andersson, D. I. (2012). Real-time evolution of new genes by innovation, amplification, and divergence. Science 338, 384–387. doi: 10.1126/science.1226521 Ohno, S. (1970). Evolution by Gene Duplication.NewYork,NY: Springer. Oldham, M. C., Horvarth, S., and Geschwind, D. H. (2006). Conservation and evolution of gene coexpression networks in human and chimpanzee brains. Proc. Natl. Acad.Sci.U.S.A. 103, 17973–17978. doi: 10.1073/pnas.0605938103 Patterson, K., Cakmak, T., Cooper, A., Lager, I., Rasmusson, A. G., and Escobar, M. A. (2010). Distinct signalling pathways and transcriptome response signatures differentiate ammonium- and nitratesupplied plants. Plant Cell Environ. 33, 1486–1501. Punta, M., Coggil, P. C., Eberhardt, R. Y., Mistry, J., Tate, J., Bournsell, C., et al. (2012). The Pfam protein families database. Nucleic Acids Res. 40, D290–D301. doi: 10.1093/nar/gkr1065 Putterill, J., Robson, R., Lee, K., Simon, R., and Coupland, G. (1995). The CONSTANS gene of Arabidopsis promotes flowering and encodes a protein showing similarities to zinc finger transcription factors. Cell 80, 885–900. doi: 10.1016/0092- 8674(95)90288-0 Rédei, G. P. (1962). Supervital mutants of Arabidopsis.Genetics 47, 443–460. Riaño-Pachon, D. M., Correa, L. G. G., Trejos-Espinosa, R., and Mueller-Roeber, B. (2008). Green transcription factors: a Chlamydomonas overview. Genetics 179, 31–39. Robson, F., Costa, M. M., Hepworth, S. R., Vizir, I., Piñeiro, M., Reeves, P. H., et al. (2001). Functional importance of conserved domains in the flowering-time gene CONSTANS demonstrated by analysis of mutant alleles and transgenic plants. Plant J. 28, 619–631. Romero, J. M., and Valverde, F. (2009). Evolutionarily conserved photoperiod mechanisms in plants. Plant Signal. Behav. 4, 642–644. doi: 10.4161/psb.4.7.8975 Ruijter, J. M., Ramakers, C., Hoogaars, W.,Bakker,O.,VanDenHoff, M. J. B., Karlen, Y., et al. (2009). Amplification efficiency: linking baseline and bias in the analysis of quantitative PCR data. Nucleic Acids Res. 37, 45. doi: 10.1093/nar/ gkp045 Rustici, G., Kolesnikov, N., Brandizi, M.,Burdett,T.,Dylag,M.,Emam,I., et al. (2013). Arrayexpress update – trends in database growth and links to data analysis tools. Nucleic Acids Res. 41, D987–D990. Saitou, N., and Nei, M. (1987). The neighbor-joining method: a new method for reconstructing phylogenetic trees. Mol. Biol. Evol.4, 406–425. Samach, A., Onouchi, H., Gold, S. E.,Ditta,G.S.,Schwarz-Sommer, Z., Yanofsky, M. F., et al. (2000). Distinct roles of CONSTANS target genes in reproductive development of Arabidopsis. Science 288, 1613–1616. doi: 10.1126/science.288.5471.1613 Schliesky, S., Gowik, U., Weber, A. P. M., and Bräutigam, A. (2012). RNA-seq assembly – are we there yet? Front. Plant Sci. 3:220. doi: 10.3389/fpls.2012.00220 Serrano, G., Herrera-Palau, R., Romero,J.M.,Serrano,A., Coupland, G., and Valverde, F. (2009). Chlamydomonas CONSTANS and the evolution of plant photoperiodic signalling. Curr. Biol. 19, 359–368. doi: 10.1016/j.cub.2009.01.044 Shimizu, M., Ichikawa, K., and Aoki, S. (2004). Photoperiod-regulated expression of the PpCOL1 gene encoding a homolog of CO/COL proteins in the moss Physcomitrella patens.Biochem. Biophys. Res. Commun. 324, 1296–1301. doi: 10.1016/j.bbrc.2004.09.194 Simon, R., Igeno, I. M., and Coupland, G. (1996). Activation of floral meristem identity genes in Arabidopsis.Nature 384, 59–62. doi: 10.1038/384059a0 Frontiers in Plant Science | Plant Evolution and Development August 2013 | Volume 4 | Article 291 |16 Romero-Campero et al. Conserved light signaling in plants Smoot, M., Ono, K., Ruscheinski, J., Peng-Liang, W., and Ideker, T. (2011). Cytoscape 2.8: new features for data integration and network visualization. Bioinformatics 27, 431–432. doi: 10.1093/bioinformatics/btq675 Sueoka,N.,Chiang,K.S.,andKates, R. J. (1967). Deoxyribonucleic acid replication in meiosis of Chlamydomonas reinhardi. J. Mol. Biol. 25, 47–66. doi: 10.1016/0022-2836(67)90278-1 Tamura, K., Peterson, D., Peterson, N., Stecher, G., Nei, M., and Kumar, S. (2011). MEGA5: molecular evolutionary genetics analysis using maximum likelihood, evolutionary distance, and maximum parsimony methods. Mol. Biol. Evol. 28, 2731–2739. doi: 10.1093/molbev/msr121 The Gene Ontology Consortium. (2000). Gene ontology: tool for the unification of biology. Nat. Genet. 25, 25–29. Tiwari,S.B.,She,Y.,Chang,H.-C., Hou, Y., Harris, A., Ma, S. F., et al. (2010). The flowering time regulator CONSTANS is recruited to the FLOWERING LOCUS T promoter via a unique cis-element. New Phytol. 187, 57–66. doi: 10.1111/j.1469-8137.2010.03251.x Tohge, T., and Fernie, A. R. (2012). Co-expression and co-responses: within and beyond transcription. Front. Plant Sci. 3:248. doi: 10.3389/fpls.2012.00248 Trapnell, C. (2010). Transcript assemblyandquantificationbyRNA- seq reveals unannotated transcripts and isoform switching during cell differentiation. Nat. Biotechnol. 28, 511–515. doi: 10.1038/nbt.1621 Trapnell, C., Patcher, L., and Salzberg, S. (2009). TopHat: discovering splice junctios with RNA-seq. Bioinformatics 25, 1105–1111. doi: 10.1093/bioinformatics/btp120 Trapnell, C., Roberts, A., Goff, L., Pertea, G., Kim, D., Kelley, D. R., et al. (2012). Differential gene and transcript expression analysis of RNA-seq experiments with TopHat and Cufflinks. Nat. Protoc. 7, 562–578. doi: 10.1038/nprot.2012.016 Urzica,E.I.,Adler,L.N.,Page,M. D.,Linster,C.L.,Arbing,M.A., Casero, D., et al. (2012). Impact of oxidative stress on ascorbate biosynthesis in Chlamydomonas via regulation of the VTC2 gene encoding aGDP-L-galactose phosphorylase. J. Biol. Chem. 287, 14234–14245. doi: 10.1074/jbc.M112.341982 Usadel, B., Obayashi, T., Mutwil, M., Giorgi, F. M., Bassel, G. W., Tanimoto, M., et al. (2009). Co-expression tools for plant biology: opportunities for hypothesis generation and caveats. Plant Cell Environ. 32, 1633–1651. doi: 10.1111/j.1365- 3040.2009.02040.x Valverde, F. (2011). CONSTANS and the evolutionary origin of photoperiodic timing of flowering. J. Exp. Bot. 62, 2453–2463. doi: 10.1093/jxb/erq449 Valverde, F., Mouradov, A., Soppe, W., Ravenscroft, D., Samach, A., and Coupland, G. (2004). Photoreceptor regulation of CONSTANS protein in photoperiodic flowering. Science 303, 1003–1006. doi: 10.1126/science.1091761 Wahl, V., Ponnu, J., Schlereth, A., Arrivault,S.,Langenecker,T., Franke, A., et al. (2013). Regulation of Flowering by Trehalose-6- Phosphate Signaling in Arabidopsis thaliana.Science 339, 704–707. doi: 10.1126/science.1230406 Wang, X. F., and Chen, G. (2003). Complex networks: small-world, scale-free, and beyond. IEEE Circuits Syst. Mag. 3, 6–20. doi: 10.1109/MCAS.2003.1228503 Wenkel, S., Turck, F., Le Gourrierec, J., Samach, A., and Coupland, G. (2006). CONSTANS and the CCAATboxbindingcomplexshare a functionally important domain and interact to regulate flowering of Arabidopsis.Plant Cell 8, 2971–2984. Wheeler, D. L., Barrett, T., Benson, D.A.,Bryant,S.H.,Canese,K., Church, D. M., et al. (2005). Database resources of the national center for biotechnology information. Nucleic Acids Res. 33, D39–D45. Xiao,L.,Wang,H.,Wan,P., Kuang, T., and He, Y. (2011). Genome-wide transcriptome analysis of gametophyte development in Physcomitrella patens. BMC Plant Biol. 15:177. doi: 10.1186/1471-2229-11-177 Yu,J.W.,Rubio,V.,Lee,N.Y.,Bai, S., Lee, S. Y., Kim, S. S., et al. (2008). COP1 and ELF3 control circadian function and photoperiodic flowering by regulating GI stability. Mol. Cell 32, 617–630. doi: 10.1016/j.molcel.2008.09.026 Zemach,A.,McDaniel,I.E.,Silva, P., and Zilberman, D. (2010). Genome-wide evolutionary analysis of eukaryotic DNA methylation. Science 14, 916–919. doi: 10.1126/science.1186366 Zhang, J. (2003). Evolution by gene duplication: an update. Trends Ecol. Evol. 6, 292–298. doi: 10.1016/S0169-5347(03)00033-8 Zobell, O., Coupland, G., and Reiss, B. (2005). The family of CONSTANS- Like genes in Physcomitrella patens.Plant Biol. 7, 266–275. doi: 10.1055/s-2005-865621 Conflict of Interest Statement: The authors declare that the research was conducted in the absence of any commercial or financial relationships that could be construed as a potential conflict of interest. Received: 29 April 2013; accepted: 13 July 2013; published online: 05 August 2013. Citation: Romero-Campero FJ, Lucas- Reina E, Said FE, Romero JM and Valverde F (2013) A contribution to the study of plant development evolution based on gene co-expression networks. Front. Plant Sci. 4:291. doi: 10.3389/fpls. 2013.00291 This article was submitted to Frontiers in Plant Evolution and Development, a specialty of Frontiers in Plant Science. Copyright © 2013 Romero-Campero, Lucas-Reina, Said, Romero and Valverde.Thisisanopen-access article distributed under the terms of the Creative Commons Attribution License (CC BY). The use, distribution or reproduction in other forums is permitted, provided the original author(s) or licensor are credited and that the original publication in this journal is cited, in accordance with accepted academic practice. No use, distribution or reproduction is permitted which does not comply with these terms. www.frontiersin.org August 2013 | Volume 4 | Article 291 |17