Architecture of gene regulatory networks controlling flower development in Arabidopsis thaliana

Abstract

Floral homeotic transcription factors (TFs) act in a combinatorial manner to specify the organ identities in the flower. However, the architecture and the function of the gene regulatory network (GRN) controlling floral organ specification is still poorly understood. In particular, the interconnections of homeotic TFs, microRNAs (miRNAs) and other factors controlling organ initiation and growth have not been studied systematically so far. Here, using a combination of genome-wide TF binding, mRNA and miRNA expression data, we reconstruct the dynamic GRN controlling floral meristem development and organ differentiation. We identify prevalent feed-forward loops (FFLs) mediated by floral homeotic TFs and miRNAs that regulate common targets. Experimental validation of a coherent FFL shows that petal size is controlled by the SEPALLATA3-regulated miR319/TCP4 module. We further show that combinatorial DNA-binding of homeotic factors and selected other TFs is predictive of organ-specific patterns of gene expression. Our results provide a valuable resource for studying molecular regulatory processes underlying floral organ specification in plants.

Fig. 1

A DNA-binding atlas of key-regulatory transcription factors in flower development. a Schema depicting the key developmental stages for flower formation, including floral transition, floral initiation (stage 0, S0), meristem specification (S2), organ specification (S4), organ differentiation (S8), and flower maturation. The width of each color bar roughly defines the time range in which the plant materials were harvested for a TF ChIP-seq experiment. Note that colors representing stages are consistently used across all the figures if applicable. b Currently available genome-wide TF-binding maps for floral regulators. The datasets are roughly assigned to the five floral stages as defined in a; the source of ChIP-seq data is labeled below the bar and also indicated in Supplementary Table 1. c Bar chart showing the number of binding sites for the investigated TFs. TFs are colored according to the development stages. d Bar chart showing the percentage of target genes as defined in three classes: transcription factors (TF; blue), miRNA-related genes (MIRs; orange) and other protein-coding genes (gray). The number of target genes is shown above bars. e Cross-regulatory interactions between 15 TFs in regulatory networks across five stages. TFs are colored and arranged in the same order along each axis. Stage-specific edges are colored as indicated, whereas regulatory interactions present in more than one stage networks are shown in gray. f Co-binding relationships and co-regulated targets by TFs. Upper triangle, highly significant co-regulation relationships are assigned a dark blue color, representing −log₁₀(P-value). The P-value was calculated from a hypergeometric test to check the significance of target overlapping between any two TFs using all the annotated genes as the background. Lower triangle, highly significant co-binding relationships are assigned a dark orange color. The significance of co-binding by any two TFs were tested by Jaccard statistics, which measures the ratio of the number of intersecting base pairs occupied by both TFs to the number of base pairs in their union. Dashed boxes indicate highly interplayed regulators in vegetative development (green) and flower development (purple). g Venn diagram of the target genes shared by AG, AP1, SEP3, AP3, and PI at stage 4

Fig. 2

Expression dynamics and functions of target genes regulated by floral regulators. a Venn diagram showing significant overlap (P-value calculated by a hypergeometric test) of target genes by floral transcription factors (TFs; green) and differentially expressed (DE) genes across flower development (red). b Clustering analysis of DE target genes. Expression values are center normalized per gene based on its log₁₀-transformed FPKM values in the four representative stages. The differentially expressed regulator genes are labeled. S0: stage 0. c Boxplot showing the distribution of expression changes (the maximum value of absolute fold change over stages) with increasing number of bound TFs. The number of genes in each category is labeled below the box. Boxes represent quartiles, center lines denote 50th percentile, and whiskers extend to most extreme values within 1.5× interquartile range (IQR). Significance codes, ****: P-value <0.0001, ***: P-value <0.001, **: P-value <0.01, *: P-value <0.05, and NS: not significant, by two-sided Mann–Whitney tests. d Expression and functional characteristics of DE genes grouped by the bound TFs. Up- and downregulated DE genes are indicated in b. BP: biological process. e Distribution of bound regulators per potential target gene. Curves illustrate the ensemble of 100 randomized networks of the same size and degree distribution. Examples of hub target genes (bound by seven regulators or more) are labeled

Fig. 3

Characterizing flowering gene regulatory networks and motifs. a Hive plot showing gene regulatory networks involved in floral regulators, miRNA genes, and transcription factors (TFs). Connections indicate potential regulation relationships. Regulator-miRNAs and regulator-TFs relationships are supported by ChIP-seq data, while miRNA-TFs relationships are based on prediction. For valorization purpose, nodes are colored according to TF or miRNA gene families or regulators, and edges are colored according to the corresponding source or target. Note that only differentially expressed TF target genes are shown. b Robustness of miRNA-mediated gene regulator networks. Plot shows the connectivity loss due to iteratively removing nodes with the maximum load (betweenness) in the networks with or without (w/o) miRNA-mediated connections. P-value was calculated based on a two-sided Kolmogorov–Smirnov test. c Systematic search of one-node, two-node, and three-node motifs from the network represented in a. The number of motif occurrences is shown below the schema. For three-node motifs, only statistically enriched (significance code denoted in parenthesis) motifs are shown. Examples of corresponding motifs are listed at the bottom panels d–g. d Autoregulation as the simplest motif consisting of the master regulators themselves. e The feedback regulation between MIR172 and AP2. f Analysis of the single-input module motifs (SIMs). A floral regulator may target two TFs that are connected via protein–protein interactions (PPIs; left schema). The right boxplot shows the percentage of target genes with PPIs. As control, the same number of target genes were sampled from the whole target gene list and the percentage of targets with PPIs were calculated. P-value was calculated based on a two-tailed Student’s t-test. Boxes represent quartiles, center lines denote 50th percentile, and whiskers extend to most extreme values within 1.5× interquartile range (IQR). g A feed-forward loop (FFL) motif mediated by MIR319a

Fig. 4

SEP3 as a direct regulator of MIR319a. a Dynamics of SEP3 binding at the loci of MIR319a (top), TCP4 (middle), and TCP10 (bottom). Arrows indicate transcriptional directions. b Pictures of single flowers and separated petals as well as other flower organs (top panel; scale bar = 1 mm) and inflorescences (bottom) of wild-type (WT) plants, enhanced expression of SEP3 (SEP3^HE4.1k), jaw-D mutant lines, and combination (F2) of jaw-D mutant and SEP3^HE4.1k. c Boxplots showing the petal width (left) and petal length (right) of different genotypes in b. The number of biologically independent plants is shown in parentheses. Boxes represent quartiles, center lines denote 50th percentile, and whiskers extend to most extreme values within 1.5× interquartile range (IQR). Significance codes, ****: P-value <0.0001, *: P-value <0.05, and NS: not significant, by two-sided Mann–Whitney tests. Bar = mean ± s.e. d Relative expression of MIR319a, TCP4, and TCP10 in the inflorescences of WT, jaw-D mutant, SEP3^HE4.1k lines, and three independent F2 lines produced by crossing jaw-D and SEP3^HE4.1k lines. e Relative expression of TCP4 and MIR319a in either whole inflorescences or stage 10-enriched flower tissues of WT and SEP3 enhanced expression lines. For all the expression analysis, data represent mean of three independent biological replicates. The Tip41 gene (AT4G34270) was used as reference. Significance codes as in c. f SEP3-MIR319a-TCP4/TCP10 coherent feed-forward loop in regulation of petal development

Fig. 5

Organ-specific gene expression controlled by floral master regulators. a Clustering analysis of domain-specific gene expression. Five major clusters correspond to five distinct spatiotemporal domains (colored in left bars): AP1-, AP3-, AG-specific, “AP1/AP3” (AP1 and AP3-common domains), and “AP3/AG” domains. Dendrograms above the heatmap illustrate the overall similarity of gene expression patterns. Example genes are indicated in the right. Domain expression data at stages 4 (S4) and 6 (S6) were obtained from ref. . b 3D t-SNE (t-distributed stochastic neighbor embedding) plot showing the overall similarity of transcription factor (TF)-binding profiles of organ-specific genes. Normalized ChIP-seq intensity for the 11 selected TFs were used in the analysis. c Predicting the contributions of promoter binding by multiple TFs to gene expression difference between two different domains by a Lasso regression model. Regression coefficients are plotted in heatmap as relative importance of the binding features. d Domain-specific genes targeted by floral quartets that determine different floral organ identities, i.e., sepals, petals, stamens, and carpels, as colored in the diagram of floral organs. The number of organ-specific genes is shown in parentheses. e Network showing domain expression patterns of target genes for the four floral homeotic tetrameric MADS TF complexes as indicated in d. Nodes as pie charts show gene expression patterns in AP1, AP3, and AG domains at stage 4. Edges represent target genes by the corresponding floral quartet complexes