Gene expression data from RNA sequencing (RNA-seq) experiments in A. thaliana, S. lycopersicum, P. patens, and C. reinhardtii were obtained from the Gene Expression Omnibus, and gene co-expression networks were inferred for each species from transcriptomic data as described in section “Materials and Methods.” All four networks are scale-free networks (Supplementary Figure S1) (Romero-Campero et al., 2013, 2016) and have comparable sizes in terms of number of nodes, but there are remarkable differences in the way they are connected (Table 1). The A. thaliana network contains more than twice as many edges than the others, the average degree of its nodes (average number of connections) is one order of magnitude higher and its average shortest path length (average number of nodes between two random nodes) is lower. Even though the number of genes of each species represented in the networks is similar, in some species they are more connected, possibly due to differences in their endogenous regulation and the availability of experimental data. For that reason, we decided to do every comparative analysis between the different species in relative terms.

To be able to compare the co-expression networks of the different species, we first identified the orthologous nodes in each of them using the OrthoMCL method (Li et al., 2003). Up to 17,053 groups of genes were obtained. Genes in the same group were considered orthologs or paralogs if they belonged to different or the same species, respectively. The four species were represented unequally, as both A. thaliana and S. lycopersicum genes were present in ca. 70% of the groups, while P. patens genes were found in little more than 50% of them, and only ca. 30% of the groups contained genes from C. reinhardtii (Figure 1B). This was already expected, given the evolutionary distance among these species and the genomic complexity of each one.

To assess the possible contribution of DELLA proteins to co-expression networks architecture, we created subnetworks based on reported DELLA interactors known to act as transcriptional regulators. First, we compiled a list of all published DELLA interactors (Supplementary Table S1), obtained their orthologs for the four species, and localized them in their respective networks. Since most of the interactions have been described for A. thaliana, the corresponding orthologs in the other species are only “putative interactors of the DELLA proteins” (PIDs), and the first neighbors of AtDELLA interactors and PIDs are their putative expression targets. Second, we built two different subnetworks using this information. The first one, called “Neighbors” subnetwork (abbreviated as AtNeigh, SlNeigh, PpNeigh, and CrNeigh), is composed of the DELLA interactors (or the corresponding PIDs) and their first neighbors (Figure 1C and Supplementary Table S2). The second one, called “Orthologs” subnetwork (abbreviated as AtOrtho, SlOrtho, PpOrtho, and CrOrtho), contains the orthologs of all the first neighbors of PIDs in all the species (Figure 1C and Supplementary Table S3). For a given species, the “Neighbors” subnetwork provides a good approximation to its actual DELLA-dependent transcriptome, while the “Orthologs” subnetwork represents the full landscape of potential transcriptional targets for DELLAs, since it includes orthologs of genes that are DELLA transcriptional targets in other species (Figure 2).

It is important to take into account a circumstance that affects the construction of subnetworks: OrthoMCL does not always retrieve orthologs for some of the genes, because either they do not exist in the other species, or the method does not provide high-confidence results. This results in a particular bias toward smaller subnetwork sizes with increasing phylogenetic distance (Table 1). However, the impact of this bias can be disregarded when analyzing relative parameters. Hence, regardless of the absolute sizes, we observed that the average degree in the Neighbor subnetworks increased dramatically in SlNeigh and AtNeigh with respect to their full networks (more than threefold and twofold, respectively), while this parameter did not change in PpNeigh, and it actually decreased in CrNeigh (Table 1). Similarly, the Orthologs subnetworks displayed an equivalent behavior as the Neighbors subnetworks: their diameter and average shortest path length decreased considerably more in SlOrtho and AtOrtho with respect to the full networks; and the same happened with the increase of the average degree. In summary, both subnetworks showed a higher compaction and interconnection of nodes in relative terms in the case of S. lycopersicum and A. thaliana compared with P. patens and C. reinhardtii, indicating that the putative interactors and targets of the DELLAs become more connected in those species presenting GA-regulated DELLAs.

A confirmation of the impact of GA regulation on the relevance of DELLA function is found in the analysis of neighborhood conservation. Figure 3A shows the percentage of genes with a significantly overlapping neighborhood in each comparison (see Materials and Methods). When comparing P. patens with the other species, there are no substantial differences between the full network and the Orthologs subnetwork. On the contrary, SlOrtho and AtOrtho contain a considerably higher proportion of genes with conserved neighborhood than their corresponding full networks (15% vs. 10%). Between S. lycopersicum and A. thaliana, the regulation of the putative DELLA targets is more conserved than for the network in general, so this group of genes seems to have a cohesive element in the two species.

Furthermore, we examined gene–gene co-expression values, as a measure of the conservation of individual edges. For every pair of linked genes in one species, if the corresponding orthologs are also linked in a second species, it is considered that gene–gene co-expression is conserved. Therefore, the calculation of conserved links between two subnetworks is a measure of functional conservation of a regulatory module. Interestingly, we observed that gene links between PpOrtho and SlOrtho were less conserved than in the full networks, and almost unaltered between PpOrtho and AtOrtho (Figure 3B). However, the gene–gene co-expression was three times more conserved between SlOrtho and AtOrtho than between their full networks (11% vs. 3.5%). In other words, these data are compatible with the proposition that the presence of GA-regulated DELLAs (in S. lycopersicum and A. thaliana) provides stronger links between transcriptional programs, not detected in an organism with GA-independent DELLAs (P. patens).

The efficiency of a transcriptional regulatory mechanism can be evaluated through two additional parameters in gene co-expression networks: shortest path length distribution and motif frequency. In network theory, average shortest-path length is defined as the average number of steps along the shortest paths for all possible pairs of network nodes. It is a measure of the efficiency of information propagation on a network, with a shorter average path length being more efficient (Vragovic et al., 2005). When we compared the distribution of shortest path lengths in full and Orthologs subnetworks, we observed a clear tendency toward shorter path lengths in the Orthologs subnetworks of organisms possessing DELLAs (S. lycopersicum, A. thaliana, and P. patens) compared with the situation in an organism without DELLAs (C. reinhardtii) (Figure 4).

Network motifs are small recurring patterns involving a few nodes that appear more frequently in biological networks than in randomized ones. They consist of a certain level of regulation which connects small sets of nodes with a particular topology. Motifs characterize a network, as some of them are useful for the regulation of determined functions, and thus conserved along evolution (Kashtan and Alon, 2005). After measuring the frequency of the eight common motifs composed of three and four nodes in the full networks, we found that there was no relative enrichment of any particular motif between species when comparing the full networks or the Orthologs subnetworks (Figure 5A). However, the AtOrtho, SlOrtho, and PpOrtho subnetworks displayed a clear enrichment in virtually every motif, compared with their respective full networks (Figure 5B). Given that the function of this sort of motifs is to allow coordinated expression of a group of genes with shared function (Alon, 2007), the increase in the proportion of small regulatory patterns among all the putative DELLA targets in species that do contain DELLAs indicates an increase in the complexity of gene regulation, in which DELLAs might mediate the coordination of transcriptional programs.

The results shown above suggest that the origin of DELLAs in land plants would be associated to an increase in the co-expression between genes that are putative targets of DELLA-interacting TFs, both in terms of size of the gene set and degree of the co-expression value. Therefore, DELLAs would have helped in the coordination of certain transcriptional circuits, and their recruitment to mediate GA signaling later in development would have further expanded their coordination capacity. To reveal the most likely functions ultimately regulated by DELLAs in the common ancestor of land plants, we carried out Gene Ontology (GO) analyses on each of the Neighbor subnetworks, with the idea that the terms shared by those in S. lycopersicum, A. thaliana, and P. patens could represent likely functions regulated by the ancestral DELLA proteins.

Not surprisingly, given the larger size of AtNeigh (Table 1), GO analysis rendered a much larger number of terms significantly enriched in this subnetwork, compared to those from the other three organisms (Supplementary Table S4). Terms referring to chloroplast function, such as plastid organization, photosynthesis, or pigment biosynthesis (including chlorophyll) were specifically enriched among the putative DELLA targets in A. thaliana only (Figure 6). This result might reflect functions whose regulation by DELLA has been acquired more recently, or it could simply be a bias of the analysis, caused by the big difference in size of the analyzed sets in the different species. On the contrary, the finding that terms comprised under general ‘response to stress’ were significantly over-represented in the subnetworks of the three land plants, but not C. reinhardtii, suggests that this function might have been the primary target of the regulation by ancestral DELLAs through their interaction with specific TFs.

The C. reinhardtii and A. thaliana networks were downloaded from the web resources of previous work (Romero-Campero et al., 2013, 2016). For the new networks, RNA-seq data were selected from equivalent experiments involving comparable tissues and environmental situations (Supplementary Table S5). The P. patens gene co-expression network was inferred from the RNA-seq data freely available from the Gene Expression Omnibus identified with accession numbers GSE19824, GSE33279, GSE36274, and GSE25237. The S. lycopersicum network was constructed based on the RNA-seq data identified with the accession numbers GSE45774, GSE64665, GSE64981, GSE68018, and GSE77340 in the Gene Expression Omnibus. In both cases, RNA-seq data was processed using the Tuxedo protocol (Trapnell et al., 2012) to obtain gene expression levels measured as FPKM. Briefly, short reads were mapped to the corresponding reference genome using Tophat, transcripts were assembled using Cufflinks and expression levels were computed using Cuffdiff. The Bioconductor R package cummeRbund (Goff et al., 2013) was used for subsequent analysis of the results generated by the Tuxedo protocol. In order to reduce noise in our analysis only genes that were detected as differentially expressed in at least one of the studies integrated in this work were considered. Differentially expressed genes were determined comparing each condition with the corresponding control within each study using a fold-change threshold of two. For each species, a matrix containing the expression levels of the selected genes was extracted. The Pearson correlation coefficient between every pair of gene expression profiles was computed using the cor function from the stats R package to generate a correlation matrix. Two genes were assumed to be co-expressed when the Pearson correlation coefficient between their expression profiles over the analyzed conditions was greater than 0.95. Following this criterion, the corresponding adjacency matrix was generated from the correlation matrix. Using the R package igraph¹ (Csardi and Nepusz, 2006), each network was constructed from its adjacency matrix and exported in gml formal for subsequent analysis.

The reference proteomes from A. thaliana TAIR10, S. lycopersicum iTAGv2.3, C. reinhardtii v5.5, and P. patens v3.3 were downloaded from Phytozome (Goodstein et al., 2012). From all the possible proteins from each locus tag only the longest protein was kept and assigned to its locus tag. These files were used to identify the orthologs among the four species with OrthoMCL (Li et al., 2003).

The networks were converted to SIF format and processed using the package igraph¹ (Csardi and Nepusz, 2006) made with R² (R Core Team, 2016). Only the edges between two non-identical nodes were conserved. If a given node was not identified in the proteome files, it was removed from the network. Afterward, components with fewer than seven elements were removed from the network to generate the complete network for each species. The orthologs for the set of manually curated DELLA interactors from A. thaliana were identified, and these nodes were selected from the complete networks. The first neighbors for all the selected nodes were identified and used to build a subnetwork. Finally, the orthologs on each species for all the genes in the previous subnetworks were identified and used to generate a new subnetwork for each species.

All networks were imported into the software package Cytoscape (Smoot et al., 2011) for their visualization using the Prefuse Force Directed layout.

The measures of network topology were calculated using both predefined and custom made functions. The gene–gene co-expression and neighborhood conservation were determined following the approach described by Netotea et al. (2014), using Fisher exact tests to check for statistical significance.

Gene Ontology analysis on Neigh subnetworks was made with AgriGO (Du et al., 2010), and represented with ReviGO (Supek et al., 2011).