Network analysis uncovers the master role of WRKY transcription factors in Arabidopsis thaliana response to N-acyl homoserine lactones

Abstract

Background

Plants can perceive bacterial molecules such as the quorum sensing signals N-acyl homoserine lactones (AHL), thus modifying their fitness in response to environmental factors. Even though the benefits conferred by AHL depend on various hormone signaling pathways, the understanding of AHL signaling, especially the response to AHL presence, remains largely unknown.

Methods

Weighted gene co-expression network analysis (WGCNA), multi-omics network analysis, and reverse transcription quantitative PCR (RT-qPCR) assays were used to identify key genes in AHL signaling.

Results

To obtain comprehensive insights on AHL signaling, we integrated available transcriptome data from Arabidopsis thaliana exposed to different single or multiple AHL molecules and performed a weighted gene co-expression network analysis. We identified several key genes regulated in plants exposed to multiple AHL molecules. Multi-omics network analysis and RT-qPCR assay revealed a potential role of WRKY transcription factors.

Conclusions

Results presented here offer good indications for exploring the mechanism of plants' response to bacterial signaling molecules, which could further support the application of AHL-producing bacteria in sustainable agriculture.

Background

N-acyl homoserine lactones (AHL) are prominent quorum-sensing molecules in bacteria-bacteria interactions. In addition, AHL molecules play an important role in interactions between bacteria and plants (Schikora et al. 2016). Although the perception of AHL molecules is not entirely understood, increasing evidence suggests that plants recognize these bacterial signals and modulate specific physiological processes such as growth or responses to biotic and abiotic stressors (Shrestha and Schikora 2020; Babenko et al. 2022). Bacterial AHL molecules consist of a hydrophilic homoserine lactone ring and an acyl chain, which can vary in length from 4 to 18 carbon atoms and may be modified on its third atom with a hydroxyl or a ketone group (Churchill and Chen 2011). Generally, enhanced primary root growth has been a result of treatments with short-chain AHL molecules including N-butyryl homoserine lactone (C4-HSL), N-hexanoyl homoserine lactone (C6-HSL), N-3-oxo-hexanoyl homoserine lactone (oxo-C6-HSL), and N-3-oxo-octanoyl homoserine lactone (oxo-C8-HSL) (von Rad et al. 2008; Liu et al. 2012). Medium-chain AHL molecules modified root architecture, increasing the number of lateral roots and root hairs, as revealed for N-decanoyl homoserine lactone (C10-HSL) (Liu et al. 2012; Cao et al. 2022) and N-3-oxo-decanoyl homoserine lactone (oxo-C10-HSL) (Pazarlar et al. 2020). Long-chain AHL molecules such as N-3-oxo-dodecanoyl homoserine lactone (oxo-C12-HSL) and 3-oxo-tetradecanoyl homoserine lactone (oxo-C14-HSL) enhanced plant resistance to biotrophic and hemibiotrophic pathogens (Schikora et al. 2011; Schenk et al. 2014; Shrestha et al. 2019; Pazarlar et al. 2020). Interestingly, root treatment with oxo-C6-HSL enhanced plant tolerance to salt stress, while promoting plant growth under elevated salt level conditions (Zhao et al. 2020). The application of oxo-C8-HSL enhanced plant resistance to the hemibiotrophic pathogen Pseudomonas syringae pv. tomato DC3000 (Pst) (Liu et al. 2020) and the necrotrophic pathogen Pectobacterium carotovorum ssp. carotovorum (Liu et al. 2022). Notably, the effects of AHL on plants depended also on their concentration (Liu et al. 2012; Palmer et al. 2014; Schenk et al. 2014; Pazarlar et al. 2020). A biphasic growth response to AHL has been shown especially for primary root elongation in Arabidopsis thaliana (Liu et al. 2012; Palmer et al. 2014): low AHL concentration (1 μM to 10 μM) promoted root growth, however, high concentration (100 μM) inhibited root growth. Furthermore, these low and high concentrations were associated with the situation in bulk soil or the rhizosphere and bacterial biofilm formation on plant roots, respectively (Danhorn and Fuqua 2007; Shrestha and Schikora 2020; Babenko et al. 2022).

In nature, plants may host various AHL-producing bacteria, often encountering multiple AHL at once. This phenomenon motivated us to explore the plant response to multiple AHL molecules. Their impact on plant growth and immunity was assessed upon exposure to four different single AHL molecules including oxo-C6-HSL, oxo-C8-HSL, oxo-C12-HSL, oxo-C14-HSL, and their combinations (Shrestha et al. 2020). Compared to responses to four single AHL molecules, their combinations, and especially the mixture of all four AHL molecules (AHL mix) had a rather low impact on the plant growth, however, it induced plant resistance to Pst. Our recent study highlighted the function of jasmonates metabolism in defense responses enhanced by AHL mix (Duan et al. 2023). Similarly, benefits conferred by different single AHL molecules were regulated by diverse phytohormones such as abscisic acid (ABA), jasmonic acid (JA), salicylic acid (SA), and auxin (Schenk et al. 2014; Liu et al. 2020, 2022; Zhao et al. 2020; Duan et al. 2023). In addition, Ca²⁺, reactive oxygen species (ROS), nitric oxide (NO) production, and phenolic compounds participated in responses of A. thaliana to different AHL molecules (Song et al. 2011; Zhao et al. 2014; Pazarlar et al. 2020; Cao et al. 2022).

Since multiple phytohormones as well as the other defense-related compounds are involved in AHL signaling, we sought to explore how plants coordinate their responses under specific conditions. The multifaceted transcription factor AtMYB44 played a positive role in oxo-C6-HSL-mediated primary root growth by regulating cytokinin and auxin-related genes (Zhao et al. 2016). The upregulation of salt‑responsive genes in plant exposed to oxo-C6-HSL could be a possible explanation for enhanced tolerance to salt (Zhao et al. 2020). The inhibitory effects on primary root growth of C10-HSL may depend on enhanced concentration of free cytosolic Ca²⁺ and ROS, as well as the activity of mitogen-activated protein kinase 6 (MPK6) and the formation of NO (Cao et al. 2022).

In addition, transcriptome analysis largely contributed to our understanding of the responses to different AHL molecules (von Rad et al. 2008; Schenk et al. 2014; Zhao et al. 2016; Liu et al. 2022; Shrestha et al. 2022; Duan et al. 2023). On one side, reanalyzing these data with an updated database could be a helpful means to understand the plant responses to AHL molecules even better. On the other side, most studies focused on a binary model based on single AHL molecule. Our recent study included a transcriptome analysis on A. thaliana response to AHL mix, which intended to mimic response of plant encountering multiple AHL molecules (Duan et al. 2023). This study revealed a rather novel response to AHL mix in A. thaliana, even though it still needs further elucidation.

Weighted gene co-expression network analysis (WGCNA) may be used to identify key genes which expression is highly correlated to particular treatments or symptoms (Langfelder and Horvath 2008). In this study, an integration of different data was used to perform following analyses: (i) functional analysis of the re-annotated data, (ii) identification of hub genes in plant responses to AHL molecules via WGCNA, (iii) identification of key genes of AHL mix-signaling via multi-omics network analysis, and (iv) direct comparison between responses to single AHL molecules and AHL mix upon same conditions via an independent RT-qPCR assay. The co-expression network analysis identified novel key (hub) genes, which were highly correlated with response to AHL mix. In addition, the multiple-omics network analysis and gene expression results from RT-qPCR assay suggested a key role of WRKY18, WRKY33, and WRKY40 in AHL signaling. In summary, results presented here improve our understanding of plant response to AHL molecules, which may facilitate the application of AHL-producing bacteria in sustainable agriculture.

Material and methods

Data collection and pre-processing

Available transcriptome data analyzed here originated from two distinct methods: microarray and RNA sequencing. The starting point of the time scale was set to 0 h, time point at which plants were treated with AHL molecule or their mixture.

Data for Arabidopsis thaliana genotype Col-0 (A. thaliana) after root inoculation with 10 µM N-hexanoyl homoserine lactone (C6-HSL) were collected from leaf tissues and root tissues separately 4 h, 24 h, and 96 h post AHL treatment (von Rad et al. 2008).

Data for A. thaliana treated with 1 µM N-3-oxo-hexanoyl homoserine lactone (oxo-C6-HSL) and 10 μM N-3-oxo-octanoyl homoserine lactone (oxo-C8-HSL) were accessed from Gene Expression Omnibus (GEO), accession number: GSE78079 and GSE197485, respectively. The whole seedlings were harvested 24 h post oxo-C6-HSL or oxo-C8-HSL treatments (Zhao et al. 2016; Liu et al. 2022).

Data for A. thaliana pretreated with 6 µM C6-HSL, N-3-oxo-decanoyl homoserine lactone (oxo-C10-HSL), and N-3-oxo-tetradecanoyl homoserine lactone (oxo-C14-HSL) for 72 h, individually, with addition of oxo-C14-HSL-primed plants challenged with 100 nm flagellin-derived 22-amino-acid peptide (flg22) for additional 2 (termed 74 h) and 24 h (termed 96 h) were collected from Schenk et al. (2014). Those data are stored in GEO, accession number: GSE52979.

Data for A. thaliana pretreated with 6 µM oxo-C14-HSL for 72 h, thereafter challenged with 100 nM flg22 for additional 2 h (termed 74 h) were extracted from Schenk et al. (2014). Raw data are stored in GEO with the number GSE156726.

Data for A. thaliana pretreated with the mixture of AHL molecules (AHL mix) including oxo-C6-HSL, oxo-C8-HSL, N-3-oxo-dodecanoyl homoserine lactone (oxo-C12-HSL), and oxo-C14-HSL at equal concentration of 6 µM for 72 h, thereafter challenged with 100 nm flg22 for additional 2 (termed 74 h) and 24 h (termed 96 h) were extracted from Duan et al. (2023). Raw data are stored in Sequence Read Archive (SRA) with the BioProject number: PRJNA865543.

Functional analysis of differentially expressed and common genes

Differentially expressed genes (DEGs) were identified using parameters described in the method chapter from examined studies (von Rad et al. 2008; Schenk et al. 2014; Zhao et al. 2016; Liu et al. 2022; Shrestha et al. 2022; Duan et al. 2023). DEGs of RNA sequencing data were identified using Deseq2 version 1.38.3 (Love et al. 2014) with the threshold (adjusted p < 0.05, fold change of feature counts > 1.5). For microarray data, the DEGs were filtered by threshold (a minimum of 1.6-fold change of fluorescent intensity). For further comparisons and data analysis, all raw RNA sequencing reads were aligned to the A. thaliana reference genome Col-CEN_v1.2 (Naish et al. 2021). The read mapping and feature counting were performed using Rsubread version 2.12.3 (Liao et al. 2019). Furthermore, genes were identified by the A. thaliana Information Resource (TAIR) ID in RNA sequencing data. Thus, the TAIR IDs of the microarray data were obtained by aligning the GenBank accession number (GB_ACC) to Platform GPL12621 in GEO. The identification of enriched gene ontology (GO) term among DEGs was performed using clusterProfiler version 4.0 (Wu et al. 2021) with the following parameters: fun = "enrichGO",ont = "BP",OrgDb = org.At.tair.db, keyType = "ENTREZID".

When preparing the common genes for weighted gene co-expression network analysis (WGCNA), genes from RNA-sequencing data (Duan et al. 2023) with minimum count (reads ≥ 10), in at least 80% of samples were selected. The count data were transformed and normalized by division by their library size using a variance stabilizing transformation (VST) function in DEseq2 version 1.38.3 (Love et al. 2014), which finally yielded in a matrix of constant variance along the range of mean values. The quantile normalized data (Schenk et al. 2014) was used for analysis of plant response to three single AHL molecules including C6-HSL, oxo-C10-HSL, and oxo-C14-HSL. Systematic biases within the two data sets were further standardized using the quantile method (Ritchie et al. 2015). Finally, 19,085 common DEGs were integrated from 20 samples, which were harvested 72 h post AHL treatment (Schenk et al. 2014; Duan et al. 2023).

Weighted gene co-expression network analysis (WGCNA)

The WGCNA package version 1.72–1 (Langfelder and Horvath 2008) was employed to construct the network, using the common gene set with the default setting, in which the desired minimum scale-free topology fitting index was set to 0.85. The detection of modules and calculation of the correlations between module and traits were performed using the function “blockwiseModules” in WGCNA package. The following parameters were used: power = 28, maxBlockSize = 16,000, minModuleSize = 30, reassignThreshold = 0, mergeCutHeight = 0.25, numericLabels = TRUE. The core (hub) genes of each module were identified using following threshold: module membership (MM) greater than 0.8 and gene significant (GS) greater than 0.4. The functional analysis of hub genes was performed using the same setting as for DEGs.

Validated predictions and interaction partners of the hub genes were extracted from (Depuydt and Vandepoele 2021), detailed information was listed in Additional file 5: Data Set 5. The visualization of network was performed in the open-source software platform Cytoscape (version 3.9.1).

Gene expression assays

Two-week-old A. thaliana were transferred from half-strength MS (Murashige and Skoog) (½ MS) medium to six-well plates with ½ MS liquid medium. Subsequently, 4 different AHL molecules (Sigma–Aldrich Chemie GmbH, Munich, Germany): oxo-C6-HSL, oxo-C8-HSL, oxo-C12-HSL, and oxo-C14-HSL at 6 µM each, were added separately. In addition, AHL combination (AHL mix) was added to the growth medium at 6 μM concentration of each four different AHL molecules. The same volume of solvent (acetone) was added to control samples. Plants were grown in six-well plates for additional 72 h, and challenged with 100 nM flg22 for another 2 h (74 h) and 24 h (96 h). Samples were harvested from six-well plates directly, each well contained 40–60 seedlings, which was defined as one biological sample (replicate). Each treatment had four independent biological replicates. The extraction and purification of RNA and analysis of gene expression using reverse transcription quantitative PCR (RT-qPCR) were performed using pipeline described in Duan et al. (2023). Primers and accession numbers of genes were listed in Additional file 6: Table S1.

Results

Plant responded to N-acyl homoserine lactones in an intricate manner

To reveal main characteristics of plant response to various AHL molecules (Fig. 1A), multiple data available from open repositories were re-annotated and compared with each other. Detailed information on data accession, processing, and re-annotation with the updated reference data was described in the method section. Differentially expressed genes (DEGs) were identified using parameters described in method sections in original studies (von Rad et al. 2008; Schenk et al. 2014; Zhao et al. 2016; Liu et al. 2022; Shrestha et al. 2022; Duan et al. 2023). In order to obtain a consistent time schedule for comparison of different results, the time point of AHL treatment was set to 0 h, which resulted in 2 h, 4 h, 24 h, 72 h, 74 h, and 96 h post treatment (hpt), regardless a possible additional challenge.

Fig. 1

Complex interactions may occur between plants and N-acyl homoserine lactones. Integration of available transcriptome results has taken into account following factors: N-acyl homoserine lactone (AHL), time point, and tissue type. The impact of different AHL molecules (A) and their mixture (AHL mix) on A. thaliana has been assessed in several samples including seedlings, root tissues, leaf tissues, and ali1 mutant (von Rad et al. 2008; Schenk et al. 2014; Zhao et al. 2016; Liu et al. 2022; Shrestha et al. 2022; Duan et al. 2023). In order to compare results from different sources, the time point of AHL treatment was set to 0 h. Notably, in some studies, plants pretreated with AHL molecules for 72 h were challenged with 100 nM flg22 for another 2 h (74 h) and 24 h (96 h). The enriched gene ontology (GO) terms with the biggest gene counts were selected as representative for plant responses to AHL molecules (B)

In the next step, enriched gene ontology (GO) terms in obtained DEGs sets were identified, those with the biggest counts were referred to as the representative GO term in the particular set (Fig. 1B, Additional file 1: Data Set 1). Exposure to AHL molecules such as oxo-C6-HSL and oxo-C8-HSL, which have similar chemical structure, seemed to enrich GO terms related to response to cold and metal ions 24 hpt, respectively (Zhao et al. 2016; Liu et al. 2022). Treatment with different AHL molecules was followed by enrichment of GO terms related to other processes at 72 hpt: (i) lipid localization in case of C6-HSL and oxo-C14-HSL, (ii) response to hypoxia in case of oxo-C10-HSL (Schenk et al. 2014). Different plant tissues such as root or leaf displayed distinct biological responses when root was exposed to C6-HSL (von Rad et al. 2008). In root tissues, the enriched GO terms were related to endomembrane system organization, generation of precursor metabolites and energy, and circadian rhythm (at 4, 24, and 96 hpt), whereas none GO term was enriched in leaf tissues 4 hpt. The GO term “response to jasmonic acid” was enriched 24 and 96 hpt (Fig. 1B).

AHL priming for enhanced resistance initiated many of the studies on plant responses to ALH molecules. Therefore, we reanalyzed also the response to a secondary flg22 challenge. Generally, plants exposed to AHL molecules revealed different responses to flg22 challenge, than naïve plants (Schenk et al. 2014). GO terms related to cell cycle and triterpenoid metabolism were enriched in oxo-C14-HSL-pretreated plant, 2 h after flg22 challenge, 74 h after AHL treatment (Shrestha et al. 2022). Secondary metabolism related GO terms were enriched in oxo-C14-HSL-primed plants 96 hpt (24 h after the flg22 challenge). Compared to the response to various single AHL molecules, plant response to AHL mix revealed a novel pattern across the time scale, from 2 to 96 h (Duan et al. 2023). Here, the GO term (response to fungus) was enriched 2 hpt and 96 hpt, whereas the GO term (defense response to bacterium) was enriched 72 hpt and 74 hpt. In addition, rRNA processing-related GO terms were enriched 24 hpt. These results revealed that the interactions between AHL molecules and plants are intricate, and may be affect by many factors such as the type of tissue (root tissues, leaf tissues, seedlings), the time point, the AHL chemical structure, and probably others.

Identification of hub genes in plant response to N-acyl homoserine lactones

Generally, AHL signaling processes were clustered into two different biological processes: i) plant response to AHL molecules, and ii) response of AHL-primed plants to a secondary challenge, such as flg22 or the other stressors. Compared to the binary model system based on AHL molecules and plant, addition of a secondary challenge significantly increase the complexity. Since we focused on plant response to AHL molecules, the following analysis was based on data from plants exposed to AHL molecules for 72 h, this duration has been successfully used in previous studies to prime plants for enhanced resistance (Schenk et al. 2014; Schikora et al. 2016; Shrestha et al. 2019, 2020, 2022; Duan et al. 2023). In total, the integration of all the data included 20 samples with 19,085 DEGs (Additional file 2: Data Set 2).

Weighted gene co-expression network analysis (WGCNA) was performed in order to identify genes with strong correlation to AHL treatment (Langfelder and Horvath 2008). In the present study, the WGCNA analysis resulted in 8 modules, with different correlation and p values for different AHL molecules and the AHL mix (Fig. 2A), including from 64 to 15,374 genes per module (Fig. 2B). Module (ME) blue (n = 468 genes) was specific in case of the solvent control (CTL) with a positive correlation coefficient (R² = 0.66). Whereas the MEred (n = 86 genes) was specific in the case of oxo-C10-HSL (C10) treatment with a negative correlation coefficient (R² = -0.55), MEgreen (n = 241 genes) and MEturquoise (n = 15,374 genes) were specific for AHL mix (MIX) with a positive correlation coefficient R² = 0.66 and R² = 0.63, respectively. Notably, MEyellow was also specific for AHL mix with a high negative correlation coefficient (R² = −0.96) and a rather small p-value (4e-11). No module was specific in the case of treatments with C6-HSL (C06) and oxo-C14-HSL (C14). The further identification of hub genes in the five correlated modules was performed using an intra-modular analysis based on the module membership (MM) and gene significance (GS). MEyellow had the highest correlation coefficient (R² = 0.73) between MM and GS among the five modules, while the other four modules had a rather lower coefficient, ranging from 0.25 to 0.41 (Fig. 2C).

Fig. 2

Identification of hub genes in plant responses to N-acyl homoserine lactones. Integrated data set containing 19,085 differentially expressed genes in 20 samples, collected from whole seedlings, 72 h post-treatment with different AHL molecules including C6-HSL (C06), oxo-C10-HSL (C10), oxo-C14-HSL (C14), and AHL mix (MIX) (Schenk et al. 2014; Shrestha et al. 2022; Duan et al. 2023), was used to perform weighted gene co-expression network analysis (A). This analysis revealed the relationships (A) between AHL (trait) and gene clusters (module), the number of genes for each module was presented in a barplot (B). Hub genes were identified using thresholds (|gene significance| > 0.4 and |module membership| > 0.8) from five modules, significantly correlating with the AHL treatment or the solvent control (CTL), (p value < 0.05). Hub gene sets identified in five modules (MEblue, MEgreen, MEred, MEturquoise, and MEyellow) were termed blueCTL, greenmix, redC10, turquoisemix, and yellowmix, respectively (C)

Furthermore, a functional enrichment analysis was performed with the identified core genes, which were filtered by threshold (MM > 0.8, |GS|> 0.4) (Figs. 2C and 3A). The GO terms related to response to organonitrogen compound and response to nitrogen compound were enriched in the core genes (termed blueCTL) filtered from MEblue. The core genes (termed yellowmix) filtered from MEyellow and blueCTL enriched in some common GO terms including response to chitin, response to hypoxia, and response to decreased oxygen levels. Notably, functions of core genes (yellowmix) were enriched in GO terms related to defense response to Gram-negative bacteria. The core genes filtered from the MEgreen (termed greenmix) and MEturqoise (termed turqoisemix) were both enriched in GO terms related to pattern specification process. The GO terms related to xylem and phloem pattern formation were enriched in the greenmix group, whereas the GO terms related to ribonucleoprotein complex biogenesis, ribosome biogenesis, and mRNA metabolic process were enriched in the turqoisemix group. Surprisingly, response to UV or UV-B were enriched in the core genes (redC10) filtered from MEred. The GO term related to response to UV-B was also enriched in the turqoisemix group. Core genes grouped in yellowmix drew our special attention since, they displayed high correlation coefficients in both module-trait relationship and MM-GS relationship. Furthermore, the GO term “defense response to Gram-negative bacteria” was enriched within this group, which may indicate that those core genes may play an essential function in the plant responses to AHL mix (Fig. 3).

Fig. 3

Functional analysis of identified hub genes. Comparisons of enriched gene ontology (GO) terms were performed in identified hub genes (A). Expression pattern of hub genes upon AHL mix priming and following challenge with 100 nM flg22 for 2 or 24 h (74 h; 96 h) was extracted from Duan et al. (2023). Differentially expressed genes were marked with symbol “ + ”, this includes genes with absolute value of fold change ≥ 1.5 between the solvent control and AHL mix, whereas genes with absolute value of fold change < 1.5 were marked with “– “. According to the gene expression pattern, only genes correlated to AHL mix were clustered into 8 groups (B). Distribution of differentially regulated genes within each hub gene group (greenmix, turquoisemix and yellowmix) was highlighted with green color (B). Within the yellowmix set, which is strongly correlated to AHL mix, genes related to the GO term: defense response to bacterium (GO: 0042742) represented the biggest group, followed by defense response to fungus (GO: 0050832) and regulation of response to external stimulus (GO: 0032101). Those three GO terms shared 55 core genes (C)

Multi-omics network analysis identified key genes in AHL mix-signaling

After identification of hub gene sets which are highly correlated with AHL mix 72 hpt, we wondered how these genes were regulated after an additional flg22 challenge. Gene expression data 72, 74, and 96 hpt (0, 2 and 24 h after flg22 challenge, respectively), were extracted from previously published RNA-Seq data (Duan et al. 2023). Gene expression patterns could be roughly clustered into eight groups (Fig. 3B and Additional file 3: Data Set 3). Briefly, most of genes in the greenmix (61.3% of total 119 genes) and turquoisemix groups (87.8% of 12,224 genes) were not regulated at these time points, Group VIII (Fig. 3B). Whereas most genes of the yellowmix group (66.4% of 250 genes) were differentially regulated before (72 hpt) and after (74 and 94 hpt) the flg22 challenge, groups I and VII (Fig. 3B). Counting the number of genes in relevant GO terms of yellowmix group revealed that the defense response to bacterium (GO:0042742; 84 genes), defense response to fungus (GO:0050832; 72 genes), and regulation of response to external stimulus (GO:0032101; 70 genes) were most prominent (Fig. 3C and Additional file 4: Data Set 4). Furthermore, 55 genes were shared between those GO terms (Fig. 3C), the predominant genes (94% of 84 genes) grouped in the GO term “defense response to bacterium” (GO: 0042742) were differentially regulated before and after the flg22 challenge (Fig. 3B). Thus, we focused on the 84 hub genes belonging to the GO term: “defense response to bacterium” (GO: 0042742), which were also related to the above analysis (Figs. 1 and 2).

To reveal even more details on the regulation of those 84 genes, annotated genes, especially from group I and II were selected and represented with their expression levels (Fig. 4A). Generally, expression of those genes was downregulated 72 hpt with AHL mix, and upregulated 2 h after the flg22 challenge (74 hpt). In addition, upregulation was present even 96 hpt in group I, while this phenomenon was not present for genes in group II. Notably, gene ZF14, which encodes a MATE efflux family protein, showed an opposite expression pattern: upregulation 72 hpt, and downregulation 74 and 96 hpt. To broaden our understanding of plant's initial responses to the AHL mix, the hub gene expression data (Duan et al. 2023) were also extracted from samples harvested 2 hpt and 24 hpt with AHL mix (Fig. 4). Most of these genes were not regulated either 2 or 24 hpt. Some few genes including GSTF6/7 from group I, CMPG2, and GSTU12 from group II were upregulated 2 hpt with AHL mix. Unlike genes from group I and II, expression of genes from group IV and VI was altered 72 and 74 hpt, respectively; whereas expression of genes in group VIII seemed to be not regulated after AHL mix treatment (Fig. 4B).

Fig. 4

Expression of genes within the GO term: defense response to bacterium. Expression patterns of the 84 hub genes were assessed from the GO term: defense response to bacterium (GO: 0042742). Those genes were identified in the module MEyellow, extracted from Duan et al. (2023). Gene expression was presented as heatmap (comparison between AHL mix and solvent control in AHL mix-primed plant and plants challenged with 100 nM flg22). Annotated genes were presented in (A) meanwhile, genes with unknown function or not regulated were presented in (B)

According to annotated protein function (Additional file 4: Data Set 4), the 84 hub genes were clustered into following groups: (i) up to 14 transmembrane proteins and receptor (-like) kinase genes, which may be involved in the perception of extracellular signals. For instance, in group I, FLS2 encodes a leucine-rich receptor-like protein kinase, which is an intensely studied plant pattern recognition receptor mediating plant response to bacteria; (ii) sixteen genes related to transcription factors and transporter, which could function in intracellular signal transduction. In group II, WRKY18 encodes a WRKY DNA-binding protein 18, which is a known pathogen-induced transcription factor (Birkenbihl et al. 2017); (iii) twenty-three secondary metabolite-related genes, which may be involved in the modification of hormone levels of signaling pathways, such as CEJ1, cooperatively regulated by ethylene and jasmonate, or BCS1 and ERF4, related to salicylic acid and ethylene, respectively; (iv) genes with less known functional characteristics.

The high complexity of plant responses to AHL molecules (Fig. 1) motivated us to explore whether key regulators in plant response to AHL molecules may be identified. To this end, we used the identified hub genes (Fig. 4) and build a model of plant response to AHL mix. Firstly, biological connections among the hub genes were identified based on known protein-DNA (PDI) and protein–protein interaction (PPI) data (Depuydt and Vandepoele 2021). Secondly, we added predicted functionally-related genes, in order to extend the hub gene matrix to known interactions (Depuydt and Vandepoele 2021). Noteworthy, changes in the expression of added genes after AHL mix treatment were also verified (Additional file 5: Data Set 5). In total, 138 genes (84 hub genes + 54 predicted genes) were used to model plant response to AHL mix (Fig. 5). The importance of genes was assessed by counting the number of their biological connection with other genes (Additional file 5: Data Set 5). The size of the node indicates the number of interconnections in the model. In addition, gene expression patterns (groups I to VIII) were indicated by colors. Three transcription factor genes WRKY18, WRKY33, and WRKY40 of group II seemed to play a central role in this model. This indicated that regulation of those WRKYs may be involved in the plant responses to AHL mix and further in the response of AHL mix-primed plant to the challenge with flg22. Another prominent node was occupied by MYBR1/MYB44 from group VIII, even though the expression of this transcription factor was not regulated at any verified time points. In summary, these results revealed that three WRKYs may function as key regulators in the AHL mix-signaling.

Fig. 5

Multi-omics network modeling of the response to AHL mix. Protein-DNA (PDI) and protein–protein interaction (PPI) data (Depuydt and Vandepoele 2021) were used to obtain the biological connections among identified hub genes (GO: 0042742) and their predicted functionally-related genes. The size of nodes represents the number of interactions between genes, the color of nodes represents the expression pattern (groups I to VIII). The relationship between genes is represented using three types of lines: solid (PPI + PDI), dots (PDI), and marquee dash dot (PPI). The network was created in Cytoscape version 3.9.1

WRKYs function in response of AHL-primed plants to flg22

The multi-omics network (Fig. 5) highlighted the potential role of WRKY transcription factors in A. thaliana response to AHL mix, which motivated us explore the regulation of WRKY using an independent approach. To this end, two-week-old plants were treated with either 6 µM oxo-C6-HSL, oxo-C8-HSL, oxo-C12-HSL, and oxo-C14-HSL or their mixture (AHL mix) for 72 h. Thereafter, plants were challenged with 100 nm flg22 for 2 h (74 hpt) and 24 h (96 hpt). The determination of gene expression (Fig. 6) was performed for four WRKY including WRKY18, WRKY33, WRKY40, and WRKY60, as well as ZF14 which revealed an oppositive expression pattern (Fig. 4), if compared to WRKY expression in AHL mix-primed plants.

Fig. 6

Regulations of WRKY in various AHL-primed plants after challenge with flg22. Two-week-old A. thaliana was treated with either 6 µM oxo-C6-HSL, oxo-C8-HSL, oxo-C12-HSL, oxo-C14-HSL, or their mixture (AHL mix) for 72 h. Total RNA was extracted prior the challenge (72 h) as well as 2 h (74 h) and 24 h (96 h) after an additional challenge with 100 nM flg22. Expression of ubiquitin ligase (At5g25760) was used for normalization. Boxes represent the interquartile range between the first and the third quartile and the middle line marks the median, whiskers indicate 1.5 × interquartile range. Each treatment and time point had four biological replicates. Comparisons are indicated by the horizontal lines above bars, statistical analysis was performed with Student’s t-test, *p < 0.05, **p < 0.01, ***p < 0.001, ****p < 0.0001

Seventy-two hours after the AHL treatment, the four WRKY and ZF14 were significantly upregulated particularly in oxo-C12-HSL, oxo-C14-HSL, and AHL mix pretreated plants (Fig. 6). Lower upregulation levels of WRKY33 and WRKY40 were observed in the oxo-C6-HSL-primed plants (Fig. 6A). After the flg22 challenge, the transcriptome reprogramming was rather dynamical (Fig. 6). For instance, WRKY33 was upregulated 74 hpt and downregulated 96 hpt in oxo-C12-HSL, oxo-C14-HSL, and AHL mix-primed plants, similar phenomenon was observed for WRKY18 and WRKY 40 (Fig. 6A). In contrast, WRKY60 was downregulated or not regulated 74 hpt (Fig. 6B). We noticed that the regulations of WRKY displayed differences in response to the four single AHL molecules and their mixture (AHL mix) (Fig. 6A), which further underlined the complexity of plant response to various AHL molecules (Fig. 1). In addition, different expression patterns between response to long-chain and to short-chain AHL molecules seemed apparent (Fig. 6). For instance, expression of WRKY33 in response to long-chain AHL molecules (oxo-C12-HSL and oxo-C14-HSL) was upregulated 74 hpt (2 h after flg22 challenge), whereas the same phenomenon was not observed in response the short-chain AHL molecules including oxo-C6-HSL and oxo-C8-HSL. Furthermore, ZF14 was upregulated in response to all four single AHL molecules but not in AHL mix-primed plants (74 hpt; 2 h after flg22 challenge), later its expression was downregulated in response to various AHL molecules, except for oxo-C6-HSL (Fig. 6B). These results demonstrated that WRKY may be involved in response to various AHL molecules and their mixture.

Discussion

Fastly accumulating evidence indicates that eukaryotes recognize bacterial quorum sensing signals such as N-acyl homoserine lactones (AHL) and adapt their own physiology (Shrestha and Schikora 2020; Wu and Luo 2021; Babenko et al. 2022). Compared to the well-characterized AHL receptors in bacteria (McCready et al. 2019; Wellington and Greenberg 2019), little is known about AHL perception in plants. Previous studies have revealed several proteins which may be involved in AHL-induced signaling cascades related to plant growth and innate immunity systems. For instance, the R2R3-type MYB transcription factor AtMYB44 was required for modulation of the primary root growth, induced by N-3-oxo-hexanoyl homoserine lactone (oxo-C6-HSL) (Zhao et al. 2016). The membrane-associated protein glucuronokinase 2/ AHL-priming protein1 (AtGlcAK2/ALI1), which belongs to the GHMP (galactokinase, homoserine kinase, mevalonate kinase and phosphomevalonate kinase) kinase superfamily, was identified as integral component of N-3-oxo-tetradecanoyl homoserine lactone (oxo-C14-HSL)-mediated enhanced resistance in A. thaliana (Shrestha et al. 2022). In addition, an acidic Ca²⁺ binding protein calmodulin (AtCaM), G-protein-coupled receptor GCR1, and the canonical Gα subunit GPA1 were reported to participate in root growth induced by oxo-C6-HSL and N-3-oxo-octanoyl homoserine lactone (oxo-C8-HSL) (Liu et al. 2012; Zhao et al. 2014). A rather different action of AHL molecules, dependent on their hydrolysis by a plant-derived fatty acid amide hydrolase to generate free L-homoserine, was proposed by Palmer et al. (2014).

Notably, all above findings investigated binary models based on the response to single AHL molecules. Comprehensive analysis of plant responses to structurally diverse single AHL molecules revealed that those molecules affect plant physiology in different aspects, including the growth and plant tolerance to biotic or abiotic stressors (Fig. 1). Our previous study suggested a model based on A. thaliana and four equivalent concentrations of AHL molecules (AHL mix) varying in the acyl chain length, mimicking plants' natural responses to multiple AHL molecules (Duan et al. 2023). Transcriptome analysis revealed a novel expression pattern in AHL mix primed-plants, compared to plant responses to single AHL molecules. Phytohormone measurements and genetic approaches highlighted the importance of jasmonates catabolism in AHL mix-signaling. However, key (hub) genes, which would drive the novel expression pattern remained unidentified. In the present study, we integrated genes differentially expressed in plants 72 h after exposure to three single AHL molecules including, N-hexanoyl homoserine lactone (C6-HSL), N-3-oxo-decanoyl homoserine lactone (oxo-C10-HSL), and oxo-C14-HSL, as well as AHL mix including oxo-C6-HSL, oxo-C8-HSL, N-3-oxo-dodecanoyl homoserine lactone (oxo-C12-HSL), and oxo-C14-HSL (Schenk et al. 2014; Duan et al. 2023). Weighted gene co-expression network analysis of the common genes identified three modules highly correlated to AHL mix, in which function of 84 genes, related to plant defense response, were roughly clustered into four groups: (i) transmembrane or receptors-like proteins to perceive extracellular stimulus; (ii) transporter related proteins in metal ion homeostasis; (iii) transcription factors to modulate the defense-related gene expression; (iv) biosynthesis, modification, and responsive related proteins to optimize plant physiology (Fig. 4). Furthermore, the comparison between responses to four single AHL molecules and the AHL mixtures demonstrated that the AHL mix appears to predominantly modulate plant immunity (Shrestha et al. 2020; Duan et al. 2023).

One of the main interests in AHL priming is the faster and stronger response after perception of various conserved microbe-associated molecular patterns (MAMPs), thus successfully inducing defense responses. The flagellin-derived peptide flg22 is a classic MAMP (Felix et al. 1999), facilitating the study on mechanism of enhanced plant resistance induced by AHL (Schenk et al. 2014; Shrestha et al. 2022; Duan et al. 2023). To further explore the regulatory network, we identified three WRKY transcription factors: WRKY18, WRKY33, and WRKY40, which seemed to play a central role in AHL mix-induced network (Figs. 4 and 5). All of them have been reported to function in an early MAMP-triggered immunity sub regulatory network (Birkenbihl et al. 2017). Protein levels of 14 WRKYs including WRKY6, 8, 11, 18, 22, 25, 28, 29, 30, 33, 40, 48, 53 and 75 were elevated more than twice upon the flg22 challenge (Birkenbihl et al. 2018). Several WRKY transcription factors and cis-regulatory elements were coregulated in the response to the flg22 challenge (Zhang et al. 2022). Additionally, many WRKYs are involved in the regulation of various phytohormone signaling pathways. For instance, WRKY33 played an essential role in salicylic acid- and jasmonic acid-related response toward Botrytis cinerea (Birkenbihl et al. 2012), and the receptor-like kinase ERECTA-mediated plant immune response (Cai et al. 2021). Plant sensitivity to ABA seemed also positively regulated by WRKY18 and WRKY60, and the overexpression of WRKY18 and WRKY60 prevented seed germination and root development but triggered plant response to abiotic stressors (Chen et al. 2010). WRKY40 and GOLDEN2-LIKE transcription factors negatively regulated ABA response in A. thaliana (Chen et al. 2010; Ahmad et al. 2019). Interestingly, the delayed upregulation of WRKY60 was observed in ABA-treated plants (Chen et al. 2010), this phenomenon was also observed in AHL-pretreated plants when challenged with flg22 (Fig. 6). Even though the multifaceted role of MYBR1/MYB44 transcription factor has been reported for responses to biotic and abiotic stress factors (Qiu et al. 2019; Zhao et al. 2022), the expression of MYBR1/MYB44 was not significantly regulated in AHL mix-primed plants at any of the analyzed time points (Fig. 5 and Additional file 6: Data Set 5). Thus, these findings further strengthen the notion that WRKY18, WRKY40, and WRKY33 may play key regulatory functions in AHL signaling.

Certainly, plants cannot distinguish between pathogenic or potentially beneficial bacteria based on the structure of produced AHL molecules. On the other side, both the invasion of pathogens and the mutualistic establishment of beneficial bacteria need to interfere with or bypass the plant's immune responses (Nishad et al. 2020; Teixeira et al. 2021). Interestingly, the downregulation of plant defense-related genes in the AHL mix-primed plants indicated that the presence of multiple AHL molecules may be correlated with a suppression of plant immune responses (Fig. 4). However, the fast conversion from suppression to activation requires further investigations.

Conclusions

N-acyl homoserine lactones (AHL) play an important role in communication between AHL-producing bacteria and their host plants. Compared to single AHL molecules, the mixture of AHL (AHL mix) regulates plant gene expression in a different manner. Since, plants may host different AHL-producing bacteria, we focused here on the response to multiple AHL. We integrated available transcriptome data and performed a weighted gene co-expression network analysis, identifying 84 key genes. Furthermore, selected key genes were used to build a network, based on their protein-DNA and protein–protein interactions, which revealed the central role of WRKY transcription factors in AHL mix signaling including WRKY18, WRKY33, and WRKY40. The exact function of the identified WRKYs in AHL signaling, which remains to be elucidated. Redundancy and epigenetic regulation could be a part of the regulation as well. Therefore, further investigations are required to reveal the exact mechanism and the potential role of WRKY transcription factors in the application of AHL-producing bacteria in sustainable agriculture.