Genome-Wide Identification of LATERAL ORGAN BOUNDARIES DOMAIN (LBD) Transcription Factors and Screening of Salt Stress Candidates of Rosa rugosa Thunb

Simple Summary The transcription factor family LBD were well-known as regulator of plant development. Several recent studies indicated LBD genes response to abiotic stresses and abscisic acid. The salt tolerant rose (Rosa rugosa) distribute wildly in coastal saline sands are ideal materials of exploring the salt -responsive LBD genes. We found 41 RrLBDs from genome of wild rose and classified them into Classes I and II. Interestingly, many plant hormone response sites and abiotic stress response sites were predicted in promoters of some RrLBDs. In them, five RrLBDs (RrLBD12c, RrLBD25, RrLBD39 and RrLBD40) were significantly induced or depressed by salt stress. We predicted the five genes as salt response candidate genes of wild rose and the sites in their promotors maybe pointcut for further study. Abstract LATERAL ORGAN BOUNDARIES DOMAIN (LBD) transcription factors are regulators of lateral organ morphogenesis, boundary establishment, and secondary metabolism in plants. The responsive role of LBD gene family in plant abiotic stress is emerging, whereas its salt stress responsive mechanism in Rosa spp. is still unclear. The wild plant of Rosa rugosa Thunb., which exhibits strong salt tolerance to stress, is an ideal material to explore the salt-responsive LBD genes. In our study, we identified 41 RrLBD genes based on the R. rugosa genome. According to phylogenetic analysis, all RrLBD genes were categorized into Classes I and II with conserved domains and motifs. The cis-acting element prediction revealed that the promoter regions of most RrLBD genes contain defense and stress responsiveness and plant hormone response elements. Gene expression patterns under salt stress indicated that RrLBD12c, RrLBD25, RrLBD39, and RrLBD40 may be potential regulators of salt stress signaling. Our analysis provides useful information on the evolution and development of RrLBD gene family and indicates that the candidate RrLBD genes are involved in salt stress signaling, laying a foundation for the exploration of the mechanism of LBD genes in regulating abiotic stress.


Introduction
Transcription factors (TFs) play key roles in plant functional regulatory maps by regulating target gene transcription [1,2]. The plant-specific TF family LATERAL ORGAN BOUNDARIES DOMAIN (LBD) was identified by the conserved DNA binding domain LATERAL ORGAN BOUNDARIES (LOB), which includes a zinc-finger-like motif composed of cysteine residues "CX2CX6CX3C," a Gly-Ala-Ser (GAS) block, and a spiral coiled structure similar to leucine (Leu) zipper related to protein dimerization, "LX6LX3LX6L" [3]. In particular, the Class II LBD protein only contains the zinc-finger-like structure in the LOB domain [3,4]. The LBD gene families of numerous plant species have been identified genome-wide. A total of 43 LBD genes were first discovered in Arabidopsis thaliana [5]. In addition, 44 LBD genes were identified in the monocotyledonous plant Zea mays [6].
(LBDs of R. rugosa) or RcLBDs (LBDs of R. chinensis), including the potential salt-responsive RrLBD members, have not been identified. This study aimed to screen the RrLBD genes involved in salt stress response by genome-wide phylogenetic analysis, chromosome location analysis, gene structure, and expression profile analysis. This study provides LBD gene candidates for further research on salt stress mechanism regulation in R. rugosa and will be helpful in the cultivation of salt-tolerant rose species by genetic engineering.

Phylogenetic Analysis of LBD Family
All the members of RcLBD family were identified from the genome of R. chinensis by the same method. AtLBDs were downloaded from the PlantTFDB (http://planttfdb. gao-lab.org/, accessed on 15 April 2021) and checked in accordance with a reference study [33]. A neighbor-joining (NJ) tree of protein sequences of RrLBDs, AtLBDs, and RcLBDs was built by MEGA-X (https://www.megasoftware.net/, accessed on 15 May 2021) with p-distance and pairwise deletion parameters [39,40]. A total of 1000 bootstrapping replications were used to verify the reliability of the phylogenies, and the online tool ITOL (https://itol.embl.de/, accessed on 15 July 2021) was used for coloring [41].

Genen Structure, Motif Composition, and Cis-Acting Elements Analysis of RrLBD Family
The motifs were predicted by the online tool MEME 5.3.3 (https://meme-suite.org/ meme/, accessed on 15 May 2021) with default parameters. The cis-acting elements of RrLBD genes were predicted from promoter regions (2000 bp upstream the gene loci) by PlantCare [42]. The gene structures, motifs, and cis-acting elements were visualized by Tbtools 1.086 [43].

Synteny Analysis of LBD Genes
The inter-species synteny analysis was conducted by reciprocal BLASTP search for potential homologous gene pairs (E < 10 −5 , top five matches). The syntenic region and homologous pairs of LBD genes were illustrated by TBtools 1.086 [43].

Expression Analysis of RrLBDs in Different Tissues under Salt Stress
The wild shrubs of R. rugosa, which are naturally distributed in the coastal beach of Western hills north village of Muping district, Yantai city, Shandong province, China (37.455 • N, 121.692 • E), were selected as plant materials. In July 2019, three wild plants were dug out from the sand, and their roots were treated with 170 mM NaCl solution immediately for 1 h (ST). Another three plants that were soaked in deionized water for 1 h were set as the control group (CK). The leaves and roots of ST and CK plants (three replications for each) were picked and frozen by liquid nitrogen immediately and stored at −80 • C. The RNA for each sample was extracted by MiniBEST Universal RNA Extraction Kit (TaKaRa, Japan), and 12 libraries, i.e., three replications for the leaves of ST (L1h), roots of ST (R1h), leaves of CK (LCK), and roots of CK (RCK), were constructed for transcriptome sequencing (RNA-seq) on an Illumina NovaSeq 6000 platform. After filtering raw reads, the fragments per kilobase of exon model per million mapped fragments (FPKM) was calculated by mapping clean reads to the R. rugosa genome using HISAT 2.2.1 (http://daehwankimlab.github.io/hisat2/, accessed on 15 May 2021). The expression profile of RrLBDs was constructed with FPKM. Foldchanges of RrLBDs were calculated on the basis of the read reads using R package DESeq2 (http://www.bioconductor.org/ packages/release/bioc/html/DESeq2.html, accessed on 15 May 2021).
Quantitative real-time polymerase chain reaction (qRT-PCR) was conducted using ChamQ SYBR Color qPCR Master Mix (Vazyme, China) on a CFX96 Real-time PCR platform (Bio-Rad, China). The PCR reaction system and program were conducted following the manufacturer's instructions. Three biological replicates with three technical replicates were prepared for each sample [44]. Table S4 lists the primers of target genes and internal reference gene.

Identification and Phylogenetic Analysis of RrLBD Family
A total of 41 RrLBD TFs were identified by PF03195 model, and their LOB domains were verified in the CDD and Pfam database. The gene names were numbered in accordance with the top one BLASTP search of known AtLBD TFs (Table S1). In addition, homologous RcLBDs were distinguished by lowercase letter suffixes, and isoforms (putative alternative splicing RcLBD mRNA) were distinguished by number suffixes. This situation only existed in RcLBD proteins. The biochemical character of RrLBD proteins changed within a minimal range. A total of 41 RrLBDs proteins ranging from 144 AA (RrLBD24) to 386 (RrLBD27a) AA had a MW ranging from 16.1 kDa (RrLBD12a) to 44.6 kDa (RrLBD27a), and their pI changed from 4.99 (RrLBD7) to 8.91 (RrLBD2b) ( Table S2).
The phylogenetic tree (Figure 1) was constructed by the full-length protein sequences of 41 RrLBDs, 43 RcLBDs, and 43 AtLBDs belonging to two categories (Class I and II), referring to the known topological structure of AtLBD family [45,46]. Seven groups, namely, Ia, Ib, Ic, Id, and Ie of Class I and IIa and IIb of Class II, were divided on the basis of the sub-topological structure. Five homologous gene pairs (LBD37, 38, 39, 40, and 42) of three species, except the missing homologs of AtBD41, indicated that the LBD of Class II was highly conserved. In Class I, Ie, containing four AtLBDs, is a unique class found in A. thaliana. A total of 13, 3, 12, and 8 RrLBDs and 18, 4, 10, and 6 RcLBDs belong to Ia, Ib, Ic, and Id, respectively, which indicated that RrLBDs have a stronger congruent relationship with AtLBDs (10, 3, 12, and 8 for each sub-class).

Gene Structure and Motif Composition Analysis of RrLBD Family
On the basis of the NJ phylogeny of RrLBDs and RcLBDs (Figure 2A), the 10 significantly enriched MEME motifs indicated that the region of the top four motifs with 2, 3-1, and 4 arrangement ( Figure 2B), which overlapped with the complete LOB domains, corresponded to the zinc-finger-like motif (CX2CX6CX3C, Figure 2D (I)), GAS block ( Figure 2D (II)), and the spiral coiled structure similar to Leu zipper (LX6LX3LX6L, Figure  2D (III)), respectively. The four motifs were conserved in all Class I LBDs ( Figure 2B,D)). Class II RrLBDs only contained motif 2 in the location of LOB domains. In addition, motifs 5 and 9 located at the C-terminal were motifs specific to Ia, whereas motifs 6 and 10 specific to Class II LOB domain replaced the location of motifs 3 and 4 of Class I LOB domains, respectively. Compared with AtLBD proteins, in addition to the important proline residue, another glycine residue in the motif 3 specific to RrLBD/RcLBD proteins is highly conserved, and it may be an important residue of the DNA binding domain. This finding showed that different AA residues in the LOB domain are indispensable to the characteristic function of family members, and the protein change in TFs will lead to the changes in its function of recognizing the promoter region [4,47].

Gene Structure and Motif Composition Analysis of RrLBD Family
On the basis of the NJ phylogeny of RrLBDs and RcLBDs (Figure 2A), the 10 significantly enriched MEME motifs indicated that the region of the top four motifs with 2, 3-1, and 4 arrangement ( Figure 2B), which overlapped with the complete LOB domains, corresponded to the zinc-finger-like motif (CX2CX6CX3C, Figure 2D (I)), GAS block ( Figure 2D (II)), and the spiral coiled structure similar to Leu zipper (LX6LX3LX6L, Figure  2D (III)), respectively. The four motifs were conserved in all Class I LBDs ( Figure 2B, 2D)). Class II RrLBDs only contained motif 2 in the location of LOB domains. In addition, motifs 5 and 9 located at the C-terminal were motifs specific to Ia, whereas motifs 6 and 10 specific to Class II LOB domain replaced the location of motifs 3 and 4 of Class I LOB domains, respectively. Compared with AtLBD proteins, in addition to the important proline residue, another glycine residue in the motif 3 specific to RrLBD/RcLBD proteins is highly conserved, and it may be an important residue of the DNA binding domain. This finding showed that different AA residues in the LOB domain are indispensable to the characteristic function of family members, and the protein change in TFs will lead to the changes in its function of recognizing the promoter region [4,47]. The gene structures of RrLBD/RcLBD homolog pairs were highly consistent ( Figure  3C). The variability degree of gene structure is related to the diversity of gene members belonging to the same group. Class I RrLBDs with exons ranging from 1 to 5 (1 to 4 exons for RcLBDs) include 7 intron-free genes, 25 genes with two exons, 3 genes with three exons, 1 gene with four exons, and 1 gene with five exons. Compared with the high gene structure diversity of Class I, Class II is relatively conservative. The coding sequence regions The gene structures of RrLBD/RcLBD homolog pairs were highly consistent ( Figure 3C). The variability degree of gene structure is related to the diversity of gene members belonging to the same group. Class I RrLBDs with exons ranging from 1 to 5 (1 to 4 exons for RcLBDs) include 7 intron-free genes, 25 genes with two exons, 3 genes with three exons, 1 gene with four exons, and 1 gene with five exons. Compared with the high gene structure diversity of Class I, Class II is relatively conservative. The coding sequence regions of most Class II RrLBDs (and RcLBDs) are split by a short intron, with the exception of RrLBD38, which has two more long introns.
The synteny analysis identified 38 LBD orthologous gene pairs of R. rugosa and R. chinensis ( Figure 3B) and 20 pairs of R. rugosa and A. thaliana ( Figure 3C). More pairs of orthologous RrLBD/RcLBD were in accordance with the LBD proteins phylogeny of three
The synteny analysis identified 38 LBD orthologous gene pairs of R. rugosa and R. chinensis ( Figure 3B) and 20 pairs of R. rugosa and A. thaliana ( Figure 3C). More pairs of orthologous RrLBD/RcLBD were in accordance with the LBD proteins phylogeny of three species (Figure 1), and this finding indicated that the lineage division of LBD genes was smaller between rosaceous plants after separation from the ancestor of three species. In particular, RrLBD37 and RrLBD38 of Chr2 and RrLBD39 and RrLBD40 of Chr6 were homologous to at least two RcLBD genes, indicating that Class II members may have functional differentiation between R. rugusa and R. chinensis.

Analysis of Cis-Acting Elements of RrLBD Promotors
A total of 19 cis-acting elements in RrLBD promoter regions were identified. These cisacting elements were classified into five categories: plant hormone response, light response, stress response, specific binding site of MYB, and endosperm expression ( Figure 4). First, the plant hormone response-related category contained the most elements, including methyl jasmonate (MeJA) responsiveness (CGTCA motif), ABA responsiveness (ABRE), salicylic acid responsiveness (TCA element), gibberellin-responsive element (GARE motif), gibberellin-responsive element (TATC box), and auxin-responsive elements (TGA element) or a part of an auxin-responsive element (TGA box). ABRE is the most widely distributed with 82 sites, whereas TGA box is distributed in 3 sites. The light response-related category contained the second highest number of elements, including light-responsive element (GT1 motif), light responsiveness (G box), a part of a module for light response (AE box), and a part of a light response element (CAG motif). Meanwhile, the stress response category including four elements, namely, wound-responsive element (WUN motif), defense and stress responsiveness (TC-rich repeats), low-temperature responsiveness, and enhancer-like elements, are involved in anoxic-specific inducibility (GC motif). The specific binding site of MYB category included the MYB binding sites (MBS) involved in drought inducibility, flavonoid biosynthetic gene regulation (MBSI), and light responsiveness (MRE). The last category included the endosperm expression element (GCN4 motif) with three sites in all promotors. These predicted binding sites or elements indicated that RrLBDs may be targeted by related TFs involved in response to plant hormone, light, abiotic stress, and endosperm development.

Salt-Responsive Expression Analysis of RrLBDs
To screen the salt stress-responsive RrLBD candidates, we screened eight RrLBDs (with one isoforms) belonging to differentially expressed genes out from roots (R1H versus RCK) and leaves (L1H versus LCK) by RNA-seq data. The expression profile based on normalized FPKM values ( Figure 5A, Table S3) manifested that RrLBD40, RrLBD25, RrLBD19, RrLBD12c, and RrLBD4.1 in roots (R1H) and RrLBD40, RrLBD39, and RrLBD38 in leaves (L1H) were induced significantly by salt. The abundance of five RrLBDs in leaves were extremely low for accurate detection (FKPM < 1), and this finding indicated that most RrLBDs especially respond to salt in roots. The relative expression level determined in qRT-PCR ( Figure 5B) proved that the profile based on transcriptome data was credible and indicated the statistical significance of four salt-responsive candidates. RrLBD40 was significantly induced in roots and leaves, RrLBD39 was significantly induced only in leaves, and RrLBD12c was significantly induced only in roots. In particular, as the only member significantly depressed by salt in roots, RrLBD25 can be a negative regulator of salt response by targeting salt-sensitive genes. The four genes can be candidates of strong response to salt stress.

Salt-Responsive Expression Analysis of RrLBDs
To screen the salt stress-responsive RrLBD candidates, we screened eight RrLBDs (with one isoforms) belonging to differentially expressed genes out from roots (R1H versus RCK) and leaves (L1H versus LCK) by RNA-seq data. The expression profile based on normalized FPKM values ( Figure 5A, Table S3) manifested that RrLBD40, RrLBD25, RrLBD19, RrLBD12c, and RrLBD4.1 in roots (R1H) and RrLBD40, RrLBD39, and RrLBD38 in leaves (L1H) were induced significantly by salt. The abundance of five RrLBDs in leaves were extremely low for accurate detection (FKPM < 1), and this finding indicated that most RrLBDs especially respond to salt in roots. The relative expression level determined in qRT-PCR ( Figure 5B) proved that the profile based on transcriptome data was credible and indicated the statistical significance of four salt-responsive candidates. RrLBD40 was significantly induced in roots and leaves, RrLBD39 was significantly induced only in leaves, and RrLBD12c was significantly induced only in roots. In particular, as the only member significantly depressed by salt in roots, RrLBD25 can be a negative regulator of salt response by targeting salt-sensitive genes. The four genes can be candidates of strong response to salt stress.

High Conservation of RrLBD Family
A total of 41 RrLBDs were located on all the seven chromosomes of R. rugosa genome. The gene number is similar to that of R. chinensis, A. thaliana, Z. mays, Solanum lycopersicum, and M. pumila. Thus, the gene number of the LBD family is highly conservative with minimal interference from ancient polyploidization events involving diploid angiosperms [48]. Several paralogous gene pairs located on RrLBD clusters of Chr2 and Chr3 indicated that these paralogous RrLBDs should evolve from tandem replication. In addition to the stable gene number, the lineage of the RrLBD family (two groups including seven subgroups) was consistent with the classification of LBD genes in other species. The

High Conservation of RrLBD Family
A total of 41 RrLBDs were located on all the seven chromosomes of R. rugosa genome. The gene number is similar to that of R. chinensis, A. thaliana, Z. mays, Solanum lycopersicum, and M. pumila. Thus, the gene number of the LBD family is highly conservative with minimal interference from ancient polyploidization events involving diploid angiosperms [48]. Several paralogous gene pairs located on RrLBD clusters of Chr2 and Chr3 indicated that these paralogous RrLBDs should evolve from tandem replication. In addition to the stable gene number, the lineage of the RrLBD family (two groups including seven subgroups) was consistent with the classification of LBD genes in other species. The homologous pairs clustered in the same branch generally have conserved LOB domains whose roles may be similar to those of the corresponding well-studied AtLBDs [49].

Salt Response of RrLBD Candidates
Soil salinity is one of the main environmental stress factors affecting plant growth and development [50]. In the past decade, TFs, including ERF, MYB, WRKY, NAC, and bHLH, have been proven to respond to salt and regulate downstream response gene expression [51,52]. The transcriptional regulation mechanism of the LBD gene family under salt stress is still unclear. Most RrLBD genes express highly in specific tissues. Except RrLBD39 and RrLBD12c, RrLBD25, RrLBD4.1, RrLBD4.2, RrLBD13, RrLBD19, RrLBD38, and RrLBD40 ( Figure 5A) are significantly more abundant in roots. These tissue-specific RrLBDs are up-or downregulated in roots and/or leaves under salt stress, and several genes should be salt-responsive RrLBD candidates.
First, for the significant depression in roots and leaves under salt stress, RrLBD25 may negatively respond to salt stress signal. Further, the defense and stress responsiveness elements in the promoter region of RrLBD25 and its homologous gene AtLBD25, which plays a transcriptional role in auxin signaling, indicate that RrLBD25 may respond to salt by cross talk with auxin signaling [53]. An abiotic stress response study of auxin-related gene families in S. bicolor indicated that SbLBD32 was highly induced by salt and IAA. This result implies that pathways cross talked between auxin and abiotic stress signaling.
The expression of RrLBD39 in leaves and RrLBD40 in roots changed acutely under salt stress. In addition to the highly induced abundance in salt-treated leaves, two MeJA responsiveness elements (CGTCA motif) of the RrLBD39 promoter region indicated that RrLBD39 is a positive regulator of salt. Abiotic stresses, such as drought and salt stresses, can promote the production of secondary metabolites in plants [54,55]. Class IIb gene RrLBD39 is homologous to AtLBD37, AtLBD38, and AtLBD39, regulating anthocyanin synthesis and nitrogen metabolism [26]. In Solanum tuberosum, StLBD1-5 is the homologous gene pair of AtLBD39, and its expression is downregulated under drought stress. This condition is due to the low expression of StLBD1-5, which promotes the accumulation of anthocyanins in leaves, thus improving the plant's drought resistance [56]. The salt response and potential metabolism regulation role of RrLBD39 indicate the balance or promoting relationship of salt tolerance and metabolism of R. rugosa. Moreover, MeJA is a kind of plant hormone that regulates defense mechanism and stress response [57]. MeJA can promote the transcription modification of several genes, thus participating in the resistance to salt stress by promoting the secondary metabolism of plants [58]. In C. sinensis, the expression analysis of CsLBDs after MeJA treatment showed that the expressions of CsLBD38 and CsLBD39_2 were upregulated, and both belong to the Class II subgroup. Combined with the research evidence of this experiment, RrLBD39 may be connected in series with the MeJA signaling pathway to participate in plant secondary metabolism and defense response to salt stress as hormone mediation. Finally, unlike the homologous gene of AtLBD41 regulating the cell specialization of the paraxial region [59], RrLBD40 has a large number of cis-acting elements involved in ABA responsiveness. ABA is an important plant stress hormone, and it plays an important role in salt stress signaling [60]. ABA in plants can interact with MYB TFs and play a role as a negative regulator of salt stress [61]. In addition, several transcriptomics studies showed that most genes of the ABA synthesis and signaling pathway are up-or downregulated under salt stress, and ABA signal pathway is involved in the regulation of salt stress-related genes [62,63]. ABA signaling pathway can also increase the activity and expression of ion transporters, thus regulating the ion homeostasis of plants [64]. The nearly 10 times increased expression of AtLBD40 in salt-treated roots may be amplified by the ABA signaling pathway to actively respond to salt stress signals. In addition, different binding sites of MYB TFs were observed in the promoter regions of RrLBD39 and RrLBD40. Thus, RrLBD39 and RrLBD40 may be target genes of MYB TFs. In Salvia miltiorrhiza, SmLBD50 interacts with SmMYB36/97 protein and participates in jasmonate signal transduction [65]. Therefore, on the basis of previous studies and this experiment, we find that RrLBDs may have a complex molecular regulatory network with other TFs, and this association is important for abiotic stress resistance, growth and development, secondary metabolism, and other important plant biological processes.
According to collinearity analysis, similar to RrLBD25 and AtLBD25, RrLBD4.1 and AtLBD4, RrLBD4.2 and AtLBD3, and RrLBD19 and AtLBD19 are homologous gene pairs. PtLBD4 is involved in the regulation of secondary growth in poplar, whereas AtLBD19 regulates callus formation; however, the specific function of LBD3 has not been clarified [66,67]. Other candidate genes, such as RrLBD12 and RrLBD13, need to be studied to determine their corresponding salt stress mechanisms.

Conclusions
In this study, we identified 41 RrLBD genes from the R. rugosa genome by bioinformatics. According to phylogenetic analysis, all RrLBD genes were divided into seven subclasses of two classes. The conserved motifs, gene structures, cis-acting elements, and collinearity analysis showed the conservation of the RrLBD gene family. The expression profiles of roots and leaves under salt stress indicated that RrLBD25, RrLBD4.1, RrLBD4.2,