To gain deeper insight into the GmLUX (Glyma.11G136600) regulatory network in soybean, two replicates of GmLUX ChIP-seq were generated from leaves of GmLUX-flag transgenic plants to characterize the binding preference genome-wide ( Figure 1A ; Supplementary Figure S1 ). Up to 55 M reads were obtained for the ChIP-seq and input libraries ( Supplementary Table S1 ). After alignment to the soybean reference genome ( Figure 1A ), a promoter and downstream enrichment pattern of GmLUX was observed ( Figure 1B ), and a total of 4,373 GmLUX binding peaks were identified by the MACS2 software ( Figure 1C ; Supplementary Table S2 ). Additionally, these peaks tended to be located within 2 kb from the transcription start site (TSS) of their corresponding nearest genes, such as the promoter and exonic regions ( Figures 1D, E ). Furthermore, the gene ontology (GO) enrichment analysis of the GmLUX peak-associated genes (3,312) ( Supplementary Table S2 ) revealed several enriched GO terms, such as nucleic acid binding transcription factor activity (GO:0001071), transcription factor activity (GO:0003700) and sequence-specific DNA binding (GO:0043565), suggesting a robust transitional regulatory role of GmLUX ( Figure 1F ; Supplementary Table S3 ).

To further narrow down the potential genes directly targeted by GmLUX and constructed a GmLUX regulatory network, a comprehensive integration of GmLUX ChIP-seq, co-expression and motif scanning data was conducted, resulting in the identification of 33 key genes ( Figures 2A, B ; Supplementary Table S4 ). These genes, harboring a GmLUX binding motif in their promoters, were found to be co-expressed with GmLUX and associated with GmLUX enriched peaks, such as Glyma.20G237200 and Glyma.17G230500 ( Figures 2C, D ). Subsequently, these genes were then used to construct a putative GmLUX regulatory network with GmLUX as the central hub ( Figure 2B ). In addition, flower related genes ( Supplementary Table S2 ) were identified in this network, such as DOF (Glyma.13G062500), AP2/ERF (Glyma.02G261700), CO-like (Glyma.10G274300) and DnaJ (Glyma.08G220000) ( Figure 2B ), which were reported to have potential roles in regulating flower and served as potential candidate genes for further investigation.

The epigenetic states of chromatin are usually closely associated with TF binding sites, serving as an additional layer of gene regulation (Cuellar-Partida et al., 2012). Hence, three histone modifications were generated by ChIP-seq, including H3K9ac, H3K27ac and H3K4me1; moreover, ChIP-seq data for four additional histone modifications (H3K27me3, H3K4me3, H4K12ac and H3K14ac) were downloaded, along with ATAC-seq data, to comprehensively understand the epigenetic states around GmLUX binding peaks ( Figure 3A ; Supplementary Figures S2, S3 ). As shown in Figure 3A , a strong enrichment of ATAC-seq, H3K4me3 and H4K12ac was observed at the center of GmLUX binding peaks, accompanied by a relatively weaker signal of H3K27ac, H3K14ac and H3K9ac, while little enrichment of H3K27me3 and H3K4me1 signals was observed at these regions. Consistent with these observations, we identified over 2,000 GmLUX binding peaks with ATAC-seq, H3K4me3 and H4K12ac modifications, and a few of peaks with H3K27me3 and H3K4me1 modifications ( Figure 3B ). Since H3K27me3 modification is a repressive mark, in contrast to the remaining active epigenetic marks, a significant difference was observed between the expression of H3K27me3-modified GmLUX binding peaks associated genes and those genes with other epigenetic mark modified peaks ( Figure 3C ). These observations indicated that GmLUX binding peaks associated with the active epigenetic marks (e.g. ATAC, H3K4me3 and H4K12ac) ( Supplementary Figure S5 ) to regulate downstream gene expression.

Furthermore, many GmLUX binding peaks with multiple modifications were observed, which could have a significant impact on the downstream gene expression ( Figure 3D ; Supplementary Table S5 ). For example, gene expression with ATAC and acetylation (Kac) modified peaks was significantly higher than that with ATAC alone, and similar situations were also found in the H3K4me3&Kac versus H3K4me3 group, as well as Kac versus H3K12ac group ( Figure 3E ).

In addition, the epigenetic modification enrichment pattern of the distal GmLUX binding peaks was distinct from the proximal peaks, except for the ATAC-seq signal ( Supplementary Figure S4 ). As the distal ATAC-seq signal was served as an indicator of potential distal enhancer-like elements, this observation suggested an enhancer role for distal GmLUX binding peaks, which required further investigation.

As a palaeopolyploid plant (Schmutz et al., 2010), soybean contains lots of duplicated genes in its genome. As previously reported, the gain and loss of cis-regulatory sequences (e.g. the ACR) could have a significant effect on the expression of homologous genes (hGenes) (Huang et al., 2021c; Fang et al., 2023). We focused on 1,011 genes with a promoter GmLUX binding site, identified their corresponding hGenes, and checked whether their hGenes shared the GmLUX binding peaks with them. Among these 1,011 paired hGenes, only 224 hGenes (22.2%) were bound by the GmLUX, while 787 hGenes (77.8%) exhibited an absence of GmLUX binding ( Supplementary Table S6 ). Furthermore, it was observed that those hGenes with both GmLUX binding peaks had similar epigenetic modifications and showed no significant impact on gene expression (p = 0.673) ( Figures 4A, B ; Supplementary Table S6 ). However, the expression of those hGenes without GmLUX binding peaks was significantly lower than their hGenes ( Figures 4A, B ; Supplementary Table S6 ). In addition, the absence of GmLUX binding peaks also caused obvious changes in the epigenetic features of hGenes ( Figure 4B ; Supplementary Table S6 ). These data suggested that the presence or absence of cis-regulatory sequences during gene duplication were important for maintaining the expression of hGenes, which may have implications for the subsequent functionalization of these genes.

Enhancer elements played an important role in gene activation and could be located in the intergenic or intronic accessible regions (Zhu et al., 2015; Meng et al., 2021). In addition, as previously reported, TF could mediate the promoter-enhancer interaction to modulate gene expression in soybean (Huang et al., 2023). In this study, a total of 1,018 and 301 GmLUX binding peaks were observed located in the intergenic and intronic regions of the soybean genome, respectively ( Supplementary Table S2 ). Of these, 837 intergenic and 160 intronic peaks overlapped with ATAC-seq peaks. Furthermore, a total of 114 genes were associated with multiple GmLUX binding peaks (41.8% associated with promoter and intronic GmLUX peaks and 58.2% associated with promoter and intergenic peaks) ( Supplementary Table S2 ). Since ACRs identified by ATAC-seq served as an indicator of enhancer-like elements, this observation indicated that GmLUX was able to mediate enhancer-like elements to regulate gene expression.

To further test our hypothesis, a flowering-related gene, CO-like (Glyma.10G274300) was selected as a candidate for further validation. This CO-like gene harbored multiple cis-regulatory regions located in its upstream region (three peaks of ATAC-seq and two peaks of GmLUX) ( Figure 5A ). Interestingly, CO-like was associated with several active histone modifications (e.g. H3K4me3, H4K12ac, H3K14ac etc.) and had relatively high expression in soybean ( Figure 5 ). Moreover, the upstream region of this CO-like gene located in a topologically associating domain (TAD, showed as a black triangle) region according to the Hi-C data, indicating a high frequency of chromatin interaction in these regions ( Figure 5A ). To confirm such an interaction, the 3C-PCR was performed using three primers (P1, P2 and P3) and an anchor primer (AP). It showed that the interaction signals of AP-P1 and AP-P2 were detected, except for that of AP-P3 ( Figure 5B ). Since AP and P1 were close to two GmLUX binding peaks, the data suggested that GmLUX mediated the AP-P1 interaction to regulate CO-like gene expression. Interestingly, in the transient experiments, only the P1 but not the AP cis-regulatory sequences had relatively high activities in protoplasts and could be activated by GmLUX overexpression ( Figure 5C ), indicating the distal binding region (e.g. P1) was necessary for GmLUX to regulate its target genes (e.g. the CO-like gene).

Glycine max (L.) Merr. cultivar Williams 82 (W82) was used as the wild type in this study, and transgenic W82 with 35S:GmLUX-Flag (Bu et al., 2021). Transgenic and wild-type seeds were sterilized with 75% ethanol and placed on sterilized vermiculite for germination, until their cotyledons were fully expanded. These seedlings were transferred to soil and grown in the greenhouse (16 h light/8 h dark, 25-28°C) to develop the healthy trifoliate leaves, which were utilized for subsequent experiments.

The construction of GmLUX ChIP-seq library was followed the previous study (Huang et al., 2021b), with minor modifications. Briefly, about 0.5 g of leaves from GmLUX-flag transgenic plants were utilized for nuclear isolation. The chromatin fragments were pull-down by the anti-Flag antibody (Sigma). The ChIPed and input DNA used for library construction were processed by Tn5 transposase (Vazyme). Histone ChIP-seq library construction was followed the method (Huang et al., 2021b; Zhang et al., 2023) and the anti-H3K27ac, anti-H3K9ac and anti-H3K4me1 antibodies (Millipore) were utilized in this study. All the ChIP-seq libraries were sequenced on the PE150 mode of the Illumina platform.

The raw data of ChIP-seq were trimmed by TrimGalore (http://github.com/FelixKrueger/TrimGalore) and the filtered reads were mapped to the soybean reference genome (V4, https://phytozome.jgi.doe.gov) by the Bowtie2 (Langmead and Salzberg, 2012). The high-quality mapped reads (MAPQ > 30) were extracted by SAMtools and used for peak calling by MACS2 (Zhang et al., 2008; Li et al., 2009). For broad peak calling, the parameters were set as ‘-trackline -extsize 147 -broad -q 0.01 -nomodel -buffer-size 500000’, while the default parameters were used for narrow peak calling. The input library was used as a control. The profile of ChIP-seq is normalized to input controls by the BamCompare function in Deeptools (Ramírez et al., 2014). Enriched peaks annotation was performed via HOMER (http://homer.ucsd.edu/homer/). The visualization of ChIP-seq data was processed using the Deeptools (Ramírez et al., 2014), pyGenomtrack (Lopez-Delisle et al., 2021) or IGV (Thorvaldsdóttir et al., 2013) software.

The RNA-seq data were obtained from a previous study (Huang et al., 2021b). Raw RNA-seq data were filtered by TrimGalore with the default parameters, and then mapped to the soybean reference genome using HISAT2 (http://daehwankimlab.github.io/hisat2). Gene expression levels were calculated by Cuffnorm (http://cole-trapnell-lab.github.io/cufflinks/cuffnorm) and represented by the fragments per kilobase per million mapped reads (FPKMs), which were utilized for subsequent analysis.

The leaf Hi-C data were downloaded from a previous report (Wang et al., 2021) and followed its analysis pipeline with minor modification. Briefly, the raw data was trimmed with the adaptors and were processed by the HiC-Pro (Servant et al., 2015) with default parameters. Seft-circle, dangling-end reads were removed, and the contacts were normalized by the iterative correction and eigenvector decomposition (ICE) method. Then the TAD domain was called by HiCexplore (Wolff et al., 2022) with default parameters and visualized by pyGenomtrack (Lopez-Delisle et al., 2021).

The GmLUX motif (Motif ID: TFmatrixID_0354) position weight matrix (PWM) was downloaded from the PlantPAN database (Chow et al., 2019) and the FIMO from the MEME suite (Bailey et al., 2009) was used for motif scanning in the GmLUX peak region.

The soybean leaf ATAC-seq raw data were obtained from previous study (Huang et al., 2021b), and filtered by the TrimGalore and mapped to the reference genome via Bowtie2 (Langmead and Salzberg, 2012) with the parameters set as: –very-sensitive -X 1000. The mapped reads with a MAPQ value over 30 were used for peak calling via Genrich software on the ATAC mode (https://github.com/jsh58/Genrich). Visualization and annotation of ATAC-seq data was performed similarly to the ChIP-seq mentioned above.

To confirm the transcriptional activation of GmLUX, the full-length CDS of GmLUX was cloned into the pRT107 vector which was derived by the 35S promoter, while the two GmLUX enriched peak regions associated with the CO-like gene were cloned into the pGreenII-0800-LUC vectors and confirmed by Sanger sequencing. Then the validated plasmids of pRT107-GmLUX and pGreenII-0800-GmCO-like-Luc were co-transformed into Arabidopsis protoplasts, isolated from the leaves of 3-week-old Col-0 by the PEG/CaCl2 method (Yoo et al., 2007). After 16 h of incubation at 22°C, the total proteins were extracted, and transcriptional activities were measured using the Dual-Luciferase Reporter Assay System (Promega, E1910) kit followed the manual instructions. The fluorescence value of firefly luciferase (LUC) and Renilla luciferase (Ren) was determined by Promega GloMax 20/20 system. The transcriptional activity of GmLUX against the CO-like associated peaks was evaluated by the ratio of LUC/Ren value. The experiments were performed in three biological replicates.

The 3C-qPCR was conducted according to the previous study (Hagege et al., 2007; Liu, 2017; Huang et al., 2023). Briefly, the isolated soybean leaf nuclei were digested with DnpII (NEB), and the sticky ends were filled with dWTP (W: A, T, G), biotin-14-dCTP (Sigma) and 40 U Klenow fragment (NEB). The blunt ends were then ligated with the T4 DNA ligase (NEB) following the manufacturer’s protocol. The biotin-labeled 3C-DNA was extracted by MyOne™ Streptavidin C1 Dynabeads (Invitrogen) according to the manufacturer’s protocol and used for qPCR using the TB Green Premix Ex II (Takara). The primers involved in 3C-qPCR analysis were listed in Supplementary Table S1 .