Jasmien Vercruysse

Michiel Van Bel

Cristina M Osuna-Cruz

Shubhada R Kulkarni

V Eronique Storme

Hilde Nelissen

Nathalie Gonzalez

Dirk Inz
email: dirk.inze@psb.vib-ugent.

Klaas Vandepoele

Comparative transcriptomics enables the identification of functional orthologous genes involved in early leaf growth

Keywords: Arabidopsis thaliana, Zea mays, orthology, data integration, networks, comparative transcriptomics, leaf development

Leaf growth is a complex trait for which many similarities exist in different plant species, suggesting functional conservation of the underlying pathways. However, a global view of orthologous genes involved in leaf growth showing conserved expression in dicots and monocots is currently missing. Here, we present a genome-wide comparative transcriptome analysis between Arabidopsis and maize, identifying conserved biological processes and gene functions active during leaf growth. Despite the orthology complexity between these distantly related plants, 926 orthologous gene groups including 2829 Arabidopsis and 2974 maize genes with similar expression during leaf growth were found, indicating conservation of the underlying molecular networks. We found 65% of these genes to be involved in one-to-one orthology, whereas only 28.7% of the groups with divergent expression had one-to-one orthology. Within the pool of genes with conserved expression, 19 transcription factor families were identified, demonstrating expression conservation of regulators active during leaf growth. Additionally, 25 Arabidopsis and 25 maize putative targets of the TCP transcription factors with conserved expression were determined based on the presence of enriched transcription factor binding sites. Based on large-scale phenotypic data, we observed that genes with conserved expression have a higher probability to be involved in leaf growth and that leaf-related phenotypes are more frequently present for genes having orthologues between dicots and monocots than cladespecific genes. This study shows the power of integrating transcriptomic with orthology data to identify or select candidates for functional studies during leaf development in flowering plants.

Introduction

Zea mays and Arabidopsis thaliana are two important model organisms for monocots and dicots, respectively, which diverged 140-200 million years ago [START_REF] Chaw | Dating the monocot-dicot divergence and the origin of core eudicots using whole chloroplast genomes[END_REF]. Although they present numerous morphological and physiological differences, many developmental processes, such as leaf growth, show remarkable similarities [START_REF] Li | The developmental dynamics of the maize leaf transcriptome[END_REF][START_REF] Nelissen | Leaf growth in dicots and monocots: so different yet so alike[END_REF][START_REF] Yu | Transcriptome dynamics of developing maize leaves and genomewide prediction of cis elements and their cognate transcription factors[END_REF]. In both organisms, leaf growth is characterized by two major consecutive phases: a cell proliferation phase, during which cells divide, and a cell expansion phase, during which cells increase their volume. The interplay between cell division and cell expansion determines final leaf size [START_REF] Hepworth | Regulation of plant lateral-organ growth by modulating cell number and size[END_REF]. However, the direction of growth is different in both species. In Arabidopsis, the cellular growth pattern is dispersed, resulting in round mature leaves with reticulate veins, while in maize, growth is directed in a longitudinal fashion, generating narrow, elongated leaves with parallel veins [START_REF] Nelissen | Leaf growth in dicots and monocots: so different yet so alike[END_REF][START_REF] Nelson | Leaf vascular pattern formation[END_REF]. As soon as the leaf emerges from the shoot apical meristem, all cells from the primordium start dividing. After a few days, cells at the tip of the leaf cease division and start to elongate, marking the beginning of cell expansion. The leaf is then composed of both dividing and expanding cells during a so-called transition period, characterized by a cell-cycle arrest front moving from tip to base until most cells are no longer dividing and only expanding [START_REF] Andriankaja | Exit from proliferation during leaf development in Arabidopsis thaliana: a not-so-gradual process[END_REF][START_REF] Donnelly | Cell cycling and cell enlargement in developing leaves of Arabidopsis[END_REF][START_REF] Karidas | The tarani mutation alters surface curvature in Arabidopsis leaves by perturbing the patterns of surface expansion and cell division[END_REF][START_REF] Nelissen | Leaf growth in dicots and monocots: so different yet so alike[END_REF][START_REF] Pyke | Temporal and spatial development of the cells of the expanding first leaf of Arabidopsis thaliana (L.) Heynh[END_REF].

Conservation of the cell division and cell expansion phases between both species suggests that the underlying molecular mechanisms are preserved. These conserved mechanisms involve orthologous genes that evolved from a gene in a common ancestor through a speciation event. Nevertheless, different events including gene duplication and functional divergence result in complex many-to-many orthology in plants, where potentially only a subset of orthologues have conserved functions whereas the other orthologues acquired different functions or are differently expressed. Sorting out the true functional orthologues which do share similar biological functions can be accomplished by combining transcriptomic data with orthology information [START_REF] Movahedi | Comparative co-expression analysis in plant biology[END_REF][START_REF] Patel | BAR expressolog identification: expression profile similarity ranking of homologous genes in plant species[END_REF]. The resulting integrated network provides a powerful means to explore transcriptome-wide patterns of conservation and divergence between species [START_REF] Movahedi | Comparative network analysis reveals that tissue specificity and gene function are important factors influencing the mode of expression evolution in Arabidopsis and rice[END_REF]. Several studies have already described orthologous genes with a conserved expression pattern and/or function in two or more species [START_REF] Hefer | Comparative interrogation of the developing xylem transcriptomes of two wood-forming species: Populus trichocarpa and Eucalyptus grandis[END_REF][START_REF] Libault | An integrated transcriptome atlas of the crop model Glycine max, and its use in comparative analyses in plants[END_REF][START_REF] Mustroph | Cross-kingdom comparison of transcriptomic adjustments to low-oxygen stress highlights conserved and plant-specific responses[END_REF][START_REF] Narsai | Comparative analysis between plant species of transcriptional and metabolic responses to hypoxia[END_REF][START_REF] Rensink | Comparative analyses of six solanaceous transcriptomes reveal a high degree of sequence conservation and species-specific transcripts[END_REF]. For example, members of the GROWTH-REGULATING FACTOR (GRF) transcription factor (TF) family play key roles in leaf growth regulation and development and have been shown to interact with the coactivator GRF-INTERACTING FACTOR (GIF) both in Arabidopsis and in maize [START_REF] Kim | A transcriptional coactivator, AtGIF1, is involved in regulating leaf growth and morphology in Arabidopsis[END_REF][START_REF] Lee | The Arabidopsis GRF-INTERACTING FACTOR gene family performs an overlapping function in determining organ size as well as multiple developmental properties[END_REF][START_REF] Nelissen | Dynamic changes in ANGUSTIFOLIA3 complex composition reveal a growth regulatory mechanism in the maize leaf[END_REF][START_REF] Vercruyssen | ANGUSTIFOLIA3 binds to SWI/SNF chromatin remodeling complexes to regulate transcription during Arabidopsis leaf development[END_REF][START_REF] Zhang | Isolation and characterization of genes encoding GRF transcription factors and GIF transcriptional coactivators in Maize (Zea mays L.)[END_REF]. EXPANSINs function in both Arabidopsis and maize as cell wall-loosening proteins facilitating cell expansion and stimulating cell enlargement when overexpressed or applied exogenously [START_REF] Cho | Altered expression of expansin modulates leaf growth and pedicel abscission in Arabidopsis thaliana[END_REF][START_REF] Geilfus | Down-regulation of ZmEXPB6 (Zea mays b-expansin 6) protein is correlated with salt-mediated growth reduction in the leaves of Z. mays L[END_REF]. In these examples, the similar function of orthologous genes hints at a conservation of the molecular mechanisms that regulate leaf growth.

To identify the molecular processes active during leaf development, several genome-wide transcriptome studies addressing changes in gene expression during early leaf growth have been realized in Arabidopsis and maize [START_REF] Avramova | Drought induces distinct growth response, protection and recovery mechanisms in the maize leaf growth zone[END_REF][START_REF] Baerenfaller | Systems-based analysis of Arabidopsis leaf growth reveals adaptation to water deficit[END_REF][START_REF] Baute | Combined large-scale phenotyping and transcriptomics in maize reveals a robust growth regulatory network[END_REF][START_REF] Dubois | Time of day determines Arabidopsis transcriptome and growth dynamics under mild drought[END_REF][START_REF] Li | The developmental dynamics of the maize leaf transcriptome[END_REF][START_REF] Nelissen | The reduction in maize leaf growth under mild drought affects the transition between cell division and cell expansion and cannot be restored by elevated gibberellic acid levels[END_REF][START_REF] Skirycz | Pause-and-stop: the effects of osmotic stress on cell proliferation during early leaf development in Arabidopsis and a role for ethylene signaling in cell cycle arrest[END_REF]. These experiments provide a bird's-eye view on the expression of all genes possibly involved in leaf growth within a species. However, comparing gene expression data across multiple species requires accurate orthology information between the species. A network including orthology allows studying expression dynamics of different organisms at the gene level. Nonetheless, the current picture of all components involved in leaf growth shared between different plants is far from complete.

Here, we aimed at identifying common and divergent leaf growth-related regulatory processes in a monocot and a dicot plant. To do so, the expression conservation between Arabidopsis and maize was studied during early leaf development. We built an integrated network based on expression data linked with complex gene orthology information to study similarities between the two species, allowing to list potential leaf growth-regulating candidates. The level of conservation among Arabidopsis and maize TFs was also explored, indicating substantial conservation of transcriptional regulation.

Results

Similar expression patterns during leaf development in Arabidopsis and maize

To discover common or divergent expression patterns between Arabidopsis and maize during leaf development, we compared gene expression profiles in both species. Two publicly available gene expression data sets were used [START_REF] Andriankaja | Exit from proliferation during leaf development in Arabidopsis thaliana: a not-so-gradual process[END_REF][START_REF] Nelissen | The reduction in maize leaf growth under mild drought affects the transition between cell division and cell expansion and cannot be restored by elevated gibberellic acid levels[END_REF]. Both sets include time-course gene expression information from samples harvested during leaf development and correspond to the three main developmental phases: cell proliferation, cell expansion and the transition between these two phases (Figure S1a). For Arabidopsis, the expression information originated from six datapoints, in which two successive datapoints corresponded to one developmental phase. For maize, the expression information from samples harvested at nine datapoints was used, with four samples corresponding to the proliferating zone, three to the transitioning zone and two to the expanding zone. In total, expression levels were measured for 29 920 Arabidopsis and 39 323 maize genes of which, respectively, 4217 and 6495 differentially expressed (DE) genes (DEGs) were selected.

To identify major expression trends in both sets, genes were clustered based on their co-expression during leaf development.

Hierarchical clustering was used to group the DEGs in a minimal set of clusters reporting coherent expression profiles within each cluster and distinct profiles between clusters (Methods, Figure S1b, S1c). To facilitate statistical data analysis, only clusters with at least 50 genes were retained. For the Arabidopsis data set, initial clustering of the DEGs yielded 11 clusters, of which six containing more than 50 genes with distinct expression trends were retained: clusters A1 to A6 that represent 97% of the DEGs. For the maize data set, the same approach resulted in the selection of seven out of eight clusters (clusters Z1 to Z7), representing 99% of the DEGs. In both sets of clusters, three trends were visible, corresponding to a peak of expression during proliferation, during transition or during expansion (Figure S1d, S1e). Additionally, in the maize data set, a fourth expression trend was visible corresponding to expression peaks during proliferation and expansion, forming an inverse profile of the transition-specific expression trend. To identify molecular processes active during leaf growth, we performed GO enrichment on all Arabidopsis and maize gene clusters (Figure S2). Clusters with decreasing expression (A1 and Z1) were enriched for terms such as cell division and DNA replication, processes that occur during early leaf development [START_REF] Andriankaja | Exit from proliferation during leaf development in Arabidopsis thaliana: a not-so-gradual process[END_REF]. In A2, however, terms such as responses to stress (wounding, hypoxia), meristem development and leaf polarity determination were enriched, of which the latter process is known to occur in very early leaf development [START_REF] Stahle | YABBYs and the transcriptional corepressors LEUNIG and LEUNIG_HOMOLOG maintain leaf polarity and meristem activity in Arabidopsis[END_REF]. In clusters with a transition-specific profile (A3, A4, Z2, Z3 and Z4), different terms were found enriched in function of the cluster. Cluster A3 was enriched for terms related to cell division, whereas A4 was enriched for galactolipid and chlorophyll biosynthesis correlating with photosynthesis onset. Clusters Z2, Z3 and Z4 included genes leading to few enriched terms that suggested divergent biological processes such as root morphogenesis, fatty acid metabolic process and protein translation, respectively. Terms including photosynthesis and chlorophyll biosynthesis were found to be enriched in clusters with increasing expression during leaf development (A5, Z5 and Z6). This is to be expected since photosynthesis begins when leaves start to expand [START_REF] Andriankaja | Exit from proliferation during leaf development in Arabidopsis thaliana: a not-so-gradual process[END_REF][START_REF] Van Dingenen | Plants grow with a little help from their organelle friends[END_REF]. In contrast, A6 included terms related to plant defence [glucosinolate biosynthesis [START_REF] Wittstock | Glucosinolate breakdown in Arabidopsis: mechanism, regulation and biological significance[END_REF]], energy production (myo-inositol hexakisphosphate biosynthesis), auxin signalling and response to brassinosteroid. For Z7, no enriched terms were found.

To conclude, we mostly found similar expression trends for Arabidopsis and maize DEGs during leaf development with coinciding gene functions, suggesting that common pathways control leaf development in both species.

Integration of orthology and expression information identifies conserved orthologues expressed during leaf development

To identify conserved gene functions during leaf development in both species, gene orthology was determined and analysed. Orthologous genes were identified with the PLAZA integrative orthology method [START_REF] Van Bel | Dissecting plant genomes with the PLAZA comparative genomics platform[END_REF] that combines three orthology prediction methods: orthologous gene families inferred through sequence-based clustering with OrthoMCL (O); reconciled phylogenetic trees (T); and multispecies best hits and inparalogs (BHI) families (B) (Methods). Whereas the latter are related to a simple BLAST-based approach, OrthoMCL and phylogenetic trees are more advanced methods better capturing complex orthology relationships. For instance, using the B evidence only, three orthologues of Arabidopsis GRF3 could be identified, whereas including O and T evidence, twelve additional maize orthologues were identified (Figure 1a). To avoid missing any potential functionally conserved genes, we considered genes as orthologues when at least one type of evidence was found using PLAZA. At the genome-wide level, we found that 83.8% (22 987 genes) and 76.0% (29 874 genes), of the Arabidopsis and maize protein-encoding genes, respectively, have orthologues in the other species (Figure 1b,1c). Different levels of complexity were observed in the constructed orthology network composed of a total of 7336 separated subnetworks (or orthologous groups): one-to-one (31%), one-to-many (28%) and many-to-many orthology (41%; Figure 1d).

To identify orthology groups showing similar expression patterns in Arabidopsis and maize, we used the integrative orthology information to compare the expression clusters from Arabidopsis and maize. To do so, the overlap between two clusters was estimated by determining the number of Arabidopsis genes of a given cluster having orthologues in a given maize cluster and by calculating the significance of this overlap (Figure 2, Methods). Of the 1811 genes of the Arabidopsis A1 cluster, there were 1709 genes with maize orthologues, and for the maize Z1 cluster that includes 2649 genes, there were 2407 Arabidopsis orthologues. The comparison of these two clusters showed that 862 Arabidopsis genes had 924 maize orthologues and that this 50.4% overlap (862/1709) was significant (P < 0.001). The same comparisons were done for all clusters, and we observed a significant overlap for eight cluster comparisons (P < 0.001).

Overall, when including significant and non-significant overlaps, four distinct scenarios could be distinguished (Figure 2). First, for five cluster pairs (A1-Z1, A5-Z5, A5-Z6, A6-Z5 and A6-Z6), a significant overlap of orthologous genes showing similar expression patterns was observed. This significant overlap suggests that in both species, gene function during leaf development is probably conserved among orthologues with similar expression. Among these five cluster pairs, 217 one-to-one orthology pairs showing conserved expression were identified, including 69 groups with 110 unknown genes, which suggest 69 new maize and 41 new Arabidopsis candidates that might regulate leaf development in both species in a similar way (Table S1). Furthermore, the 217 one-to-one pairs with conserved expression corresponded to 91.2% of all DE one-to-one orthology pairs. Second, in three cases (A2-Z3, A3-Z6 and A4-Z6), there was a significant overlap between clusters with a divergent expression profile. For instance, in the A2-Z3 overlap (25 Arabidopsis and 22 maize genes), Arabidopsis genes showed high expression during proliferation, whereas for maize, the peak of expression occurred in the transition zone. The difference in expression was even more clear in the A3-Z6 and A4-Z6 significant overlaps. Clusters A3 and A4 showed high expression during transition, whereas in Z6 high expression occurred during cell expansion. This divergence could indicate that either the orthologous genes are involved in different processes or that they are part of similar processes, however regulated in a different manner during leaf development. One-to-one orthology was only found in the A3-Z6 and A4-Z6 overlaps, corresponding to, respectively, 0.8% (two groups) and 1.3% (three groups) of all DE one-to-one orthology pairs.

Third, seven of all cluster comparisons did not show any significant overlap although the respective clusters showed a conserved expression profile (A2-Z1, A3-Z2, A3-Z3, A3-Z4, A4-Z2, A4-Z3 and A4-Z4). Interestingly, six of these comparisons were between clusters that shared a transition-specific expression profile. The absence of significant overlap agrees with the poor GO enrichment overlap found for the five clusters (Figure S2). Finally, the other cluster comparisons between clusters with divergent profiles did not show any significant overlap.

Detection of conserved processes involved in leaf development

To identify the biological processes in which orthologues with conserved expression during leaf development take part, a GO enrichment analysis was done using the Arabidopsis shared Functional orthology by comparative transcriptomics 555 orthologues from the five significant overlaps A1-Z1, A5-Z5, A5-Z6, A6-Z5 and A6-Z6 that exhibited conserved expression (Figure 3).

A significant overlap of orthologues was observed for the A1-Z1 comparison. The enriched terms were related to different steps of the cell cycle including the M phase (cell plate formation and mitotic cytokinetic process) and the S phase (regulation of DNA duplication). For example, we observed 16 members of the MINI-CHROMOSOME MAINTENANCE (MCM) complex and 19 members of the histone family. Interestingly, RNA methylation suggesting modification of the RNA metabolism and histone H3-K9 methylation resulting in gene silencing [START_REF] Hu | Epitranscriptomic RNA methylation in plant development and abiotic stress responses[END_REF] were also found enriched in the orthologues showing a high expression during the proliferation phase.

The enriched terms in the A5-Z5 and A5-Z6 overlaps were mainly related to photosynthesis and energetic/carbohydrate metabolism. Within the genes part of the A5-Z5 overlap, 12 genes encoding glycoside hydrolases (including AtBGLU8 and AtBLGU9), which are involved in cell expansion [START_REF] Cosgrove | Growth of the plant cell wall[END_REF][START_REF] Minic | Plant glycoside hydrolases involved in cell wall polysaccharide degradation[END_REF], were present. Furthermore, 11 members of the glutathione S-transferase (GST) gene family (including AtGSTU17 and AtGSTU29) were also found in the GO enrichment analysis. In leaves, GST transcripts have been reported to accumulate in older rather than in younger tissues [START_REF] Gong | Expression of glutathione-S-transferase and its role in plant growth and development in vivo and shoot morphogenesis in vitro[END_REF], validating their allocation to the A5-Z5 overlap. In the shared Arabidopsis and maize orthologues from the A5-Z6 (389 Arabidopsis and 367 maize orthologues) overlap, aquaporins (11 genes all of which are either members of PLASMA MEMBRANE INTRINSIC PROTEIN or TONOPLAST MEMBRANE INTRINSIC PROTEIN families) were present. The aquaporin families have been shown to have an important role in expanding cells, causing turgor pressure by water uptake in Arabidopsis and maize [START_REF] Chaumont | Characterization of a maize tonoplast aquaporin expressed in zones of cell division and elongation[END_REF][START_REF] Ludevid | The expression pattern of the tonoplast intrinsic protein c-TIP in Arabidopsis thaliana is correlated with cell enlargement[END_REF].

In the A6-Z5 and A6-Z6 overlaps, instead of genes related to metabolism execution, we found genes related to regulatory processes such as hormone signalling, response to auxin and to The six selected Arabidopsis clusters are depicted horizontally (A1-A6, diamonds) and the seven selected maize clusters are ordered vertically (Z1-Z7, triangles), all of them sorted according to their visual trend. Clusters A1, A2 and Z1 have high expression in the early stages, corresponding to the proliferation phase, and a decreasing expression as leaf development advances (red). Clusters A3, A4, Z2, Z3 and Z4 show high expression during transition and low expression at the beginning and later stages of leaf development (blue). Clusters A5, A6, Z5 and Z6 have an upward gene expression trend coinciding with the onset of the expansion phase (green). Cluster Z7 shows high expression during proliferation and expansion, and low expression during the transition phase (yellow). The total number of genes in each cluster and the number of DEGs with an DE orthologue are shown under or right of the cluster name for Arabidopsis or maize, respectively. For each comparison between two clusters, the overlap is shown for both the Arabidopsis (upper number) and maize genes (lower number). Eight significant overlaps (P < 0.001) are shown in dark grey, and two significant overlaps (P < 0.01) are shown in light grey. At the bottom right, two representative clusters of the Arabidopsis (left) and the maize (right) DEG set are shown. Normalized gene expression is depicted on the x-axis according to the different datapoints on the y-axis, which corresponds to 8-13 DAS for the Arabidopsis data set and to distance to the leaf base for the maize data set. Ath., Arabidopsis thaliana; Zma., Zea mays; Nr., number; sp., species; DE, differentially expressed; DAS, days after stratification.

brassinosteroids, and secondary metabolite biosynthesis (phenylpropanoid metabolic process). For example, the AUXIN-INDOLE-3-ACETIC ACID (AUX/IAA) protein family, among which SHORT HYPOCOTYL 2 (SHY2) and SOLITARY ROOT (SLR) were found in these overlaps. Additionally, genes encoding eight members of the SMALL AUXIN UP RNA (SAUR) protein family, of which some are known to regulate plasma membrane H + -ATPase activity resulting in cell elongation, were found in the A6-Z6 overlap [START_REF] Spartz | SAUR inhibition of PP2C-D phosphatases activates plasma membrane H + -ATPases to promote cell expansion in Arabidopsis[END_REF]. Several genes related to cell wall biosynthesis were also found, such as LACCASE17 (LAC17) and CHITINASE PROTEIN-LIKE2 (CTL2), involved in lignin biosynthesis, and XYLOGLUCAN ENDOTRANSGLUCOSYLASE/HYDROLASE7 (XTH7) [START_REF] Berthet | Disruption of LACCASE4 and 17 results in tissuespecific alterations to lignification of Arabidopsis thaliana stems[END_REF][START_REF] S Anchez-Rodr Iguez | Chitinase-like1/pom-pom1 and its homolog CTL2 are glucan-interacting proteins important for cellulose biosynthesis in Arabidopsis[END_REF].

In Arabidopsis, although the overlaps including A5 or A6 genes showed a general increased expression during leaf growth, it was obvious that the enriched categories are different in both overlaps with A5 containing genes related to metabolism while A6 contains genes that are more related to regulatory processes (hormone signalling) and secondary metabolism. Interestingly, in maize, these two types of categories seem to be distributed in both clusters Z5 and Z6.

With this analysis, we were able to identify 1967 Arabidopsis genes with 2050 maize orthologues showing conserved expression and function, being part of a significant overlap. Conserved biological processes related, but not limited, to the cell cycle were found for 862 Arabidopsis and 924 maize orthologous genes with proliferation-specific expression in the A1-Z1 overlap. Interestingly, RNA and histone methylation were also enriched. A large part of the 1105 Arabidopsis and 1126 maize orthologous genes with expansion-specific expression (A5-Z5, A5-Z6, A6-Z5 and A6-Z6) were clearly related to photosynthesis and carbohydrate metabolism, but also to hormonal signalling. To conclude, we were able to show large functional conservation among orthologues in dicots and monocots with conserved expression. Functional orthology by comparative transcriptomics 557

Phenotypic analysis of Arabidopsis genes showing conserved and divergent expression profiles

To further explore the impact of conserved expression on gene function during leaf development, publicly available phenotypic data were mapped onto the orthology and expression information we combined. TAIR [START_REF] Berardini | The Arabidopsis information resource: making and mining the "gold standard" annotated reference plant genome[END_REF] and RARGE II [START_REF] Akiyama | RARGE II: an integrated phenotype database of Arabidopsis mutant traits using a controlled vocabulary[END_REF] were consulted to generate a list of 738 Arabidopsis genes that if mutated, lead to a leaf-related phenotype such as an increased leaf rosette size (Methods). The genes in A1, A5 and A6 that had a significant overlap were used in this analysis. For these clusters, three groups were defined as follows: (i) a conserved group containing Arabidopsis genes with at least one orthologue showing conserved expression, (ii) a divergent group in which Arabidopsis genes only had orthologues with divergent expression and (iii) a group of Arabidopsis genes lacking maize orthologues, thus representing species-specific genes. For example, cluster A5 included 119 genes with no maize orthologues and 756 which did have a DE maize orthologue. Of these 756 genes, 673 belonged to the conserved group and 83 to the divergent groups (Table S2). Phenotypic data were mapped onto these three groups for the three clusters. Globally, 4.0% (738/18 369) of the Arabidopsis genes lead to a leaf phenotype when mutated, while for A5 conserved and divergent groups, percentages of, respectively, 6.24% (42 genes; P-value 2.92EÀ02) and 7.23% (six genes; Pvalue 0.53) were found. For the A5 species-specific group, 3.36% of leaf-related phenotypes were found (four genes; P-value 0.87). These results suggest that DE A5 genes with DE maize orthologues during leaf development have a higher possibility to cause a leafrelated phenotype when mutated. For A1 and A6, the percentage of leaf-related phenotypes in the conserved group (3.58% and 2.86%, respectively) and the divergent group (4.84% and 5.88%, respectively) was greater as compared to the species-specific group (2.94% and 2.47%, respectively).

Over all groups, there were only nine genes that result in an increased leaf size when mutated, of which three could be confirmed [START_REF] Brauner | Exaggerated root respiration accounts for growth retardation in a starchless mutant of Arabidopsis thaliana[END_REF][START_REF] Cheng | Functional association of Arabidopsis CAX1 and CAX3 is required for normal growth and ion homeostasis[END_REF][START_REF] Gong | Expression of glutathione-S-transferase and its role in plant growth and development in vivo and shoot morphogenesis in vitro[END_REF][START_REF] Storozhenko | c-Glutamyl transpeptidase in transgenic tobacco plants. Cellular localization, processing, and biochemical properties[END_REF]. In contrast, 49 genes were found that decreased leaf size when mutated, of which 13 could be confirmed [START_REF] Chiang | Isolation of the Arabidopsis GA4 locus[END_REF][START_REF] Hobbie | The axr4 auxin-resistant mutants of Arabidopsis thaliana define a gene important for root gravitropism and lateral root initiation[END_REF][START_REF] Kim | Pathogen-induced Arabidopsis WRKY7 is a transcriptional repressor and enhances plant susceptibility to Pseudomonas syringae[END_REF][START_REF] Lintala | Structural and functional characterization of ferredoxin-NADP + -oxidoreductase using knock-out mutants of Arabidopsis[END_REF][START_REF] Liu | Loss-of-function mutation in SCY1 triggers chloroplast-to-nucleus retrograde signaling in Arabidopsis thaliana[END_REF][START_REF] Mattioli | Modulation of intracellular proline levels affects flowering time and inflorescence architecture in Arabidopsis[END_REF][START_REF] Ohkama-Ohtsu | Characterization of the extracellular c-glutamyl transpeptidases, GGT1 and GGT2, in Arabidopsis[END_REF][START_REF] Peng | LOW PSII ACCUMULATION1 is involved in efficient assembly of photosystem II in Arabidopsis thaliana[END_REF][START_REF] Simkov A | Integration of stress-related and reactive oxygen species-mediated signals by Topoisomerase VI in Arabidopsis thaliana[END_REF][START_REF] Uematsu | Increased fructose 1, 6-bisphosphate aldolase in plastids enhances growth and photosynthesis of tobacco plants[END_REF][START_REF] Xia | The Arabidopsis eukaryotic translation initiation factor 3, subunit F (AteIF3f), is required for pollen germination and embryogenesis[END_REF][START_REF] Xiao | DAC is involved in the accumulation of the cytochrome b6/f complex in Arabidopsis[END_REF][START_REF] Yoshizawa | BPG3 is a novel chloroplast protein that involves the greening of leaves and related to brassinosteroid signaling[END_REF]. For both phenotypes, a larger number of genes were found in the conserved expression groups compared with the divergent expression groups, suggesting that expression conservation is an important property to identify genes involved in leaf development (Table S3).

In the A5 conserved group, five genes causing an increased leaf size through overexpression, corresponded to AT1G04680, GLU-TAMATE:GLYOXYLATE AMINOTRANSFERASE1 (GGT1), 1-AMI-NOCYCLOPROPANE-1-CARBOXYLIC OXIDASE (ATACO2), RESPONSIVE TO DESICCATION22 (RD22) and STARCH-FREE1 (STF1). Apart from GGT1 and STF1, which are known to decrease leaf size when knocked out [START_REF] Brauner | Exaggerated root respiration accounts for growth retardation in a starchless mutant of Arabidopsis thaliana[END_REF][START_REF] Gong | Expression of glutathione-S-transferase and its role in plant growth and development in vivo and shoot morphogenesis in vitro[END_REF], these genes were only described in RARGE II and are thus interesting candidates for further characterization. Of these genes, we know that at least one maize orthologue showed conserved expression during leaf development and could thus also be involved in leaf size regulation. Additionally, in the A6 conserved group, TF DNA-BINDING WITH ONE ZINC FINGER2 (DOF2) was found to increase leaf size when overexpressed and is interestingly known to be expressed at the tip of the leaf in very early developmental stages [START_REF] Gardiner | Expression of DOF genes identifies early stages of vascular development in Arabidopsis leaves[END_REF]. With this analysis, we conclude that for orthologues with conserved expression, a significant fraction is associated with known leaf phenotypes.

Apart from conservation analysis, also cases suggesting expression divergence of orthologues were investigated. For example, in the A5 divergent group, we identified GIBBERELLIN 3-OXIDASE 1 (GA3ox1), considered the last and rate-limiting enzyme in the gibberellin biosynthetic pathway leading to the production of bioactive gibberellins, that causes a semi-dwarfed phenotype in Arabidopsis when knocked out [START_REF] Chiang | Isolation of the Arabidopsis GA4 locus[END_REF][START_REF] Fleet | Overexpression of AtCPS and AtKS in Arabidopsis confers increased ent-kaurene production but no increase in bioactive gibberellins[END_REF]. Together with the expression of GA3ox1, which increases as leaf growth advances with a peak during cell expansion, the ga3ox1-3 mutant phenotype demonstrates a function for this gene in the control of leaf expansion [START_REF] Davies | Plant Hormones: Biosynthesis, Signal Transduction, Action![END_REF][START_REF] Mitchum | Distinct and overlapping roles of two gibberellin 3oxidases in Arabidopsis development[END_REF]. In contrast, the maize GA3ox1 orthologue exhibited the inverse expression pattern with a peak in expression during cell proliferation and a decrease in expression as leaf growth advances coinciding with the accumulation of gibberellin at the transition zone, controlling the shift towards cell expansion at this boundary [START_REF] Nelissen | A local maximum in gibberellin levels regulates maize leaf growth by spatial control of cell division[END_REF]. Although both orthologues may metabolically carry out the same biochemical function, the regulation of this function is different in Arabidopsis and maize due to the spatial organization of the leaf. In total, we could identify eight Arabidopsis genes of which four, three and one belonged to the A1, A5 and A6 divergent groups, respectively, that had divergent expression patterns compared with their maize orthologues and that exhibited an increase (one gene) or decrease (seven genes) in leaf size when mutated (Table S3). There were only four genes in the A5 (two genes) and A6 (two genes) groups with no orthologues that displayed an increase (one gene in A6NO) or decrease (two genes in A5NO and one gene in A6NO) in leaf size when mutated, suggesting an Arabidopsis-specific leaf-related function. These results reveal the potential to identify possible candidates that are either Arabidopsis-specific or gene candidates that have an altered expression profile, suggesting a change in regulation or function.

Visualization of the integrated orthologous expression network

To study the complexity of the genome-wide orthologous relations between clustered Arabidopsis and maize DEGs, a visual representation was generated using Cytoscape [START_REF] Shannon | Cytoscape: a software environment for integrated models of biomolecular interaction networks[END_REF]. The complexity of the network was limited by only featuring the orthologous relations between DEGs of clusters A1-A6 and Z1-Z7 (Supplemental network file).

The filtered network included 2371 Arabidopsis and 2675 maize DEGs, visualized as nodes, that were connected through 4829 edges, representing orthology. Genes were coloured according to their expression profile trend. Additionally, edges were coloured when the orthologues were part of a significant overlap. This approach allowed studying transcriptional conservation between both species at both the gene level and the genome-wide level.

The network was divided into 1487 orthologous groups that can be seen as (parts of) gene families in which node colour visualizes conserved/divergent expression. There were 421 orthologous groups with conserved expansion-specific expression, of which 67.2% had a one-to-one orthology relation. Conserved proliferation-specific expression was observed for 501 orthologous groups, of which 64.7% had one-to-one orthology relation. Only four groups containing ten genes were found with conserved transitionª 2019 The Authors. Plant Biotechnology Journal published by Society for Experimental Biology and The Association of Applied Biologists and John Wiley & Sons Ltd., 18, 553-567 specific expression, of which three displayed one-to-one orthology. These data therefore revealed that 62.7% of the integrated orthology network showed expression conservation among orthologues. In contrast, 37.3% (561 groups) of the groups contained orthologues with divergent expression. Within these groups, only 161 corresponded to one-to-one orthology (28.7%), showing that orthologues with divergent expression have more complex orthology compared with completely conserved expression orthology groups. It seems thus more likely to observe orthologues with divergent expression in complex orthology groups with many-tomany orthology, suggesting neofunctionalization and subfunctionalization within such group. For example, among the 197 groups displaying one-to-two orthology (one Arabidopsis gene with two maize orthologues), for 59 groups one maize gene showed divergent expression (Table S4), suggesting a different role for the two orthologues. For instance, for Arabidopsis ALCOHOL DEHYDROGENASE (ADH), with high expression during proliferation, only ZmADH1 had conserved expression, whereas ZmADH2 had high expression during expansion. This example illustrates the duplication of a progenitor gene, after which the expression, and possibly the function, of ZmADH2 could have deviated [START_REF] Dennis | Molecular analysis of the alcohol dehydrogenase 2 (Adh2) gene of maize[END_REF].

Next, the conservation between DEGs acting as TFs or involved in chromatin-remodelling (CR), and their DE orthologues, which are expected to regulate the expression of numerous downstream genes, was analysed. There were 467 orthologues (encoding 211 Arabidopsis TFs and 256 maize TFs), which formed 106 orthologous groups corresponding to 47 TF and transcriptional regulator (TR) families (Table S5). These 47 families could be divided into four groups based on the number of orthologues with conserved expression: complete conservation; more than 50% conservation; less than 50% conservation; and complete divergence (Table S6). There were 19 TF and TR families, including 45 Arabidopsis genes and 56 maize orthologues, showing complete expression conservation. Among these families, two TF and TR groups were noticeable: a first group was related to chromatinmediated transcription regulation [among which HIGH MOBILITY GROUP (HMG) genes, GNAT family members, SNF2 genes, SET genes and the ULTRAPETALA1 (ULT1) TF], while a second group consisted of more specific regulators of gene expression such as YABBY, GRF, E2F and NAC TF families (Figure 4). YABBY TFs are involved in the determination of the adaxial/abaxial leaf polarity and correspondingly, all DE TFs of this family showed high expression during proliferation [START_REF] Stahle | YABBYs and the transcriptional corepressors LEUNIG and LEUNIG_HOMOLOG maintain leaf polarity and meristem activity in Arabidopsis[END_REF]. The major cell-cycle regulators, the E2F/DP TFs, which regulate the G1-to-S phase transition [START_REF] Berckmans | Transcriptional control of the cell cycle[END_REF][START_REF] De | Control of proliferation, endoreduplication and differentiation by the Arabidopsis E2Fa-DPa transcription factor[END_REF], and the GRF TFs [START_REF] Gonzalez | Leaf size control: complex coordination of cell division and expansion[END_REF] involved in the regulation of cell proliferation, had similar expression profiles corresponding to high expression during proliferation.

We observed 15 gene families, among which members of the TEOSINTE BRANCHED1/CYCLOIDEA/PCF (TCP) TF family, for which the expression conservation was higher than 50% (Figure S3). In Arabidopsis, class I and class II work antagonistically by activating and suppressing growth, respectively, balancing leaf developmental processes [START_REF] Nicolas | The role of TCP transcription factors in shaping flower structure, leaf morphology, and plant architecture[END_REF]. TCP orthologues with conserved expression included the class II AtTCP3, AtTCP4 and ZmTCP23, and class I AtTCP21 and ZmTCP4 with high expression during expansion (Figure 5a). Divergent expression was observed in the group comprising AtTCP19, ZmTCP19 and ZmTCP21. AtTCP19 and ZmTCP19 had conserved expression, whereas AtTCP19 and ZmTCP21 differed in expression, which identifies ZmTCP19 as the most likely functional orthologue of AtTCP19 (Figure 5a). Other TF families with more than 50% of the orthologues with conserved expression included, but were not limited to, the CCAAT, AUX/IAA and the C2H2 families.

Members of the MYB and bHLH family were categorized in the group with less than 50% conserved expression between orthologues. The rather low degree of conservation, 36.36% and 41.67%, respectively, can be explained by multiple duplication events in both species followed by functional divergence [START_REF] Carretero-Paulet | Genome-wide classification and evolutionary analysis of the bHLH family of transcription factors in Arabidopsis, poplar, rice, moss, and algae[END_REF][START_REF] Dubos | MYB transcription factors in Arabidopsis[END_REF][START_REF] Feller | Evolutionary and comparative analysis of MYB and bHLH plant transcription factors[END_REF]. Only three TF families, i.e. LIM, TAZ and BRI1-EMS-SUPPRESSOR 1 (BES1), included orthologues lacking expression conservation, suggesting that these proteins might operate differently during leaf development in both species.

We conclude that TF families playing key roles in cellular processes required for leaf growth frequently have similar expression profiles in both conserved and divergent orthology groups, indicating that they might also regulate the expression of similar target genes in both organisms.

Integration of transcription factor binding site information to identify conserved regulators involved in leaf development

The presence of enriched TF binding sites (TFBSs) in the promoters of orthologues with conserved expression suggests the preservation of the mode of action of the TFs and their targets during leaf growth in dicots and monocots. Therefore, we investigated whether known Arabidopsis TFBSs linked with cell cycle (E2F) or leaf maturation (TCP) were overrepresented in the promoters of Arabidopsis DEGs and maize orthologues with conserved expression. The studied binding sites were previously determined using protein binding microarrays [START_REF] Weirauch | Determination and inference of eukaryotic transcription factor sequence specificity[END_REF]. First, the proof of concept was established using the well-characterized Arabidopsis E2F TFs known to bind the degenerative WTTSSCSS motif and of which many Arabidopsis targets are already known and described [START_REF] Vandepoele | Genome-wide identification of potential plant E2F target genes[END_REF]. One specific instance of the degenerative E2F TFBSs, TTTCCCGC, was enriched in the promoters (500 bp upstream from the start codon) of A1 (47 genes, enrichment fold 2.04, P-value 2.2EÀ06) and Z1 (71 genes, enrichment fold 2.36 and P-value 1.26EÀ11; Figure S4), both containing genes with high expression during proliferation. Within this pool of genes, 12 Arabidopsis genes had 16 maize orthologues, indicating the presence of conserved E2F target genes (Table S7). Among these common targets were MCM2 and MCM5, involved in chromatin organization, as well as REPLICATION PROTEIN A 2B (RPA2), involved in DNA replication, repair and recombination [START_REF] Elmayan | Arabidopsis RPA2: a genetic link among transcriptional gene silencing, DNA repair, and DNA replication[END_REF]. These results were benchmarked using two previously published studies that identified E2F-regulated genes [START_REF] Naouar | Quantitative RNA expression analysis with Affymetrix Tiling 1.0R arrays identifies new E2F target genes[END_REF] and E2F-bound genes [START_REF] Verkest | A generic tool for transcription factor target gene discovery in Arabidopsis cell suspension cultures based on tandem chromatin affinity purification[END_REF]. The overlap between both studies resulted in 393 possible E2F targets. All 12 Arabidopsis genes were confirmed as previously described E2F targets, which validates our approach to identify putative, new TF targets in both species jointly.

To identify additional regulators and target genes involved in leaf development, the conservation of the known TCP TFBSs in both Arabidopsis and maize DEGs was studied. To do so, eight experimentally identified Arabidopsis TFBSs representing three class I TCPs (TCP16, TCP19 and TCP20) and five class II TCPs (TCP2, TCP3, TCP4, TCP5 and TCP24) were employed (Figure 5b; [START_REF] Weirauch | Determination and inference of eukaryotic transcription factor sequence specificity[END_REF]. After motif enrichment analysis, two TFBSs representing TCP3 and TCP24, both classified as JAW-TCPs responsible for leaf curvature [START_REF] Efroni | A protracted and dynamic maturation schedule underlies Arabidopsis leaf development[END_REF] Functional orthology by comparative transcriptomics 559 A1 with proliferation-specific expression (Table S8). Among the set of 229 A1 and 272 A5 genes with an enriched TCP TFBS in their promoters, 34 genes shared 34 orthologues from a total of 279 Z5 genes, of which 25 Arabidopsis genes and 25 maize orthologues showed conserved expansion-specific expression (Figure 5c). We can thus consider these 68 genes (i) to be differentially expressed during leaf growth with a peak in expression during cell proliferation or expansion, (ii) to be In each panel, the top image represents the integrated orthology network in which Arabidopsis TFs (diamonds) and maize orthologues (triangles) are linked. The colour of the nodes corresponds to an increasing (green) or a decreasing (red) expression trend during leaf development. The thickness of the edge corresponds to the number of evidences found for that particular orthology link (B: best hits and inparalogs families, O: OrthoMCL and T: reconciled phylogenetic trees). The colour of the edges represents a significant orthology overlap between the respective Arabidopsis and maize genes. The heat map in each panel represents the normalized expression values of each DE TF throughout leaf development. Ath., Arabidopsis thaliana; zma., Zea mays;DAS, days after stratification.

orthologues in Arabidopsis and maize, and (iii) to have a TCP binding site in their first 500 bp that is enriched in their respective clusters. From these three observations, 68 genes could thus be identified as putative TCP TF targets. Among the set of targets, LOX2, a known target of TCP4 and TCP20 [START_REF] Danisman | Arabidopsis class I and class II TCP transcription factors regulate jasmonic acid metabolism and leaf development antagonistically[END_REF] and its maize ortholog, GRMZM5G822593, were identified as TCP3 and/or TCP24 targets (Figure 5c).

In conclusion, through the integration of conserved TFs active during leaf development and TF binding site information, detailed gene regulatory subnetworks could be identified offering new insights into target genes for specific TFs in both species. For the E2F TFs, we could show strong conservation of the TFBSs in both Arabidopsis and maize, while for the TCP family, the TFBS conservation was not that pronounced.

Discussion

Leaf growth has been studied extensively both in Arabidopsis and in maize and there are remarkable similarities, at cellular and molecular levels, between the seemingly very different outcomes (for a review, see [START_REF] Nelissen | Leaf growth in dicots and monocots: so different yet so alike[END_REF]. Here, we aimed at answering the question to what extent there is a likely functional conservation between all orthologous genes expressed during leaf development in these distinct species.

With our analysis, we show that comparative analysis of orthologous gene expression profiles can help identifying potential new players involved in leaf growth. We found 217 one-toone orthology groups with conserved expression, of which the gene function is likely to be similar in both species. Maize and/or Arabidopsis candidates could therefore be selected for functional characterization. Proof of concept is demonstrated with LOW PSII ACCUMULATION 1 (LPA1), involved in photosystem II and chloroplast maintenance [START_REF] Peng | LOW PSII ACCUMULATION1 is involved in efficient assembly of photosystem II in Arabidopsis thaliana[END_REF]. It seems that the LPA1 expression is conserved as it takes part in a one-to-one orthology with GRMZM2G383154, previously described as a chloroplast maintenance gene as well [START_REF] Huang | Construction of plastid reference proteomes for maize and arabidopsis and evaluation of their orthologous relationships; the concept of orthoproteomics[END_REF]. Among these 217 orthology groups, 41 with unknown genes in both species were identified which could be potential leaf growth regulators (Table S1).

From the one-to-many orthology combined with the expression information, possible sub-or neofunctionalization of orthologous genes in both species could be shown by studying orthologous groups with divergent expression. This is best shown in one-to-two orthology, in which one of the two maize orthologues changed expression as a consequence of a duplication event, causing reduced selection to maintain the ancestral expression pattern and/or function, resulting in a tailored or new function for this orthologue. We were able to identify 59 one-totwo orthology groups, of which one of the two maize orthologues with divergent expression could have a new function (Table S4). Clearly, correctly sorting out these orthologues showing different levels of expression conservation is of key importance to select valid gene candidates when translating trait information from Arabidopsis to maize.

Furthermore, integration of transcriptomic data could help untangling many-to-many orthology and narrow down the genuine functional orthologues. For instance, REGULATORY PARTICLE AAA-ATPASE 2A (RPT2a), involved in the regulation of cell expansion duration [START_REF] Ueda | The HALTED ROOT gene encoding the 26S proteasome subunit RPT2a is essential for the maintenance of Arabidopsis meristems[END_REF], is part of a complex orthologous group with 59 Arabidopsis genes and 77 maize orthologues. By including DEGs during leaf development into this orthology network, the most likely functional maize orthologue of RPT2a, GRMZM2G056569, could be identified since it exhibits the same expression profile. Another example involves CYCLIN D3;1 (CYCD3;1) regulating the proliferation duration [START_REF] Dewitte | Arabidopsis CYCD3 D-type cyclins link cell proliferation and endocycles and are rate-limiting for cytokinin responses[END_REF]. CYCD3;1 is part of an orthologous group including ten Arabidopsis genes and 23 maize orthologues. After integration of the expression data, two maize orthologues, GRMZM2G140633 and GRMZM2G133413, with the same expression pattern, were found. These two genes could be more interesting candidates to study their role in cell proliferation in maize, compared with the other 21 maize orthologues (Supplemental network file). Narrowing down the orthologues to a set of potential functional orthologues by integrating gene expression data helps prioritizing and selecting potential functional orthologues simultaneously in Arabidopsis and maize for further functional studies.

To further aid the selection of potential new growth-regulating genes, we conducted a phenotypic analysis that revealed a stronger involvement of genes with orthologues in leaf-related phenotypes compared with specific-specific genes. Furthermore, this analysis yielded nine putative Arabidopsis candidates that are likely to increase leaf size when overexpressed (Table S3). Additionally, the phenotypic analysis allowed us to list four potential Arabidopsis-specific leaf-regulating genes that do not have a maize ortholog, but which do alter the leaf phenotype when mutated. Moreover, differences in leaf size and shape between both species could be partially shown as we identified eight Arabidopsis genes with nine maize orthologues exhibiting divergent expression involved in leaf growth.

Transcription factors and regulators that tightly control gene expression play an important role in plant development including leaf growth. Here, for 19 TF/TR families with DEGs, orthologues showed expression profiles that are similar during leaf growth in Arabidopsis and maize, suggesting that their transcriptional regulation is conserved. Some of the TFs have a well-known role in leaf development, such as the previously mentioned GRF, YABBY and E2F family members. Other genes having orthologues with conserved expression are involved in CR such as the Arabidopsis ATXR6 gene and its maize ortholog, PHD7, members of the SET regulators [START_REF] Baumbusch | The Arabidopsis thaliana genome contains at least 29 active genes encoding SET domain proteins that can be assigned to four evolutionarily conserved classes[END_REF]. ATXR6 interacts with the PROLIFERATING CELL NUCLEAR ANTIGEN (PCNA) and is up-regulated during the S phase of the cell cycle [START_REF] Raynaud | Two cell-cycle regulated SET-domain proteins interact with proliferating cell nuclear antigen (PCNA) in Arabidopsis[END_REF]. The expression profile of ATXR6 and PHD7, showing high expression during the cell proliferation phase, coincides therefore with that of the known SET genes.

In addition, we show evidence for expression conservation among putative orthologous targets of TFs, based on similar, enriched TFBSs in the promoters of these orthologues. This level of conservation is interesting, knowing the evolutionary distance between both species and suggests possible conservation of TFBSs in other plant species. For the E2F TFs, 12 Arabidopsis and 16 maize genes were identified as putative targets, of which all of the Arabidopsis targets were previously described as E2F targets. Even in more complex TF families such as the TCP family, conservation of the TFBS could be shown, albeit less evident. This low conservation might be because the TFBSs used in this study are Arabidopsis-centred and maize promoters might harbour adapted or even completely different TFBSs recruiting conserved TFs. Secondly, the size of the studied promoters (500 bp upstream from the start codon) was kept constant in both species, while maize promoters are considered larger and more diffuse as compared to Arabidopsis promoters since the majority of the maize genome is covered by transposable elements [START_REF] Lisch | How important are transposons for plant evolution?[END_REF]. In conclusion, through the cross-species analysis of gene expression profiles functional orthologous candidates with conserved expression could be identified for further functional characterization in both species during leaf development, bridging the knowledge between different model species for flowering plants.

Methods

Selection of differentially expressed genes and clustering

Normalized expression values from the transcriptome analysis of Arabidopsis [START_REF] Andriankaja | Exit from proliferation during leaf development in Arabidopsis thaliana: a not-so-gradual process[END_REF] and maize [START_REF] Nelissen | The reduction in maize leaf growth under mild drought affects the transition between cell division and cell expansion and cannot be restored by elevated gibberellic acid levels[END_REF] were used to select DEGs. The data sets for both species were processed in a similar manner. First, fold changes and false discovery rates (FDRs) for each gene were calculated by comparing each datapoint to every other datapoint. For the Arabidopsis data set, with six datapoints, 15 comparisons (6 9 (6 À 1)/2) were done per gene, and for the maize data set including nine datapoints, 36 comparisons were done. Next, DEGs were selected based on À2 ≥ fold change ≥2 and FDR ≤ 0.05. The resulting DEGs in both data sets corresponded to the union of all comparisons. Both DE data sets were clustered according to hierarchical clustering in R version 3.4.0 [START_REF] Core | R: A Language and Environment for Statistical Computing[END_REF]. Distances between genes were calculated with the Euclidean distance measure, and complete linkage was used as agglomeration method. All cluster plots were generated with the ggplot2 R-package [START_REF] Wickham | ggplot2: Elegant Graphics for Data Analysis[END_REF].

GO enrichment

GO enrichment analysis was done using the PLAZA 3.0 Dicots platform. GO enrichment results were filtered based on P value ≤ 1E À3 and fold change ≥2. Only biological process (and molecular in case of separate clusters) terms were considered. The heat maps that display the GO enrichment results were created in R version 3.4.4 with the ComplexHeatmap [START_REF] Gu | Complex heatmaps reveal patterns and correlations in multidimensional genomic data[END_REF]. The rows with GO terms were clustered with hierarchical clustering (Euclidean distance measure, clustering method complete linkage). Finally, the top ten GO terms with enrichment in at least one cluster or cluster overlap were displayed.

Orthologue identification between Arabidopsis and maize

PLAZA was used to determine the orthologous genes between Arabidopsis and maize [START_REF] Proost | PLAZA 3.0: an access point for plant comparative genomics[END_REF][START_REF] Van Bel | Dissecting plant genomes with the PLAZA comparative genomics platform[END_REF]. This platform includes nine species: four dicots [Arabidopsis thaliana (TAIR10), Arabidopsis lyrata (JGI v1.0), Populus trichocarpa (JGI assembly release v3.0, annotation v3.0) and Solanum lycopersicum (ITAG 2.4 on genome build SL2.50)], four monocot species [Zea mays (AGPv3.0), Sorghum bicolor (v3.1), Oryza sativa ssp. Japonica (MSU RGAP 7) and Brachypodium distachyon (v3.1)] and, additionally, an out-group species, Amborella trichopoda (Amborella v1.0), to root the phylogenetic tree of those species. This PLAZA version allows to identify orthologues with four methods, corresponding to orthologous gene families inferred through sequence-based clustering with OrthoMCL [START_REF] Li | OrthoMCL: identification of ortholog groups for eukaryotic genomes[END_REF], phylogenetic trees, multispecies best hits and inparalogs families [START_REF] Linard | OrthoInspector: comprehensive orthology analysis and visual exploration[END_REF] and synteny information, although we did not use the latter method. In order to establish an orthologous relation between two genes, we used a threshold of one prediction method.

Significant overlap between Arabidopsis and maize clusters

We used a random sampling analysis to calculate the significant overlap. Each Arabidopsis cluster was compared to each maize cluster. The overlap between two clusters was determined as the group of orthologous Arabidopsis and maize genes present in those clusters. In order to determine significant overlaps between Arabidopsis and maize clusters, the P-value was calculated by comparing randomly sampled maize clusters to each existing Arabidopsis cluster and counting the overlap again. This was done 1000 times so the P-value could be determined at the 0.001 level.

In silico leaf-related phenotypic analysis TAIR [START_REF] Berardini | The Arabidopsis information resource: making and mining the "gold standard" annotated reference plant genome[END_REF] and RARGE II [START_REF] Akiyama | RARGE II: an integrated phenotype database of Arabidopsis mutant traits using a controlled vocabulary[END_REF] were consulted in order to perform the phenotypic analysis. The TAIR germplasm list contained phenotypic information for 1841 Arabidopsis mutant lines. This gene list was filtered using the description of the phenotype based on the presence of the 'leaf' or 'leaves' keywords, and then manually curated for overexpressed or knocked out genes that lead to a leaf-related phenotype. This process resulted in a curated list of 291 genes. The RARGEII database included 17 809 Arabidopsis genes of which 2726 showed a phenotype. 453 genes were selected since they were annotated with rosette leaf terms according to the Plant Ontology structure: wrinkled, upturned, epinastic, coiled, increased or decreased in length, size and width, present in greater or fewer numbers in organism and abnormal. Genes involved in a pale, dark or variegated colouring of the rosette leaves were left out of the analysis. Combining curated data from both databases resulted in a list of 738 unique genes with leafrelated phenotypes when overexpressed or knocked out. For each gene group, the enrichment and significance of the leaf-related S8). Four different colours represent the four bases (yellow: G, red: T, blue: C and green: A), and the relative size of the bases corresponds to the frequency of occurrence. (c) Visual representation of the 34 Arabidopsis and 34 maize putative TCP3 and/or TCP24 targets that were identified based on the TFBS enrichment analysis. Genes (nodes) are coloured according to their expression trend: increasing profile (green) or decreasing profile (red) during leaf development, while the colour of the node edge shows the gene to be a target of TCP3 (blue), TCP24 (yellow) or both (purple). The thickness of the edge corresponds to the number of evidences found for that particular orthology link (B: best hits and inparalogs families, O: OrthoMCL and T: reconciled phylogenetic trees). The colour of the edges represents a significant orthology overlap between the respective Arabidopsis and maize genes. Ath., Arabidopsis thaliana; Zma., Zea mays; DE, differentially expressed. Functional orthology by comparative transcriptomics 563 phenotypes was calculated using the hypergeometric score and corrected for multiple testing with the Benjamini and Hochberg method [START_REF] Benjamini | Controlling the false discovery rate: a practical and powerful approach to multiple testing[END_REF].

Network visualization

Cytoscape 3.1 [START_REF] Shannon | Cytoscape: a software environment for integrated models of biomolecular interaction networks[END_REF] was used to visually represent the integrated comparative network. Orthologous information was used to create the network: genes were represented as nodes and their orthology relation as edges. Gene expression data together with species information and gene aliases were added as attribute format and visually represented by differently coloured nodes in the network. The .cys file is made available as supporting file.

Transcription factor binding site conservation

Cis-BP was consulted for all directly inferred TCP binding sites in position weight matrix (PWM) format of Arabidopsis [START_REF] Weirauch | Determination and inference of eukaryotic transcription factor sequence specificity[END_REF]. All TFBSs were mapped on the first 500 bp (including the 5 0 UTR) of the promoter using the Cluster-Buster algorithm with settings: background model based on 3-kb input sequences, cluster score threshold 0 and pseudo-count 0.375; [START_REF] Frith | Cluster-Buster: finding dense clusters of motifs in DNA sequences[END_REF]. All mappings were ranked according to P-value, and the top 3000 hits for unique genes were used for calculating the enrichment. The enrichment of all TFBSs in each Arabidopsis and maize cluster was calculated employing the hypergeometric score and corrected for multiple testing with the Benjamini and Hochberg method [START_REF] Benjamini | Controlling the false discovery rate: a practical and powerful approach to multiple testing[END_REF]. S1 Overview of the 217 one-to-one orthology groups with conserved expression and part of one of the significant overlaps (A1-Z2, A5-Z5, A5-Z6, A6-Z5, A6-Z6) Table S2 Calculation of percentage, enrichment, and FDR of leafrelated phenotype data for the groups with conserved or divergent expression, and the Arabidopsis group with no orthologues of clusters A1, A5, and A6 Table S3 Group, ID, name, TF information, and phenotype description resulting from a knockout (KO) or an overexpression (OE) described in TAIR and or RARGE and confirmed by a publication for 49 genes, 10 genes, and 48 genes possibly involved in a decrease or increase in leaf size or another leafrelated phenotype, respectively Table S4 Overview of the 59 one-to-two orthology groups exemplifying possible neo -or subfunctionalization of the maize orthologue having divergent expression Table S5 Overview of 211 DE Arabidopsis TFs and their 256 DE maize orthologues either with conserved or divergent expression throughout leaf development Table S6 Overview of the DE TF families including Arabidopsis genes with their maize orthologues that are conserved, partially conserved, or divergent in expression throughout leaf development Table S7 Overview of the identified E2F targets with conserved expression in Arabidopsis (12 genes) and maize (16 genes) Table S8 Overview of the Arabidopsis and maize clusters in which TCP3 and TCP24 TF binding sites are enriched based on the hypergeometric score (HGS) Supplemental network file with the full integrated network as well as the case-studies handled in this manuscript.

Figure 1

 1 Figure 1 Overview of the properties of the orthologous network constructed between Arabidopsis and maize using the integrative orthology method in PLAZA. (a) Visual representation of Arabidopsis GROWTH-REGULATING FACTOR 3 (GRF3) (diamond) and its orthologues in maize. Three maize orthologues (upper three triangles) are confirmed by multispecies Best-Hits-and-Inparalogs families (B) and reconciled phylogenetic trees (T), and twelve additional orthologues (lower 12 triangles) are predicted by T. Arabidopsis genes are represented by diamonds and maize genes by triangles. (b) and (c) Venn diagrams depicting the number and the percentage of Arabidopsis (b) and maize (c) genes sharing orthologues with the other species through the three evidence sources in PLAZA. (d) Pie chart showing the percentage of orthologous groups with one-to-one, one Arabidopsis gene with many maize orthologues (one-to-many Arabidopsis), one maize gene with many Arabidopsis orthologues (one-to-many maize) or many-to-many orthologous links.

Figure 2

 2 Figure 2 Clustered expression profiles of DEGs in Arabidopsis and maize during leaf development together with the orthology overlap between clusters.The six selected Arabidopsis clusters are depicted horizontally (A1-A6, diamonds) and the seven selected maize clusters are ordered vertically (Z1-Z7, triangles), all of them sorted according to their visual trend. Clusters A1, A2 and Z1 have high expression in the early stages, corresponding to the proliferation phase, and a decreasing expression as leaf development advances (red). Clusters A3, A4, Z2, Z3 and Z4 show high expression during transition and low expression at the beginning and later stages of leaf development (blue). Clusters A5, A6, Z5 and Z6 have an upward gene expression trend coinciding with the onset of the expansion phase (green). Cluster Z7 shows high expression during proliferation and expansion, and low expression during the transition phase (yellow). The total number of genes in each cluster and the number of DEGs with an DE orthologue are shown under or right of the cluster name for Arabidopsis or maize, respectively. For each comparison between two clusters, the overlap is shown for both the Arabidopsis (upper number) and maize genes (lower number). Eight significant overlaps (P < 0.001) are shown in dark grey, and two significant overlaps (P < 0.01) are shown in light grey. At the bottom right, two representative clusters of the Arabidopsis (left) and the maize (right) DEG set are shown. Normalized gene expression is depicted on the x-axis according to the different datapoints on the y-axis, which corresponds to 8-13 DAS for the Arabidopsis data set and to distance to the leaf base for the maize data set. Ath., Arabidopsis thaliana; Zma., Zea mays; Nr., number; sp., species; DE, differentially expressed; DAS, days after stratification.

Figure 3

 3 Figure3Heat map showing the top 10 enriched GO terms (biological process) per overlap with enrichment in at least one overlap for the four significant overlaps (P < 0.001) with conserved expression showing a peak in expression during cell expansion (green) and for the one significant overlap (P < 0.001) with conserved expression with an expression peak during cell proliferation (red).

Figure 4

 4 Figure 4 Overview of the DEGs from four TF families, YABBY, GRF, E2F and NAC, showing conserved expression patterns in both Arabidopsis and maize. In each panel, the top image represents the integrated orthology network in which Arabidopsis TFs (diamonds) and maize orthologues (triangles) are linked. The colour of the nodes corresponds to an increasing (green) or a decreasing (red) expression trend during leaf development. The thickness of the edge corresponds to the number of evidences found for that particular orthology link (B: best hits and inparalogs families, O: OrthoMCL and T: reconciled phylogenetic trees). The colour of the edges represents a significant orthology overlap between the respective Arabidopsis and maize genes. The heat map in each panel represents the normalized expression values of each DE TF throughout leaf development. Ath., Arabidopsis thaliana; zma., Zea mays;DAS, days after stratification.

ª

  2019 The Authors. Plant Biotechnology Journal published by Society for Experimental Biology and The Association of Applied Biologists and John Wiley & Sons Ltd., 18, 553-567

Figure 5

 5 Figure 5 Overview of the conserved expression among TCP TFs, their orthologues and their putative targets. (a) Overview of the Arabidopsis (diamonds) class I and class II TCP DEGs with their DE maize (triangles) orthologues during leaf growth. The heat map represents the normalized expression values of each DE TF throughout leaf development in both species. (b) Visual representation of eight experimentally verified TF binding sites (TFBSs) corresponding to eight TCP TFs, three class I TFs and five class II TFs. The TFBSs of the underlined TFs, TCP3 and TCP24, were found enriched in promoters of Arabidopsis and maize gene clusters (TableS8). Four different colours represent the four bases (yellow: G, red: T, blue: C and green: A), and the relative size of the bases corresponds to the frequency of occurrence. (c) Visual representation of the 34 Arabidopsis and 34 maize putative TCP3 and/or TCP24 targets that were identified based on the TFBS enrichment analysis. Genes (nodes) are coloured according to their expression trend: increasing profile (green) or decreasing profile (red) during leaf development, while the colour of the node edge shows the gene to be a target of TCP3 (blue), TCP24 (yellow) or both (purple). The thickness of the edge corresponds to the number of evidences found for that particular orthology link (B: best hits and inparalogs families, O: OrthoMCL and T: reconciled phylogenetic trees). The colour of the edges represents a significant orthology overlap between the respective Arabidopsis and maize genes. Ath., Arabidopsis thaliana; Zma., Zea mays; DE, differentially expressed.

ª

  2019 The Authors. Plant Biotechnology Journal published by Society for Experimental Biology and The Association of Applied Biologists and John Wiley & Sons Ltd., 18, 553-567

Figure S1

 S1 Figure S1 Visual representation of the two data sets and overview of the hierarchical clustering of the Arabidopsis and maize DEGs.Figure S2 Heat map showing the top 10 enriched GO terms (BP and MF) per cluster with enrichment in at least one cluster.Figure S3 Visual representation of the Arabidopsis class I (left) and class II (right) TCP gene family with the maize orthologues.Figure S4 Schematic overview of the identified putative E2F targets benchmarking. TableS1Overview of the 217 one-to-one orthology groups with conserved expression and part of one of the significant overlaps (A1-Z2, A5-Z5, A5-Z6, A6-Z5, A6-Z6) TableS2Calculation of percentage, enrichment, and FDR of leafrelated phenotype data for the groups with conserved or divergent expression, and the Arabidopsis group with no orthologues of clusters A1, A5, and A6 TableS3Group, ID, name, TF information, and phenotype description resulting from a knockout (KO) or an overexpression (OE) described in TAIR and or RARGE and confirmed by a publication for 49 genes, 10 genes, and 48 genes possibly

Figure S2

 S2 Figure S1 Visual representation of the two data sets and overview of the hierarchical clustering of the Arabidopsis and maize DEGs.Figure S2 Heat map showing the top 10 enriched GO terms (BP and MF) per cluster with enrichment in at least one cluster.Figure S3 Visual representation of the Arabidopsis class I (left) and class II (right) TCP gene family with the maize orthologues.Figure S4 Schematic overview of the identified putative E2F targets benchmarking. TableS1Overview of the 217 one-to-one orthology groups with conserved expression and part of one of the significant overlaps (A1-Z2, A5-Z5, A5-Z6, A6-Z5, A6-Z6) TableS2Calculation of percentage, enrichment, and FDR of leafrelated phenotype data for the groups with conserved or divergent expression, and the Arabidopsis group with no orthologues of clusters A1, A5, and A6 TableS3Group, ID, name, TF information, and phenotype description resulting from a knockout (KO) or an overexpression (OE) described in TAIR and or RARGE and confirmed by a publication for 49 genes, 10 genes, and 48 genes possibly

Figure S3

 S3 Figure S1 Visual representation of the two data sets and overview of the hierarchical clustering of the Arabidopsis and maize DEGs.Figure S2 Heat map showing the top 10 enriched GO terms (BP and MF) per cluster with enrichment in at least one cluster.Figure S3 Visual representation of the Arabidopsis class I (left) and class II (right) TCP gene family with the maize orthologues.Figure S4 Schematic overview of the identified putative E2F targets benchmarking. TableS1Overview of the 217 one-to-one orthology groups with conserved expression and part of one of the significant overlaps (A1-Z2, A5-Z5, A5-Z6, A6-Z5, A6-Z6) TableS2Calculation of percentage, enrichment, and FDR of leafrelated phenotype data for the groups with conserved or divergent expression, and the Arabidopsis group with no orthologues of clusters A1, A5, and A6 TableS3Group, ID, name, TF information, and phenotype description resulting from a knockout (KO) or an overexpression (OE) described in TAIR and or RARGE and confirmed by a publication for 49 genes, 10 genes, and 48 genes possibly

Figure S4

 S4 Figure S1 Visual representation of the two data sets and overview of the hierarchical clustering of the Arabidopsis and maize DEGs.Figure S2 Heat map showing the top 10 enriched GO terms (BP and MF) per cluster with enrichment in at least one cluster.Figure S3 Visual representation of the Arabidopsis class I (left) and class II (right) TCP gene family with the maize orthologues.Figure S4 Schematic overview of the identified putative E2F targets benchmarking. TableS1Overview of the 217 one-to-one orthology groups with conserved expression and part of one of the significant overlaps (A1-Z2, A5-Z5, A5-Z6, A6-Z5, A6-Z6) TableS2Calculation of percentage, enrichment, and FDR of leafrelated phenotype data for the groups with conserved or divergent expression, and the Arabidopsis group with no orthologues of clusters A1, A5, and A6 TableS3Group, ID, name, TF information, and phenotype description resulting from a knockout (KO) or an overexpression (OE) described in TAIR and or RARGE and confirmed by a publication for 49 genes, 10 genes, and 48 genes possibly

ª

  2019 The Authors. Plant Biotechnology Journal published by Society for Experimental Biology and The Association of Applied Biologists and John Wiley & Sons Ltd., 18, 553-567

ª

  2019 The Authors. Plant Biotechnology Journal published by Society for Experimental Biology and The Association of Applied Biologists and John Wiley & Sons Ltd., 18, 553-567

  , were enriched in clusters A5 and Z5 with high expression during expansion and ª 2019 The Authors. Plant Biotechnology Journal published by Society for Experimental Biology and The Association of Applied Biologists and John Wiley & Sons Ltd., 18, 553-567

ª 2019 The Authors. Plant Biotechnology Journal published by Society for Experimental Biology and The Association of Applied Biologists and John Wiley & Sons Ltd., 18, 553-567

ª 2019 The Authors. Plant Biotechnology Journal published by Society for Experimental Biology and The Association of Applied Biologists and John Wiley & Sons Ltd., 18, 553-567 Jasmien Vercruysse et al.

Acknowledgements

We thank Annick Bleys for proofreading and preparing the manuscript. This work was supported by the European Research Council under the European Community's Seventh Framework Programme [FP7/2007-2013] under ERC grant agreement no [339341-AMAIZE]11 and by Ghent University ('Bijzonder Onderzoeksfonds Methusalem project' no. BOF08/01M00408 and Multidisciplinary Research Partnership 'Biotechnology for a Sustainable Economy' Grant 01MR0510W).

Author contributions

J.V., K.V., N.G. and D.I. designed the research methodology, performed data analysis and wrote the article. M.V.B. assisted in computing orthologues. S.R.K. helped with analysis of the PWMs and C.O.-C. generated the heat maps. V.S. preprocessed the transcriptomic data sets. H.N. assisted in the biological interpretation of the maize genes.

Conflict of interest

The authors declare no conflict of interest.

Supporting information

Additional supporting information may be found online in the Supporting Information section at the end of the article.