bri1–5/bak1–1D, bri1–5/brs1–1D and bri1–5/bri1–1D partially reconstitute bri1–5 gene expression
To better understand the molecular mechanisms of key BR signaling genes, we performed a phenotypic screening and expression analysis of bri1–5 and its three activation-tag suppressors along with their corresponding wild-type, WS2. Two suppressor strains bri1–5/bak1–1D and bri1–5/brs-1D were obtained from [12]. An additional bri1–5/bri1–1D mutant was generated in the framework of the current study (see Methods). Sequencing the BRI1 flanking region from the suppressor bri1–5/bri1–1D showed that the activation tag was inserted 534 bp downstream of the BRI1 gene (Supplementary Fig. 1-A). All suppressor mutants were shown to indeed overexpress the activation tagged gene as confirmed by Real-Time qPCR (RT-qPCR) (Table S1). Phenotypically, all bri1–5 suppressors (bri1–5/bak1–1D, bri1–5/brs1–1D, and bri1–5/bri1–1D) lines displayed larger seedlings than the bri1–5 mutant, but still significantly smaller than the WS2 (Fig. 2). Of all suppressor mutants, the bri1–5/bak1–1D line best approximated the growth phenotype of the WS2, and its larger seedling seemed to be mainly the effect of its larger root length and to a lower extent of its larger hypocotyl length (both of which were significantly larger than the bri1–5 mutant). The contribution of the epidermal cell length in recovering the bri1–5 is marginal in the bri1–5/bak1–1D (line with the largest seeding) but seems much more pronounced in the bri1–5/brs1–1D (Fig. 2, Supplementary Fig. 1: B-F). This indicates that in the bri1–5/brs1–1D mechanisms other than those in the bri1–5/bak1–1D line play a role in alleviating the bri1–5 phenotype.
To gain insight into which pathways in each of the studied lines were responsible for recovering the bri1–5 growth phenotype to wild-type level, we performed gene expression analysis. All suppressor lines, together with the wild-type (WS2) and the bri1–5 background were sampled at a 7-day seedling stage. To assess the reproducibility of the expression analysis, we measured the extent to which the expression profiles of replicate samples were similar using Principal Component Analysis (PCA): PCA indeed showed that the largest fraction of the variation in gene expression between the samples could be assigned to differences in genetic background and not to differences between replicates of the same genetic background, confirming the reproducibility (Fig. S2). In addition, microarray results were confirmed using RT-qPCR for a randomly selected set of differentially expressed genes (Fig. S3).
We determined for each mutant line its differential expression versus the same common reference i.e. the expression state in WS2, resulting in a total of 1413 differentially expressed genes (Additional file 1). The Venn diagram represented in Fig. 3 shows to what extent the different lines share the same differentially expressed genes (aberrantly expressed versus the WS2 control). Fig. 3 and the scatter-plots in Fig. S4 (A-C) show that of all suppressor lines, bri1–5/bak1–1D could restore the largest number of genes that were affected in expression in bri1–5 (about two-thirds of the genes that were differentially expressed in bri1–5 were no longer differentially expressed in bri1–5/bak1–1D). This is in line with its observed phenotypic behavior as indeed bri1–5/bak1–1D seems to also phenotypically best compensate for the bri1–5 mutation.
The Venn diagram in Fig. 3 also shows that the bri1–5/brs1–1D and bri1–5/bri1–1D lines share the largest fraction of similarly affected genes. The latter is also illustrated in Fig. S4 panel D-F which shows that from all pairwise comparisons between suppressor lines, the level of differential expression relative to the mutant bri1–5 is most correlated between the suppressor lines bri1–5/brs1–1D and bri1–5/bri1–1D (i.e. R2 = 0.20). This suggests a similar role for BRI1 and BRS1 in the BR signaling pathway. Note that in Fig. S4 D-F, rather than performing a direct correlation analysis of the expression between two mutant lines, we performed correlation analysis with the expression of each mutant line relative to the same reference (expression in bri1–5). In this way, the correlation analysis is driven by the expression of the genes that change their expression relative to bri1–5. Although this results in lower correlation values than when directly comparing the expression values of the mutant lines, it better reflects the consistency between mutant lines in restoring genes affected in the bri1–5 mutant.
To confirm the extent to which the different suppressor strains molecularly restore the defects in the bri1–5 mutation, we compiled a list of marker genes representative of downstream pathways affected by BR signaling (Additional file 2). This consisted of 233 marker genes that were according to literature regulated by BR signaling (genes that became up or down-regulated upon treatment with exogenous BRs or by overexpressing the BR signaling genes). Of those marker genes, only those that were significantly affected in the bri1–5 line were retained in order to identify the mutant line that best suppresses the bri1–5 mutation (96 marker genes). Fig. 4 and Fig. S5 show how the expression of these genes is, as compared to WS2 affected in the bri1–5 mutant and how some of those genes got restored in the suppressor mutants. These results confirm what we observed based on the global expression analysis, i.e. that the bri1–5/bak1–1D restored the bri1–5 affected marker genes to the largest extent, and that molecularly the bri1–5/brs1–1D and bri1–5/bri1–1D mutant tend to behave more similarly in restoring the same marker genes.
Identifying compensatory and restoring pathways
Pathway analysis (see Methods) unveiled the pathways overrepresented amongst the differentially expressed gene sets in each of the mutant lines. Fig. S6, S7 and S8 and Table S2 show a number of pathways that are differentially expressed in both bri1–5 and all of the suppressor lines. These represent the pathways that are responsible for the aberrant growth phenotype in the bri1–5 mutant and that could not entirely be restored or compensated for in the suppressor lines. Among others, pathways related to cell wall synthesis (cell wall cellulose synthesis), protein and lipid metabolism can explain the residual discrepancy between the WS2 growth phenotype and the suppressors.
We assumed that if the suppressor strains alleviate the phenotype of the bri1–5 mutant, they could do so because they either restore the pathways disrupted in the bri1–5 mutant to wild-type levels or they induce genes that compensate for the bri1–5 affected pathways. Both mechanisms are reflected in the expression data. Processes that are aberrantly expressed in the bri1–5 mutant, but not in any of the suppressor lines, represent pathways that are restored to WS2 levels in all of the suppressors. This seems to be the case for some genes related to cytochrome P450 oxidase (Table S2). The fact that they are restored (or not significantly affected) in any of the suppressors indicates they might be essential for the recovery of the WS2 phenotype. Interestingly, the genes related to “glutathione S transferases” (Fig. S7, Table S2) are largely down-regulated in the bri1–5 mutant, restored to normal in bri1–5/bri1–1D and bri1–5/bak1–1D and up-regulated as compared to WS2 levels in the bri1–5/brs1–1D suppressor, indicating that some overcompensation is needed for this pathway in the bri1–5/brs1–1D background in order to restore the bri1–5 phenotype. In addition, ABA-related metabolism (Fig. S8 and Table S2) seems to have been affected by all suppressors, but at least not to a significant level in the bri1–5 mutant. Therefore, ABA signaling seems to represent a compensatory pathway, i.e. a pathway that needs to be triggered in the suppressor strains in order to restore the bri1–5 affected pathways and phenotype.
As less than 5000 genes can be mapped using pathway analysis, we performed a more elaborate analysis using a network-based approach. Network analysis provides an intuitive way of combining expression data with prior information on known molecular interactions or already available functional data [23, 24]. This approach first maps candidate genes, that are identified through expression analysis, on an integrated molecular interaction network. Then, it identifies subnetworks that connect as many candidate genes as possible [24]. By leveraging candidate genes identified through expression analysis with known interaction information, spuriously identified candidate genes can be removed as they will not be part of the subnetworks. In addition, genes relevant to the process of interest that are themselves not regulated at the level of expression are indirectly identified by being part of a connected component/subnetwork to which also many of the candidate genes belong. Such an integrated analysis provides a more comprehensive view of the process of interest. Here, we applied such an integrated network-based strategy to gain a more in-depth insight into the molecular mechanisms through which the suppressor lines can restore the bri1–5 phenotype to WS2 levels. To perform this network analysis, we started from the gene sets depicted in Fig. 3.
Involvement of hormone signaling in alleviating the bri1–5 phenotype
To study the interaction between non-restored, restored, and compensatory pathways in more depth, we combined the following gene sets for network analysis (see Fig. 3): i) genes that were most likely restored in the suppressors (genes of group A i.e. the genes with altered expression in the bri1–5 mutant, but restored to WS2 level in at least 2 suppressors), ii) genes that were compensatory in most of the suppressors (genes of group B i.e. the genes not differentially expressed in the bri1–5 mutant, but differentially expressed in at least two suppressor strains) and iii) genes altered in the bri1–5 mutant that most likely were not restored in the suppressors (genes of group C i.e. the genes, differentially expressed in the bri1–5 mutant and at least two of the suppressor strains). This combined set of genes (789 genes) is referred to as the set of seed genes or the genes we want to maximally connect on the interaction network.
Network analysis (see Material and Methods) identified 8 sub-networks (Fig. 5) containing the set of seed genes that could be connected through the interaction network. These subnetworks contain not only seed genes, but also connector genes. These are genes that are not differentially expressed themselves, but that are still recovered by the network analysis, because of their high connectivity with seed genes. As they are needed to connect seed genes in the network, they are most likely involved in the same processes as the seed genes. The subnetworks were annotated based on their enrichment in known GO functions (being enriched in respectively negative regulation of ABA, response to auxin, fatty acid metabolism process, developmental process, oligopeptide transport, response to ROS, BR homeostasis, and Ethylene activated signaling (Fig. 5)). This indicates that these are the pathways that contribute to alleviating bri1–5 signaling deficiency in the suppressor strains.
In-depth analysis shows that the subnetwork enriched in ABA signaling (Fig. 5, subnetwork 1) contains several known negative regulators of ABA signaling: HAI1, HAB1 ABI1, ABI2, and PP2CA acted as compensatory genes: these were are up-regulated in at least two bri1–5 suppressor lines compared to wild-type, but were not affected in the bri1–5 mutant (HAI1 and HAB1 being significantly up-regulated in all suppressor mutants; ABI1, ABI2, PP2CA, being significantly up-regulated in two suppressors, see Fig. S9); In addition, HAI2 was affected in the bri1–5 line, but could not be restored in at least two suppressors (non-restored gene), and HAB2 was identified as a connector node.
Interestingly, several targets of the ABA signaling pathway (DTX50, HVA22D, PUB19, COR15B, next to HAI1, HAB1) were identified as differentially expressed in all three suppressors (identified based on a GO enrichment of the core of group B, but not in bri1–5). This indicates that ABA signaling has indeed been affected in the suppressor strains to compensate for the bri1–5 signaling deficiency. Of these, DTX50, HVA22D, PUB19, COR15B could not be connected by PheNetic on the interaction network, implying they are either not annotated in the interaction network (COR15B) or quite distantly located from each other in the network.
The aforementioned negative regulators of ABA signaling in subnetwork 1 belong to the protein phosphatase 2C (PP2C) gene family which has nine members in total (HAI2, HAB2, HAB1, HAI3, PP2CA, ABI1, AHG1, ABI2, AHI1). PP2C is known to indirectly repress ABI5, the main activator of ABA signaling [25]. PP2C is also known to repress BIN2 activity [4, 5]: as BIN2 activates ABI5 by phosphorylating SnRKs [4, 5], repressing ABA signaling by PP2C via blocking SnRKs phosphorylation seems to compensate for the deficiency in BRI1 mediated signaling (Fig. 1). The subnetwork enriched in ABA signaling (network 1) also contains members of the PYR/PYL/RCAR family as connector genes (RCAR5, RCAR6, RCAR7, RCAR10, RCAR14, PYL8). The PYR/PYL/RCAR family constitutes the receptor of ABA signaling and promotes the activation of SnRKs by repressing PP2C [15, 26]. The fact that the SnRKs (SnRK2.5) and PYR/PYL/RCAR genes were identified as connector genes implies that they are likely involved in the pathways that connect the affected, restored, and compensatory genes of subnetwork 1. They are most likely not primarily regulated at the expression level, given their role in phosphorylation-mediated signaling [5, 27]. This explains why they were detected as connector genes and not retrieved by differential expression analysis.
We could not find any link in the literature to explain how PP2C can be up-regulated by BR signaling in order to repress ABA signaling. It seems that there exist some missing links between BZR1/BES1 (or downstream TFs) and the PP2C gene family. By partially recovering BR signaling in the suppressor lines, one would expect that PP2C gene expression levels restore to WS2 expression level. However, they appear to become up-regulated in suppressors, indicating that further compensatory repression of the ABA signaling is required in order to restore the bri1–5 phenotype.
Other than the ABA subnetwork, subnetworks related to other hormone signaling processes like auxin signaling (subnetwork 2), ROS signaling (subnetwork 6), and ethylene signaling (subnetwork 8) were also detected. It is well known that crosstalk between these phytohormone signaling pathways exists [28, 29]. Hence, interfering with one pathway e.g. ABA signaling through BR signaling might affect other phytohormone pathways as well. Based on the pathway analysis using MapMan and the network analysis we can conclude that ABA signaling is mostly a compensatory pathway (genes indicated in red color in Fig. 5 represent compensatory genes), whereas auxin, ethylene, and ROS signaling pathways are at least partially restored to the WS2 level (genes indicated in green color in Fig. 5 represent restored genes). However, restoring those pathways to the WS2 level seems to also depend on the presence of at least some compensatory genes (genes indicated in red color in Fig. 5 represent compensatory genes). The remaining subnetworks are enriched for fatty acid metabolism (subnetwork 3), developmental processes (subnetwork 4), and oligopeptide transport (subnetwork 5) confirming that all suppressors could partially restore some affected primary metabolic pathways. This is in line with the MapMan results and the recovered phenotypes.
Negative feedback of BR signaling on BR biosynthesis
The network result shows that the BR biosynthesis (subnetwork 7) is not well recovered in any of the suppressors. This subnetwork contains 4 genes (CYP90C1, CYP90B1/DWF4, CYP85A2, CYP90A1/DWF3) that are differentially expressed in the bri1–5 lines and all of the suppressors. Those genes, belonging to the cytochrome P450 superfamily play a role in BR biosynthesis by converting the sterol “campesterol” to BRs [30]. CYP708A3, another BR biosynthesis gene was found to be differentially expressed in all suppressors and the bri1–5 mutant (Fig. 6). However, as this gene was not present in the interaction network, it was missed by the network analysis. Unlike the four aforementioned cytochrome P450 genes, the molecular function of the CYP708A3 gene is still unknown. Interestingly it shows an expression pattern that is anticorrelated to that of the other cytochrome P450 genes (Fig. 6). This result along with the fact that CYP708A3 is known to be up-regulated by exogenous brassinolide (BL) treatment [31, 32] suggests it acts as an inhibitor of BRs biosynthesis.
The fact that the expression of BR biosynthetic genes is affected by mutations in BR signaling genes points towards the existence of negative feedback of BR signaling on BR biosynthesis. If indeed negative feedback exists between BR signaling and biosynthesis, this feedback should be reflected in quantitative differences in overexpression of the BR signaling and biosynthesis genes in the bri1–5 and suppressor mutants. The better the signaling can be restored in the suppressors (as reflected by the phenotype), the less we expect the expression of the BR biosynthesis to be aberrant. We indeed found that the expression of the BR-biosynthesis genes (CYP90C1, CYP90A1, CYP85A2, CYP90B1, CYP708A3) are less affected in the strains that better mimic the wild-type phenotype (see Fig. 6, the best suppressor of bri1–5, bri1–5/bak1-D, shows the lowest expression change of the biosynthesis genes). This further supports the existence of negative feedback from BR regulation on BR biosynthesis: a more sustained BR signaling results in decreased BR biosynthesis, whereas suboptimal BR signaling is compensated for by higher transcriptional activity of BR biosynthetic genes.
Link between stress response and BR signaling
Next, we analyzed the genes that were uniquely altered in each of the suppressors, i.e. the suppressor-specific compensatory genes (genes of group D, group E, and group F, respectively). Based on GO analysis we mainly found stress related processes to be overrepresented in each of the groups. Genes involved in these stress related processes seem to be scattered over the interaction network as they could not be recovered as well delineated subnetworks, indicating that different stress-related genes are induced in the different lines. According to the literature, there is crosstalk between BR signaling and signaling by other hormones in response to stress, especially via ABA and auxin signaling [28]. In the absence of BRs (or low amounts of BRs), BIN2, affects ABA and auxin signaling, resulting in the induction of stress response genes [4, 5, 15]. On the other hand, some stress-response genes are known to be targets of BZR1 and BES1 [5, 15] (Fig. 1), indicating that also when BRs levels are high, stress response genes can be activated. These observations show that balanced BRs levels are needed for normal growth and that deviation from the optimal levels (either too high or too low) would activate stress response mechanisms. We observed that by partially recovering bri1–5 signaling deficiency by suppressors mutants, the transcript level of some stress-response genes is restored to normal, but other stress-response genes become induced (Additional file 3: GO enrichment for genes exclusively differentially expressed in each suppressor, “GO_only_bri1–5”, “GO_only_bri1–1D”, “GO_only_bak1–1D”, “GO_only_brs1–1D”). This observation is in line with this complex effect of BRs and BR signaling on stress response pathways.
Iron ion homeostasis, ferroxidase, and glutathione transferase activity are identified as compensatory mechanisms unique to the bri1–5/brs1–1D suppressor
Unlike for BRI1 and BAK1, much less is known about the role of BRS1 in BR signaling. Therefore we had a closer look at genes of group D which are exclusively differentially expressed in bri1–5/brs1–1D mutant and hence comprise compensatory pathways specific for bri1–5/brs1–1D. GO enrichment showed that the genes of this group (group D, 333 genes) are not only overrepresented in stress related processes (see above) but also in glutathione transferase (up-regulated), (Fig. 7). This overrepresentation in glutathione transferase is in line with the MapMan results. These results showed how in the bri1–5/brs1–1D mutant the expression of the glutathione transferase was not only restored as compared to the other suppressor lines, but even overcompensated as compared to WS2 levels (Fig. S7). We also found that several members of the CCAAT-binding factor complex (CBC) (NFYA2, NFYA3, NFYA6, NFYA10) were uniquely up-regulated in bri1–1D/brs1–1D (Fig. 7 and Fig. S9). Members of this complex have been associated with the control of iron homeostasis in Candida glabrata [33]. In addition, iron ion homeostasis/ferroxidase activity was also found to be down-regulated, specifically in the bri1–5/brs1–1D. In the ferroxidase reaction, four H+ are used to catalyzes the oxidization of Fe2+ to Fe3+, repressing this reaction results in the accumulation of H+ which can be transported to the apoplast via plasma-membrane pumping mediated by ATPase (H + -ATPase transporters) [34]. Accordingly, we also found that the main inhibitor of H+-ATPase transporters, CBC1, was significantly down-regulated in bri1–5/brs1–1D (fold change − 1.7, adj p-value 8.36e-06), but not in the other suppressors. This implies that H+-ATPase transporters are more active in bri1–5/brs1–1D to export H+ from cytosol into apoplast, making the apoplast more acidic (Fig. 1). In line with this hypothesis, the up-regulated glutathione transferase activity in the brs1–1D mutant (Fig. S9) might be essential to compensate for the more acidic environment in the bri1–5/brs1–1D and would be required for maintaining redox homeostasis. In addition, we hypothesize that the observed acidification could generate a cellular environment that improves BRI1-BRs binding or BRI1-BAK1 dimerization and hence contributes to restore the bri1–5 mutant phenotype.