Genome-wide mapping of putative PTI-related distal enhancers in Arabidopsis thaliana

To screen for PTI-related enhancers in Arabidopsis, we first generated an ATAC-Seq library using leaf tissues from wild-type Col-0. Using the common cutoff method, we identified 4460 peaks localized in intergenic regions, which were defined as putative distal enhancers (Additional file 1: Table S1), and the remaining peaks were located in gene-related regions (Additional file 2: Fig. S1a). Given that chromatin accessibility may vary under different experimental conditions, we compared our enhancers with the enhancers identified in Arabidopsis seedling leaves from two previous studies to refine our distal enhancer library [9, 33]. Based on different sequencing methods, the enhancers from the three different sources were merged separately, and we obtained 4702 ATAC-seq-based enhancers (hereafter referred to as enhancersA) and 7515 DNase-seq-based enhancers (hereafter, enhancersD) (Additional file 3: Table S2). The average lengths of enhancersA and enhancersD after merging were 283 bp and 349 bp, respectively (Additional file 2: Fig. S1b). Approximately 61.5% of the enhancers obtained from ATAC-seq overlapped with distal DHS regions (Additional file 2: Fig. S1b). Thus, we constructed a putative enhancer library of Arabidopsis in the vegetative growth stage based on distal ACRs and DHSs.

To identify putative distal enhancers involved in PTI responses and immune gene reprogramming, we examined whether enhancer transcripts could be used as active markers. rRNA depletion RNA-seq, which has been proved to be a practical and effective method to identify transcribed enhancers [34], was performed to obtain a global view of eRNA transcripts to assess enhancer expression levels with/without flg22 treatment. We found that active transcription occurred in ∼50% of the genome, corresponding to 336 million reads (>14 Gb), 12% of which were localized in intergenic regions for each sample. We defined potential enhancer transcript regions (see “12” for details) and mapped the rRNA depletion RNA-seq reads. The average expression levels of enhancer transcript regions were calculated as transcripts per million (TPM) values. These enhancers with detectable expression in our analysis were defined as transcribed enhancers, which comprised approximately 25% of total putative distal enhancers. A total of 1079 enhancersA and 1547 enhancersD were expressed without flg22 treatment, while 1102 enhancersA and 1509 enhancersD were expressed 1 h post-flg22 stimulation (Fig. 1a, Additional file 4: Table S3). The transcribed enhancers displayed higher chromatin accessibility than the nontranscribed enhancers, promoters, and random sequences (Additional file 2: Fig. S1c). Then, we calculated the flg22-induced expression differences in individual enhancer transcript regions and identified 697 or 936 upregulated transcribed enhancers from the two libraries (Fig. 1b). These data indicate that a portion of enhancers are actively transcribed and regulated by flg22-mediated immune signaling.

Fig. 1

Distal enhancer transcription were responsive to flg22 treatment in Arabidopsis. a Numbers of transcribed distal enhancer before and after flg22 signaling activation. Trans-enhancersA, ATAC-seq-based transcribed enhancers; Trans-enhancersD, DNase-seq-based transcribed enhancers. b Volcano Plots show differentially expressed transcribed enhancersA and transcribed enhancersD upon flg22 treatment. The numbers of up (x-axis > 0, warm orange) and down (x-axis < 0, blackish green) transcribed enhancers are shown, respectively. c Heatmaps of transcript levels at transcribed enhancers, transposons, promoters (proximal ATAC-seq peaks), and random sequences. Transcript level with and without flg22 treatment were calculated at and around (±0.5 kb) the sequence candidates. The color scales are in BPM for transcript level. d Distribution of the mean distance from enhancers to the nearest RNA Pol II peaks. e, f Distribution of the signals (BPM scale) derived from GRO-seq (e) and pNET-seq (f) at transcribed enhancers, nontranscribed enhancers, and random sequences. Left and right panels represent ATAC-seq-based enhancers and DNase-seq-based enhancers, respectively

Global characterization of Arabidopsis enhancer transcripts

To investigate the characteristics of enhancer transcripts in Arabidopsis, a sliding window approach was used to identify candidate intergenic regions enriched for eRNA expression. We identified 454 and 543 consecutive transcripts as intact enhancer eRNAs at 0 and 1 h, respectively (Additional file 2: Fig. S1d). The mean length of the identified eRNAs was approximately 200 bp. Only a small portion of the intergenic transcripts and eRNAs could be found in the poly(A) RNA-seq datasets (Additional file 2: Fig. S1d), which suggested that most of the transcribed enhancers were nonpolyadenylated. To study the directionality of Arabidopsis enhancer transcription, we analyzed the distribution of RNA-seq reads on both sides of transcribed enhancers. The strand-specific RNA-seq data obtained in this study provided potentially valuable information for mapping sense and antisense transcripts. Heatmap analysis showed that the expression levels of transcribed enhancers were higher than those of transposons (TEs) and random sequences but lower than those of promoters (Fig. 1c). Approximately 65% of the transcribed enhancers showed bidirectional transcription characteristics, while others were only detected on one strand. These results also showed that flg22 treatment influences the eRNA transcription level of a subset of distal enhancer regions.

To further confirm the authenticity of the enhancer transcription events that we detected, we also cited some published data for analysis [35, 36]. Although most of the identified distal enhancers were located less than 5 kb from the nearest putative transcription start sites (TSSs), the expressed enhancers were closer to the RNA Pol II peaks than the nontranscribed enhancers (Fig. 1d, Additional file 2: Fig. S1e). In particular, approximately 20% of the transcribed enhancers overlapped with RNA Pol II binding peaks, which was not the case for nontranscribed enhancers. Furthermore, the transcription signals derived from GRO-seq and plant native elongating transcript sequencing (pNET-seq) were significantly stronger around transcribed enhancers than around nontranscribed enhancers and random sequences (Fig. 1e, f). The pausing of RNA Pol II and production of nascent RNAs around enhancers further provided direct and reliable measurements of transcribed enhancer activity. Additionally, four reported jasmonate-related enhancers [37] showed open chromatin and increased transcription in response to flg22 (Additional file 2: Fig. S1f). Taken together, our data identified transcribed enhancers, which were further verified by multiple lines of evidence.

Global characterization of transcribed enhancers in Arabidopsis

The observation of increased chromatin accessibility at transcribed enhancers prompted us to evaluate other chromatin and epigenetic modification characteristics. We first investigated the relative enrichment of a list of histone modifications and the histone variant H2A.Z at enhancers based on public datasets were generated from Arabidopsis vegetative tissues [9, 38,39,40,41]. The H3K9ac, H3K23ac, H3K27ac, and H3K56ac levels of transcribed enhancers were lower than those of promoters but higher than those of nontranscribed enhancers and random sequences (Fig. 2a). Among these marks, H3K9ac and H3K27ac are associated with both promoters and enhancers in metazoans and are conserved marks allowing the prediction of proximal DHSs in Arabidopsis [42]. Interestingly, transcribed enhancers were enriched with H3K14ac, which showed even higher levels than in promoter regions (Fig. 2a). In contrast, transcribed enhancers were less associated with the inactive chromatin-associated modification H3K27me3 than nontranscribed enhancers and random enhancers (Fig. 2a). H3K4me1 is specifically enriched at enhancers in metazoans [18] but shows low levels in both transcribed and nontranscribed enhancers in Arabidopsis. H3K9me2 and H3K27me1 levels in the central region of enhancers were higher than those in the flanking sequences and promoters but lower than those in random sequences. In addition, a slightly reduced level of H2A.Z was related to an increase in enhancer transcription activity. Further insight into chromatin characteristics was obtained by computing chromatin state predictions for Arabidopsis intergenic regions based on the relative enrichment levels of these ten histone modification marks by using ChromHMM [43]. The results showed that enhancers were highly related to a silenced state, while promoters were more closely related to a highly active state (Additional file 2: Fig. S2a). Importantly, compared to nontranscribed enhancers, the transcribed enhancer-associated chromatin state was more active.

Fig. 2

Transcribed enhancers represent a more active and conserved class of distal cis-regulatory elements. a Distribution of different histone modifications and the histone variant H2A.Z at transcribed enhancers, nontranscribed enhancers, promoter (proximal ATAC-seq peaks), and random sequences. b DNA methylation level (mC, weighted average) profiles around transcribed and nontranscribed enhancers. c Correlation between DNA methylation level (mC, weighted average) and chromatin accessibility of all transcribed enhancers is represented graphically by a scatterplot. d Distribution of the unmethylated regions (UMRs) in Arabidopsis genomic regions. The total number of UMRs and genome proportion of each part are shown. e Overlap of the total length of enhancers relative to UMRs. The percentage of enhancers overlapping with the UMRs was listed. I, transcribed enhancersA; II, nontranscribed enhancersA; III, transcribed enhancersD; IV, non-transcribed enhancersD. f Conversation analysis of transcribed and non-transcribed enhancers by PhastCons. g, h Mean count of conserved noncoding sequences (CNSs) (g) and single-nucleotide polymorphisms (SNPs) (h) per one hundred transcribed or non-transcribed enhancers. Significant differences among groups were analyzed using the one-tailed Student’s t test

Recent studies in plants have shown that open chromatin regions are usually accompanied by low levels of DNA methylation. We detected methylation within 1 kb on both sides of the enhancer midpoints. The DNA methylation level of central enhancer regions was approximately 0.01 (Fig. 2b), which was significantly lower than the average methylation level of the whole genome of Arabidopsis [44]. DNA methylation around transcribed enhancers was lower and was more negatively correlated with chromatin accessibility than that around nontranscribed enhancers (Fig. 2c). The plot of methylation and chromatin accessibility levels across 5 chromosomes showed that enhancers, especially transcribed enhancers, were enriched in euchromatin and depleted in heterochromatin (Additional file 2: Fig. S2b). Large constitutively hypomethylated regions in the genome are usually referred to as DNA methylation valleys (DMVs) or unmethylated regions (UMRs), which have unique chromatin characteristics and may contain functional genes and regulatory elements [20, 45]. Therefore, we identified a total of 71221 UMRs in the Arabidopsis genome, among which 15.43% were located in intergenic regions, 42.39% were located in promoter regions, and 42.19% were located in gene body regions (Fig. 2d). To determine whether the identified UMRs and enhancers represent the same regulatory elements, we analyzed the overlap between enhancers and UMRs. Not surprisingly, approximately 99% of the transcribed enhancers were located in UMRs, while nontranscribed enhancers did not show such significant overlap (Fig. 2d). Thus, these results revealed that transcribed enhancers were more likely to be hypomethylated.

The evaluation of chromatin characteristics related to enhancer transcriptional activity impelled us to explore whether enhancer sequences exhibit uneven evolutionary conservation. Hence, we quantified the sequence conservation of enhancers in Arabidopsis and other cruciferous species based on the PhastCons conservation score [46]. Enhancer regions were shown to be more conserved than their flanking regions, and the conservation score of transcribed enhancers was higher than that of nontranscribed enhancers (Fig. 2f). Furthermore, evolutionarily conserved noncoding sequences (CNSs) of Arabidopsis [46, 47] showed higher enrichment within transcribed enhancers (Fig. 2g), which further suggested that transcribed enhancers are more evolutionarily conserved. EnhancerA3732, which presented the highest induction fold among the enhancers whose transcription was upregulated by flg22 treatment, was highly conserved among four related species (Additional file 2: Fig. S2e). In addition, the single-nucleotide polymorphisms (SNPs) significantly correlated with lesion size identified within the Arabidopsis/Botrytis cinerea pathosystem by a genome-wide association study (GWAS) [48] were more enriched in transcribed enhancers (Fig. 2e). Previous studies using transgenesis assays focusing human and zebrafish developmental enhancers or human and mouse heart enhancers showed a high degree of functional conservation despite sequence divergence [49]. We mapped read coverage across the Arabidopsis genome based on the data described above (Supplementary Fig. 2d). Taken together, the results showed that the transcribed enhancers displayed more potential features related to conservation in evolution and function in plant disease resistance.

Immune TF binding sites are enriched within PTI-induced enhancers

To understand the potential mechanisms of enhancers regulated by flg22, we investigated the enrichment of TF binding motifs in transcribed enhancers relative to that in nontranscribed enhancers. Intriguingly, the majority of the differentially enriched TF motifs (P-value < 0.05) in upregulated enhancers were WRKY TF binding sites (Fig. 3a, Additional file 2: Fig. S3a), which are mostly involved in the regulation of plant defense-associated gene expression [50]. In contrast, the differentially enriched TF motifs located in downregulated enhancers were the binding sites of diverse TF families. We further analyzed the motifs enriched in upregulated enhancers using downregulated enhancers as controls. Similarly, the enrichment of WRKY family TF motifs in upregulated enhancers was apparent (Fig. 3b, Additional file 2: Fig. S3b). Additionally, the W-box WRKY TF binding motif, a DNA motif with the core sequence TTGAC(T/C), appeared more frequently in upregulated enhancers (Fig. 3c). Most upregulated enhancer-related TFs exhibited significant upregulation post-flg22 treatment (Additional file 2: Fig. S3c), indicating that these TFs may be implicated in plant innate immunity.

Fig. 3

WRKY family transcription factors are enriched in up-transcribed enhancers. a Significant enrichment (P-value < 0.05) of transcription factors binding motifs in transcribed enhancersA relative to that in non-transcribed enhancersA. b Significant enrichment (P-value < 0.05) of TF binding motifs in upregulated enhancersA relative to that in downregulated enhancersA. The percentage of sequences in the target group versus the background group are displayed to the left of the genes. Enrichment P-values are listed to the right of the genes as −log10 transformed values. c, e–g Average number of W-box motif (TTGACC/T) (c) and ChIP-seq peaks of transcription factors (e–g) in up- or downregulated enhancers, and non-transcribed enhancers. WRKYs binding regions with flg22 treatment for 2 h (e); other differentially enriched TFs binding regions (ASL18, mock; LBD18, mock; NLP4, mock) (f); and other immune TFs binding regions (HD2B, 30 min after flg22 treatment; IDD4, 1 h after flg22 treatment; SARD1, 24 h after Pseudomonas syringae pv. maculicola (P.s.m.) ES4326 treatment) (g). Different letters denote significant differences by the one-way ANOVA test (P-value < 0.05). d Distribution of distances of peaks for ChIP-seq using the anti-all-WRKY antibody at control condition to the nearest enhancer midpoint

To further confirm the enrichment of WRKY TF binding sites, we integrated the chromatin immunoprecipitation-sequencing (ChIP-seq) results of WRKY TFs obtained with/without flg22 treatment [51]. Relative to downregulated enhancers and nontranscribed enhancers, upregulated enhancers were located closer to the binding peaks identified following ChIP-seq performed with an anti-all-WRKY antibody under normal conditions (Fig. 3d). Upon flg22 treatment, the number of binding peaks for both all-WRKYs and immune-related WRKYs, such as WRKY18, WRKY33, and WRKY40, in upregulated enhancers was almost three times greater than those in other types of enhancers, while no significant difference could be observed between downregulated enhancers and nontranscribed enhancers (Fig. 3e). For instance, EnhancerA1671 and EnhancerA1286, which are located upstream of the WRKY18 and WRKY40 loci, respectively, showed a dramatic increase in WRKY binding after flg22 treatment (Additional file 2: Fig. S3e). These results suggested that multiple WRKY TFs, inducing the key immune regulators WRKY18, WRKY33, and WRKY40, could preferentially bind to upregulated enhancers, and their binding activity was further promoted by flg22-mediated immune signaling.

We continued to analyze the binding of other differentially enriched TFs to up- and downregulated enhancers. Surprisingly, we found a very low number of binding regions for ASYMMETRIC LEAVES2-LIKE 18 (ASL18), LOB DOMAIN-CONTAINING PROTEIN 18 (LBD18), (NODULE INCEPTION)-LIKE PROTEIN 4 (NLP4), TANDEM ZINC FINGER PROTEIN 9 (TZF9), BASIC HELIX LOOP HELIX PROTEIN (bHLH64), and AT2G41835 TFs in the identified enhancers, and there was no significant difference between those located in up- or downregulated enhancers (Fig. 3f, Additional file 2: Fig. S3d). The ChIP-seq results obtained for HISTONE DEACETYLASE 2B (HD2B), INDETERMINATE-DOMAIN 4 (IDD4) and SARD1 during immune response progression were also included in the analysis [38, 52, 53]. Unlike HD2B and IDD4, the binding peaks of SARD1, an activator of plant immunity that promotes the production of SA and the activation of defense gene expression, were significantly enriched in upregulated enhancers relative to downregulated and nontranscribed enhancers (Fig. 3g). This indicated that upregulated enhancers may recruit specific WRKYs and other immune-related TFs to activate their own transcription and subsequently regulate target immune gene expression.

Interactome of enhancers and immune genes with PTI-induced transcription

To dissect the regulatory network between immune-related enhancers and gene expression, we chose different methods to construct the comprehensiveness and accuracy of the interactome of enhancers and genes showing PTI-induced transcription. We first connected enhancers to genes by choosing the closest annotated protein-coding gene to each enhancer (Additional file 5: Table S4). Gene Ontology analysis of the target genes of upregulated enhancers showed that they are mainly linked to the responses to organic substances, stimuli and stresses, regulation of defense responses, innate immune system, and negative regulation of cell death (Fig. 4a). In addition, downregulated enhancer-targeted genes were involved in the response to a hormone stimulus, regulation of biological processes, biosynthetic processes, and metabolic processes (Fig. 4a). To obtain the informative features of target gene transcription dynamics, we referenced time-series RNA-seq data obtained from Arabidopsis leaves in response to flg22 treatment [54], and the expression dynamics of target genes were illustrated in a density map. The closest genes to upregulated enhancers showed significant upregulation in the early stage (1, 2, 3, and 5 h) and presented gradually decreased expression at the middle and late stages (9 and 18 h), while more downregulated enhancer-targeted genes showed decreased expression at all time points (Fig. 4b). The functional analysis and gene expression results indicated that upregulated enhancers play important as-yet-undiscovered roles in the regulation of immune gene expression.

Fig. 4

Induction of immune responsive genes are related with the upregulated transcribed enhancers. a Enrichment of the closest genes of upregulated enhancers (left panel) and of downregulated enhancers (right panel) with GO terms. b Frequency distribution of difference in gene expression, which were from the closest genes of upregulated (left panel) and of downregulated enhancers (right panel), post flg22 treatment compared with mock water treatment. c Representative example of immune genes interacted with upregulated enhancers in the loops. Pale blue bars indicate eRNA transcription regions. Plus and minus signs indicate Arabidopsis seedlings treated with flg22 and those that were mock-treated, respectively. The raw read counts of each gene are shown. Frw, forward strand; rev, reverse strand. d Integrated gene regulatory networks (iGRNs) among upregulated enhancersA, WRKY transcription factors and regulatory genes

To facilitate gene regulation, distal enhancer elements interact physically with promoter elements via chromatin looping at the TSSs of their target genes [55]. We further studied the spatial relationship between enhancers and target gene expression using statistically significant INT-Hi-C (combining the isolation of nuclei tagged in specific cell types and Hi-C) interaction data with a 2-kb resolution from Arabidopsis leaves [56]. The results showed that enhancersA and enhancersD interacted with 2069 and 3033 genes, respectively, and there were more interaction pairs between transcribed enhancers and genes (Additional file 2: Fig. S4a). Upregulated enhancers were associated with more flg22-treated upregulated genes than downregulated enhancers at 1 h. Some of the genes associated with upregulated enhancers have been proven to play important roles in plant immunity, including important immune genes encoding WRKY15, WRKY48, MYB70, and ETHYLENE RESPONSIVE ELEMENT BINDING FACTOR 4 (ERF4) TFs; immune signaling-related genes such as MITOGEN-ACTIVATED PROTEIN KINASE KINASE KINASE 2 (MAPKKK2), JASMONATE-ZIM DOMAIN 9 (JAZ9), AT3G46710 (an NLR gene); and secondary metabolism-related genes such as CYP71B22 and CYP71B23 (Fig. 4b, Supplementary Fig. 4b). We also found that some de novo-upregulated enhancers interacted with the target upregulated genes at 0 h, but no enhancer transcripts were detected at this time point, suggesting that the transcription of these enhancers may be induced in a signal-dependent manner. Taken together, the results of the identification of enhancer-gene interactional loops provide an overview of immune gene regulation at the 3D genome level.

Considering the significant binding of WRKYs to upregulated enhancers, we constructed potential integrated gene regulatory networks (iGRNs) based on the assessment of the coexpression intensity between upregulated enhancers carrying W-box motifs with WRKY genes and other regulatory genes to enhance the global understanding of upregulatory interactions. In total, the iGRNs covered 58 enhancers and 252 target genes (954 interactions) (Fig. 4c, Additional file 2: Fig. S4c). Among these sequences, 24 WRKY genes that may play a regulatory role were identified and specifically labeled. In the iGRNs, most enhancers are associated with more than five genes. EnhancerD3701, on chromosome 3, showed the highest frequency of gene association. There were 32 genes coexpressed with this enhancer, including 14 WRKYs. Similarly, some genes were associated with multiple enhancers. Based on the networks, the enhancer-TF gene connections were further assessed. Thus, we preliminarily established a causal link between flg22-induced enhancer activation and corresponding immune gene activity on a genome-wide scale based on the perspective of the nearest neighbor strategy, physical interaction, and coexpression.

Comparison of different immune elicitors regulated enhancers

To further reveal the importance of enhancers in plant innate immunity, we analyzed the enhancers regulated by several other patterns from different source organisms, including chitin (an oligosaccharide fragment released from fungal cell walls), nlp20 (a 20-amino-acid fragment of NECROSIS AND ETHYLENE-INDUCING PEPTIDE 1-LIKE PROTEINS produced by oomycetes and bacterial and fungal microbes), and pep2 (a 23-amino-acid DAMP peptide from Arabidopsis endogenous peptide 2). The analysis also included INA, which is a functional analog of salicylic acid (SA) and can amplify defense response signals as a secondary signal molecule [57, 58]. Using the same workflow and pipeline employed for flg22 data analysis, we identified a similar number of transcribed enhancers during immune responses downstream of these elicitors (Additional file 6: Table S5). In contrast to flg22 elicitation, the numbers of genes and enhancers whose upregulation was induced by these four elicitors were relatively lower. A total of 841, 920, and 948 upregulated enhancers were identified at 1 h after chitin, nlp20, or pep2 treatment, respectively, while more enhancers were detected after flg22 treatment (Fig. 5a). The plant defense hormone INA induced the transcription of 833 enhancers (Fig. 5a, Additional file 7: Table S6). Importantly, the upregulated transcribed enhancers induced by these elicitors similarly showed the enrichment of immune-related SARD1 and WRKY TFs (Fig. 5b), reflecting the convergence between early immune signaling triggered by different patterns.

Fig. 5

Regulated transcribed enhancers show common characteristic and specificity during multiple PTI signaling. a Volcano Plots shown significant differential expressed transcribed enhancers identified upon the treatments of different immune elicitors. The number of upregulated enhancers is shown. b Average number of WRKYs and SARD1 binding regions in upregulated enhancers (bright color) and downregulated enhancers (light color). c Venn plots showing the overlap number of upregulated enhancers (left panel) and upregulated genes (right panel) that identified upon different patterns treatments. d Box plots of sum of fold change in expression (log2 scale) of core pattern-induced enhancers (CPIEs) (upper panel) and core pattern-induced genes (CPIGs) (lower panel), each pattern treatment in long time series (5, 10, 30, 90, and 180 min) in Col-0 wild-type (WT) and cognate receptor mutant. Note that wak1 mutants are not viable, and thus the OG treatment was paired with a mock water treatment. e Heatmaps of fold change in expression (log2 scale) of 4 patterns co-regulating upregulated enhancers (left) and upregulated genes (right) in each elicitor treatments

To further investigate the similarity of transcribed enhancers in different immune events, the commonly upregulated enhancers and upregulated genes induced by all these patterns were screened. There were 129 enhancers that could be upregulated by all tested patterns, which were defined as core pattern-induced enhancers (CPIEs), while different patterns commonly induced the expression of 90 genes that were considered core pattern-induced genes (CPIGs) (Fig. 5c, Additional file 8: Table S7). To demonstrate the core enhancers and genes included in immune signaling pathways triggered by different patterns, we cited gene expression data obtained after stimulation with different several patterns reported in recent publications [38, 59]. The expression of CPIE-linked genes and CPIGs was clearly upregulated after seven different pattern treatments in wild-type plants, but their induction was abolished in the corresponding receptor mutant (Fig. 5d). Gene annotation of the closest CPIE genes revealed that many of them are important plant immune-related genes, such as ERF4 and LecRK-IX.2 (Additional file 8: Table S7). The corresponding CPIEs also showed evolutionary conservation in five cruciferous species (Additional file 2: Fig. S5a). We found that the functions of CPIE-related genes were enriched in the categories of response to stimulus, response to stress, and immune response (Additional file 2: Fig. S5b). The heatmaps showed the expression patterns of upregulated enhancers and genes induced by different patterns (Fig.