Soybean LBD genes (GmLBDs) have been previously identified (Yang et al., 2017). To further understand the functions of GmLBDs in their origination, classification, and even the evolutionary relationship with other species in Leguminosae, we first identified LBD family members in the Phaseolus vulgaris genome (P. vulgaris, common bean), a species with a relatively close evolutionary relationship with soybean in Leguminosae, and in the Cucumis sativus genome (C. sativus, cucumber), an eudicot species. A local BLASTP search was carried out using 42 known Arabidopsis LBD proteins as the query in common bean and cucumber genomes in the NCBI database. Subsequently, all potential LBD protein sequences were further verified by domain analysis using Pfam (LOB domain, DUF260, Pfam number: Pfam03195). As a result, a total of 42 CsLBDs in cucumber and 50 PvLBDs in common bean were finally identified. CsLBDs and PvLBDs were named according to the order of locations on the chromosomes (Supplementary Tables 1, 2).

To further improve the understanding of the phylogenetic classification and evolution of the LBD family in the soybean genome, a comprehensive phylogenetic tree was constructed using 788 amino acid sequences of LBD protein from 16 plant species, including one species each of green alga, moss, fern, and basal angiosperm; eight species in eudicots; and four species in monocots (Supplementary Table 3). The 788 amino acid sequences of LBD protein contained 696 known LBD proteins from 14 species and 92 LBD proteins that were newly identified in this study (Yang et al., 2017; Zhang et al., 2020).

Phylogenetic analysis showed that all LBD proteins were classified into two classes (Class I and Class II); Class I is further divided into five subgroups, namely, Class IA, Class IB, Class IC, Class ID, and Class IE, whereas Class II is divided into Class IIA and Class IIB (Figure 1), which is consistent with the previous results (Tang, 2013; Chanderbali et al., 2015; Zhang et al., 2020). Among the 90 GmLBDs, 74 were clustered into Class I, with 19 in Class IA, 25 in Class IB, 19 in Class IC, 4 in Class ID, and 7 in Class IE, while 16 GmLBDs were clustered into Class II.

Given that soybean and common bean have been demonstrated with a close genetic relationship (Vlasova et al., 2016) and both of them are important cash crops. We, therefore, constructed a phylogenetic tree between soybean and common bean to explore the evolutionary relationship of LBD genes in these two genomes using full-length protein sequence. Phylogenetic analysis showed that the homologous relationships between GmLBDs and PvLBDs were obviously observed since almost all clades were included by one PvLBDs and one or two GmLBDs (Figure 2). The homologous relationships were inspected by checking GmLBDs and PvLBDs in the same clades and summarized (refer to Table 1). Intriguingly, we found that a total of 38 homologous gene groups were detected, including all PvLBDs and 91% of GmLBDs (82/90), suggesting that the gene duplication in soybean LBDs is another character. In around 70% of clades (25 in 42 clades), one PvLBD and two GmLBDs were closely clustered into one new clade. In the homologous gene group summary, 25 (65%) groups showed the gene ratio of PvLBDs:GmLBDs as 1:2, i.e., one PvLBD has two GmLBDs orthologs.

Soybean experienced one-time independent whole-genome duplication (WGD) compared with common bean, and they diverged only 19.2 million years ago, a relatively short time compared with other legume sister species. To verify the mechanism of gene duplication between soybean and common bean LBD genes, syntenic maps of LBD homologs in these two genomes were built; GmLBD89 and GmLBD90 were excluded in syntenic analysis because of the unassembled genome locations (Figure 3). As a result of 90 GmLBDs and 50 PvLBDs, 112 collinear gene pairs were detected and merged into 38 collinear groups. The collinear groups contain 93% (82/88) of GmLBDs, and 100% (50/50) of PvLBDs genes, which is perfectly consistent with the homologous gene groups summary (Table 1). Interestingly, all these collinear gene pairs contained about 1.62 pairs of conserved genes on average, which is also consistent with the summary of paralog ratios between these two species (Table 1). The syntenic analysis results again proved the close evolutionary relationship of these two species and suggested that WGD might be the main force for LBDs expansion in the soybean genome.

A total of 18 AtLBD ancestries in angiosperms have been previously retraced based on the gene collinearity investigation and phylogenetic relationships (Kong et al., 2017). To further characterize the evolutionary patterns, we also tried to trace the ancestries in GmLBDs and PvLBDs. Interestingly, our data exhibited that GmLBDs and PvLBDs can be detected in all 18 AtLBDs ancient lineages (Table 2), suggesting that no ancestor genes were lost in soybean and common bean genomes. In each AtLBD ancient lineage, 2–10 GmLBD and 1–6 PvLBD paralogs were presented. In most ancient lineage, the number of GmLBDs was much more than that of AtLBDs, such as in lineage 2, AtLOB and AtLBD25 vs 10 GmLBDs and in lineage 11, AtLBD3 vs 6 GmLBDs, which indicated the extensive expansion of GmLBDs in these ancient lineages. However, in the common bean genome, no obvious gene expansion was found except in lineages 2, 9, and 11. The decrease of gene number was also found in some ancient lineages for both GmLBDs and PvLBDs, such as in ancient lineage 8, 4 AtLBDs vs 2 GmLBDs vs 1 PvLBD and in ancient lineage 15, 6 AtLBDs vs 3 GmLBDs vs 2 PvLBD. The ancestry retracement in soybean and common bean demonstrated that LBD is reluctant to be lost and the similar expansion and decrease patterns in some ancient lineages between soybean and common bean mean that they might suffer from parallel evolution.

LBD proteins have been reported to play important roles in controlling plant growth and development and also in responding to stress, such as pathogen infection (Iwakawa et al., 2002; Shuai et al., 2002; Feng et al., 2012). To further examine the potential roles of GmLBD proteins in plant immunity, especially in response to P. sojae infection, some GmLBD candidate genes were first identified. Based on previous studies, there are some LBD proteins, which have previously been characterized to involve in plant immune response or upregulated by pathogen infection. AtLBD20 in Arabidopsis and CsLOB1 in Citrus sinensis were found to be involved in plant immunity response to the pathogen (Thatcher et al., 2012b; Hu et al., 2014). So a BLASTP search was performed against the soybean genome database¹ using AtLBD20 and CsLOB1 as queries and identified 6 GmLBD homologs. In total, 13 GmLBDs were reported to be highly induced in responses to biotic stresses (Yang et al., 2017). Accordingly, a total of 19 GmLBDs genes were selected as candidates for further functional characterization (Supplementary Table 4).

To determine whether these candidate genes play roles in plant defense response against P. sojae attack, the expression patterns of 19 GmLBDs genes upon P. sojae infection were examined. The quantitative reverse-transcription PCR (qRT-PCR) was performed using RNA that was extracted from hairy roots of soybean susceptible species Huachun 6 and collected at different time points [0, 1.5, 3, 6, 12, 20, and 24 h after infection (hpi)] after P. sojae infection. The results showed that 16 out of 19 genes were successfully amplified by qRT-PCR, and 15 GmLBDs were found to be induced in the early infection period except for GmLBD37 when compared with the uninfected samples (Figure 4). Among them, 9 GmLBD genes (GmLBD9, GmLBD16, GmLBD23, GmLBD88, GmLBD30, GmLBD55, GmLBD90, GmLBD43, and GmLBD70) were considered as highly upregulated genes, since the highest expression levels were increased at least 20-fold, especially for GmLBD90, which was induced more than 50-fold from 6 to 24 hpi and reached the highest expression level at 20 hpi for up to 150-fold (Figure 4A). The remaining six GmLBD genes (GmLBD31, GmLBD45, GmLBD51, GmLBD49, GmLBD59, and GmLBD63) were upregulated approximately 2–6 times (Supplementary Figure 1).

To further confirm the expression profiles of GmLBDs upon P. sojae infection, four GmLBD genes (GmLBD9, GmLBD16, GmLBD23, and GmLBD88), which showed highly induced expression patterns in Huachun 6, were selected to investigate their expressions in other soybean cultivars. We chose two soybean cultivars (Williams 82 and Tianlong 1) for further determining the expression pattern of these four genes. Williams 82, with a resistance gene (Rps1k) (Mideros et al., 2007), is known to be resistant to P. sojae strain P6497, and Tianlong 1 showed moderate resistance against P6497 than susceptible species Huachun 6 (personal communication). Compared to that in Huachun 6, GmLBD88, GmLBD16, and GmLBD23 genes also showed the upregulated expression patterns in Williams 82 and Tianlong 1 upon P6497 inoculation, while no clear expression change was detected in GmLBD9 (Figures 4B,C).

To further elucidate the function of these genes in soybean, their tissue-specific expression patterns were also examined in soybean roots, stems, and leaves by qRT-PCR. The results showed that these four genes were ubiquitously expressed in all plant organs tested, with the highest expression level in roots for GmLBD9, GmLBD16, and GmLBD23 genes and the highest expression level in stems for the GmLBD88 gene (Supplementary Figure 2). Given that most of the candidate GmLBD genes showed induced expression upon P. sojae infection, we inferred that they might play important roles in early soybean defense response.

To investigate the functions of GmLBD genes in soybean immunity, four GmLBD genes (GmLBD9, GmLBD16, GmLBD23, and GmLBD88) showed remarkably upregulated expression in Huachun 6 were chosen to reveal how they regulate soybean immunity through the transient overexpression and knockdown assays. GmLBD9, GmLBD16, GmLBD23, and GmLBD88 were first transiently overexpressed in soybean hairy roots by Agrobacterium-mediated transformation, and then those transformed hairy roots were inoculated with P. sojae strain P6497-RFP. We discovered that more P. sojae oospores can be observed in hairy roots expressing GmLBD9 and GmLBD88 at 48 hpi in relative to those hairy roots expressing empty vector (EV) by microscope observation (Figure 5A). Consistently, oospores and biomass of P. sojae were much higher in soybean hairy roots inoculated with GmLBD9 and GmLBD88 than in roots inoculated with the EV control (Figures 5C,D); these data indicate that the expression of GmLBD9 and GmLBD88 could promote the colonization of P. sojae in soybean hairy roots. Interestingly, we found that the individual expression of GmLBD16 and GmLBD23 inhibited P. sojae infection in soybean hairy roots, showing fewer oospores and lower relative biomass of P. sojae in transiently expressing hairy roots (Figures 5A, 6C,D). Moreover, the expression of those four recombinant proteins was confirmed by western blot, respectively (Figures 5B, 6B). Overall, these results suggest that GmLBD9 and GmLBD88 may be the negative immune regulators for soybean resistance against P. sojae, while GmLBD16 and GmLBD23 may be the positive immune regulators.

These data prompted us to further verify their roles in soybean Phytophthora root and rot. Two GmLBDs (GmLBD9 and GmLBD23) were selected to perform transient silencing in soybean hairy roots by RNA interference (RNAi) technique. GmLBD9 and GmLBD23 have been proved to be negative and positive immune regulators in soybean immunity, respectively. To specifically silence GmLBD9 or GmLBD23, approximately 200–300 bp of 5′- or 3′-end UTR fragments were cloned into pK7GWIWG2D vector to generate an RNAi recombinant construct. And then, these constructs were introduced into soybean hairy roots by Agrobacterium-mediated transformation. Each RNAi construct successfully silenced these targets as shown by qRT-PCR analysis and revealed that the expression of GmLBD9 and GmLBD23 were obviously reduced by 70–80% in silencing hairy roots (Figure 7B). Subsequently, these silenced hairy roots were challenged with P. sojae strain P6497-GFP. The results displayed that GmLBD9-silenced hairy roots showed less oospores (Figures 7A,C) and lower relative biomass of P. sojae than those hairy roots induced by EV (Figure 7D). Together with the above data, our results further indicated that the GmLBD9 gene negatively regulates soybean immunity against P. sojae infection. As such, the infection on GmLBD23-silenced roots exhibited that much more oospores were observed through fluorescence microscope observation (Figure 7A). Quantification of oospore number and biomass also proved that silencing of GmLBD9 promoted colonization of P. sojae in soybean hairy roots (Figures 7C,D), indicating that the GmLBD9 gene positively manipulates soybean immunity against P. sojae infection. Taken together, our data suggest that GmLBD9 and GmLBD88 are two negative immune regulators and GmLBD16 and GmLBD23 are two positive regulators of plant immunity against P. sojae infection.

Soybean cotyledons (Huachun 6, Williams 82 and Tianlong 1) were grown in a greenhouse at 25°C with a 16/8 h (light/dark) photoperiod. P. sojae strains, namely, P6497, P6497-RFP, and P6497-GFP, were routinely maintained on a 10% vegetable (V8) juice medium at 25°C in darkness.

To obtain cucumber and common bean LBD protein sequences, all known 43 Arabidopsis LBD protein sequences were used as a query to perform BLASTP with an e-value of 1 × e^––10 against the protein sequences database of C. sativus and P. vulgaris (NCBI²). Redundant sequences that are partial or alternatively spliced sequences from the same locus were removed. Then conserved domain of LBDs (LOB domain, DUF260, Pfam number: Pfam03195) acquired from Pfam³ was used for a blast to identify CsLBDs and PvLBDs with DUF260 as a query (El-Gebali et al., 2019). Finally, each gene was named based on its location on the chromosome.

To construct a phylogenetic tree of LBDs in 16 species, 788 full-length LBD protein sequences were aligned using the multi-sequence alignment program ClustalW. A phylogenetic tree was constructed with ML (maximum-likelihood) method in MEGA X and 100 times of bootstrap replicates. The phylogenetic tree was further manipulated by the program Interactive Tree of Life (iTOL⁴) (Zhu et al., 2019).

The chromosomal length and locations of each GmLBD and PvLBD were retrieved from the soybean genome database in SoyBase⁵ and common bean genome database in NCBI. Advanced Circos program in TBtools (Chen et al., 2020a) was used for collinearity analyses.

For overexpression assay in soybean hairy roots, fragments containing full-length CDS sequences were amplified with gene-specific primers (Supplementary Table 5) and then ligated into vector PFGC5941 by homologous recombination (Vazyme, C112-02-AB), which adds a C terminal FLAG tag. For gene silencing assay, fragments derived from the 5′ or 3′ UTR regions with 200–300 bp in length were amplified and then cloned into pK7GWIWG2D vectors with Gateway technology (Thermo Fisher Scientific, 12538120).

Soybean seeds (Huachun 6) were sanitized with a mixture of 84 disinfectants and concentrated hydrochloric acid (96:4) and then grown on germination medium (Gm medium). After around 6 days of growth, cotyledons were removed from soybean seedings and cut a wound (around 0.3 cm in diameter and 0.2–0.3 cm in depth) close to the petiole with a sterile knife. Agrobacterium rhizogenes (A. rhizogenes) K599 cell suspensions were inoculated on the wound, and then cotyledons were continued growing on Murashige and Skoog medium (MS medium). Around 3–4 weeks later, soybean hairy roots overexpressing or silencing GmLBDs were observed by fluorescence microscopy and further confirmed by western blot or qRT-PCR. Selected overexpressing hairy roots and silencing hairy roots were infected with P6497-RFP and P6497-GFP in wet and dark conditions at 25°C for around 48 h, respectively. Around 3- to 5-day-old P. sojae hyphae grown in 10% of V8 medium were used for soybean hairy roots infection.

Total RNA was isolated from hairy roots using TRIzol Reagent (TaKaRa, 9109). Isolated RNA samples were quantified using a NanoDrop spectrophotometer (Thermo Fisher Scientific, NanoDrop One) and then treated with DNase I (Thermo Fisher Scientific, AM2222) to remove any residual DNA contamination. In total, 1 μg of DNA-free RNA samples were converted to cDNA using a cDNA synthesis kit (Vazyme Biotech, R212-02-AF). First-strand cDNA was synthesized using HiScript II 1st Strand cDNA Synthesis Kit (+gDNA wiper) with Oligo-(dT) 23VN (Vazyme Biotech, R212-02-AF).

For measuring the transcript level of GmLBD genes from soybean during P. sojae infection, 1- or 2-week-old soybean (Huachun 6, Williams 82 and Tianlong 1) secondary roots was infected with R6497 and collected at different time points after infection. For measuring the transcript level of GmLBDs in different tissues, leaves, roots, and stems were sampled from 3- to 4-week-old soybean plants. Total RNA was extracted and used as a template for reverse transcription. To determine gene silencing efficiency, RNA was extracted from 3- to 4-week-old hairy roots induced by K599 containing PK7GWIWG2D (II) and PK7GWIWG2D (II)-GmLBDs. In all the quantitative PCR (qPCR) reactions, gene-specific primers (Supplementary Table 5) were designed specifically in 5′- or 3′-end UTR of genes avoiding the fragments for gene silencing. Then, qPCR was carried out using Maxima SYBR Green/ROX qPCR Master Mix (Vazyme Biotech, Q711-02-AA) (Zhang et al., 2019; Zhu et al., 2020). GmCYP2 gene (Hu et al., 2009) was used as an internal reference. Three independent biological replicates were conducted for each treatment with similar results.

Cetyltrimethylammonium bromide (CTAB) method was used for genomic DNA extraction (Li et al., 2017). Briefly, 500 μl of CTAB buffer was added into crushed soybean hairy roots to lyse plant cells in a water bath at 60°C for 1 h, and an equal volume of chloroform-isoamyl alcohol (24:1) was added while vigorously shaking for 30 s. After centrifugation at 12,000 rpm at room temperature, the supernatant was transferred into a new Eppendorf tube, an equal volume of ice-cold 100% ethanol was added to precipitate DNA at –20°C for more than half an hour, and then 75% ethanol was used to wash DNA. After drying in the hood, 100 μl of sterilized H₂O was added to dissolve DNA at 55°C. The extracted DNA was also quantified by Nanodrop and then diluted into the same concentration for biomass assay. Primers of GmCYP2 (Hu et al., 2009) from soybean and PsActin from P. sojae were used for biomass assay (Shi et al., 2020; Supplementary Table 5).

Proteins from the sample lysate were fractionated by SDS-PAGE gel. The separated proteins were transferred from gels to PVDF blotting membrane (GE, A10203127) (pretreated with methanol for 15 s) using Transfer Buffer (BIO-RAD, Cat. #10026938). The membrane was then blocked using 5% non-fat dry milk dissolved in TBST (TBS with 0.1% Tween 20) (also called TBSTM) for 1 h at room temperature with 30–40 rpm shaking, followed by three washes with TBST, and then primary antibody anti-Flag (1:5,000; MBL, M185-3L) and secondary antibody goat anti-mouse antibody (1:5,000; MBL, Lot 366) were applied to the membranes for 1 h in order. Finally, the membrane was visualized using the Western Blotting Substrate kit (Thermo Fisher Scientific, 34580) by multifunctional fluorescent molecular imager (GE, Amersham Imager 600) at 780 and 800 nm excitation (Chen et al., 2020b).

The statistical analysis in biomass assay, gene expression profile, and oospore numbers quantification was performed using an unpaired T-test in Graphpad Prism5.