The Soybean Laccase Gene Family: Evolution and Possible Roles in Plant Defense and Stem Strength Selection

Laccase is a widely used industrial oxidase for food processing, dye synthesis, paper making, and pollution remediation. At present, laccases used by industries come mainly from fungi. Plants contain numerous genes encoding laccase enzymes that show properties which are distinct from that of the fungal laccases. These plant-specific laccases may have better potential for industrial purposes. The aim of this work was to conduct a genome-wide search for the soybean laccase genes and analyze their characteristics and specific functions. A total of 93 putative laccase genes (GmLac) were identified from the soybean genome. All 93 GmLac enzymes contain three typical Cu-oxidase domains, and they were classified into five groups based on phylogenetic analysis. Although adjacent members on the tree showed highly similar exon/intron organization and motif composition, there were differences among the members within a class for both conserved and differentiated functions. Based on the expression patterns, some members of laccase were expressed in specific tissues/organs, while some exhibited a constitutive expression pattern. Analysis of the transcriptome revealed that some laccase genes might be involved in providing resistance to oomycetes. Analysis of the selective pressures acting on the laccase gene family in the process of soybean domestication revealed that 10 genes could have been under artificial selection during the domestication process. Four of these genes may have contributed to the transition of the soft and thin stem of wild soybean species into strong, thick, and erect stems of the cultivated soybean species. Our study provides a foundation for future functional studies of the soybean laccase gene family.


Introduction
Laccases, p-diphenol:oxygen oxidoreductases (EC 1.10.3.2), represent the largest subgroup of multicopper-containing oxidases (MCOs), also known as "blue" oxidases [1][2][3]. They are widespread not only in fungi and plants, but also in bacteria and insects [4][5][6][7]. Laccase was initially discovered in Japanese lacquer tree by Yoshida in 1883 [8]. Subsequently, it was defined as a metal-containing oxidase [9]. Commonly, the laccase protein sequence has three conserved domains (Cu-oxidase, Cu-oxidase_2, and Cu-oxidase_3) that are used to identify canonical laccases. The active center of arabidopsis.org). In total, 17 laccase sequences were obtained according to published research [3]. Arabidopsis laccase members contain Cu-oxidase, Cu-oxidase_2, and Cu-oxidase_3 (PF00394, PF07731, and PF07732) domains. The three domains were searched in the soybean genome using HAMMER 3.0. We took the intersection of the three BLAST results and manually eliminated genes with no open reading frame and multi-transcripts from the same gene. Additionally, identified genes were verified using the SMART domain database (http://smart.embl-heidelberg.de).

Sequence Alignment and Phylogenetic Analysis
An alignment of 93 GmLac and 17 AtLac peptides was performed using ClustalW in MEGA 7.0 [49]. All conserved sites were used to construct evolutionary relationships with the neighbor-joining (NJ) method (1000 replicates). The final phylogenetic tree was visualized and edited in iTOL (http://itol.embl.de/) [50].

Gene Duplication and Expansion Analysis
Intra-species synteny blocks were detected by the MCScanX algorithm [51], with following settings (five genes required to call a collinear block based on the previous all-to-all BLASTP result E-value ≤ 1 × 10 −5 ). The ratio of nonsynonymous substitutions per nonsynonymous site (Ka) to synonymous substitutions per synonymous site (Ks) was computed using KaKs_Calculator 2.0 [52]. Then, the estimated date (mya, million years ago) of each duplication event was calculated using the mean Ks values (T = Ks/2λ), assuming clocklike rates (λ) of 6.1 × 10 −9 [53]. Tandem repeat genes were calculated through the physical location on the chromosome, and segmental duplication genes were defined as those located in the same synteny blocks.

Gene Structure and Motif Compositions Analysis
A phylogenetic tree based only 93 GmLac peptides was constructed according to same method in Section 2.3. The exon/intron organizations of soybean laccase genes were identified with the Gene Structure Display Server, GSDS 3.0 (http://gsds.cbi.pku.edu.cn/) [54]. Conserved motifs of laccase proteins were identified statistically using MEME (http://meme-suite.org/tools/meme). Discovery mode was set as discriminative mode and we used 17 Arabidopsis laccase protein sequences as a control. The maximum number of motifs to find was set at 10 [55]. Visualization of motif compositions was executed with TBtools V0.52 (https://github.com/CJ-Chen/TBtools).

Expression Pattern of Soybean Laccase Family Genes
Nine soybean-tissue/organ RNA sequencing (RNA-Seq) data, FPKM (Fragments per kilobase of transcript per million fragments mapped) values, were obtained from the Phytozome database [56]. Average expression valued were calculated, and a three-fold increase compared to the average value was defined as a high-expression gene. Color gradations (light to dark green representing different expression abundances) were appended to the expression data.
2.8. qRT-PCR Analysis of GmLacs Expression Following Phytophthora sojae Infection P. sojae strain (P7076) was posted on the roots of 72-h-old seedings (grown in darkness), collecting infected samples at 0, 0.5, 3, 6, and 12 h. Total RNA of these samples was extracted using the Easy-Pure Plant RNA kit (Transgen, ER301). The RNA quality was detected through an ultra-microspectrophotometer (IMPLEN, P330). RNA samples were reverse-transcribed with EasyScript One-Step gDNA Removal and complementary DNA (cDNA) Synthesis SuperMix (Transgen, AE311) strictly following the manufacturer's instructions. The gene expression pattern was analyzed through qRT-PCR using SYBR ® Green Master (Roche, 4913914001) in the Agilent Technologies Stratagene Mx3005P Detection Device. Three replicates were investigated to confirm the result. The relative expression levels were calculated via the 2 −∆∆Ct method after taking GmPDF (serine/threonine-protein phosphatase 2A regulatory subunit A gene, Glyma.20G114000) as an internal control. The qRT-PCR analysis primers are listed in the Table S8.

Selective Pressure Analysis
SNP (Single nucleotide polymorphism) data gathered from resequencing of 302 soybean genomes were used for the genetic diversity analysis of GmLac genes in soybean [57]. The SNPs with missing data >10% or MAF (Minor Allele Frequency) <5% were filtered. The soybean accessions were divided into two populations: Glycine soja and G. max (landrace and improved cultivar). Nucleotide diversity (π) was calculated using a 20-k-10-k sliding window and VCFtools [58]. The ratio of diversity π G.max /π G.soja for each window was calculated. The π G.max /π G.soja ratio of a target region was compared to the mean π G.max /π G.soja ratio of the corresponding entire chromosome to determine whether the target region was under selection pressure or not [59]. The 5% windows with the lowest π G.max /π G.soja ratios values were compared with the target region π G.max /π G.soja ratios to determine the genes that were highly selected [60].

Identification and Classification of Laccase Genes in Soybean
To identify the laccase genes in soybean, we downloaded the reference genome of Glycine max from Phytozome and searched the genome with HAMMER 3.0 software for Cu-oxidase, Cu-oxidase_2, and Cu-oxidase_3 domains (PF00394, PF07731, and PF07732). After elimination of the redundant genes with only one or two typical domains on the peptide sequence or with no integral open reading frames, 93 putative soybean laccase-coding genes GmLac1 through GmLac93 were identified (Table S1). Furthermore, we obtained the physical locations of these laccase genes and mapped them onto 19 of the 20 soybean chromosomes (Figure 1). Chromosome 18 had the greatest number of laccase genes, which was 12 (GmLac70-GmLac81). On the contrary, chromosomes 3, 9, 15, and 16 contained only one gene each. GmLac91 and GmLac92 mapped onto scaffold_27, while GmLac93 mapped onto scaffold_614. Basic information of all soybean laccases including gene name, gene locus, amino-acid length, signal peptide length, molecular weight, pI value, and subcellular localization was also determined and is presented in the Table S1. The amino-acid length ranged from 315 to 637 aa, with molecular weight ranging from 34,958.38 to 71,749.27 Daltons, and isoelectric points (pI) ranging from 4.95 to 9.96. Most soybean laccase pI values were above 7.00, whereas only 19 laccase pI values were below 7.00. The predicted subcellular locations of 76 soybean laccases included the extracellular space. Thirty-four of the laccases were predicted to localize only to the extracellular space. Thus, most laccases are probably secretory proteins. Laccases may also commonly be located in the plasma membrane and less commonly in lysosomes. Finally, very few laccases localized to other organelles including the mitochondria, peroxisomes, and vacuole. Specifically, GmLac21, GmLac39, and GmLac86 may be located in the mitochondria and peroxisomes, whereas GmLac77 may be located only in the mitochondria.
GmLac90 was the only one predicted to be located in the vacuole. Transmembrane topology and signal peptide analysis indicated that nine laccases contain transmembrane domains. GmLac29 contains two transmembrane regions. The hydrophobic region (H-Region) of signal peptides of most laccase members ranged from 7-20 aa (Table S1).
Genes 2019, 10, x FOR PEER REVIEW 5 of 20 extracellular space. Thus, most laccases are probably secretory proteins. Laccases may also commonly be located in the plasma membrane and less commonly in lysosomes. Finally, very few laccases localized to other organelles including the mitochondria, peroxisomes, and vacuole. Specifically, GmLac21, GmLac39, and GmLac86 may be located in the mitochondria and peroxisomes, whereas GmLac77 may be located only in the mitochondria. GmLac90 was the only one predicted to be located in the vacuole. Transmembrane topology and signal peptide analysis indicated that nine laccases contain transmembrane domains. GmLac29 contains two transmembrane regions. The hydrophobic region (H-Region) of signal peptides of most laccase members ranged from 7-20 aa (Table S1).

Phylogenetic Analysis of the Soybean Laccase Gene Family
To further investigate the phylogenetic relationships of laccases between soybean and Arabidopsis, we constructed a neighbor-joining phylogenetic tree based on multiple sequence alignments of 93 putative soybean laccase and 17 Arabidopsis laccase peptide sequences ( Figure 2). According to classification standard of Arabidopsis laccases and location on the phylogenetic tree [61], 93 soybean laccases were divided into five groups. The distribution of the laccase genes in each group was rather uneven; 15 members were classified as class I, 13 members as class Ⅱ, eight members as class Ⅲ , eight members as class Ⅳ, and 49 members as class Ⅴ. Unlike in Arabidopsis, soybean laccase members in class Ⅳ classified into two subgroups: IVa subgroup (Gm18 and Gm19) and IVb subgroup (the six other genes). None of the soybean laccases clustered into class Ⅵ. The specific gene list of each category is presented in the Table S2a-c.

Phylogenetic Analysis of the Soybean Laccase Gene Family
To further investigate the phylogenetic relationships of laccases between soybean and Arabidopsis, we constructed a neighbor-joining phylogenetic tree based on multiple sequence alignments of 93 putative soybean laccase and 17 Arabidopsis laccase peptide sequences ( Figure 2). According to classification standard of Arabidopsis laccases and location on the phylogenetic tree [61], 93 soybean laccases were divided into five groups. The distribution of the laccase genes in each group was rather uneven; 15 members were classified as class I, 13 members as class II, eight members as class III, eight members as class IV, and 49 members as class V. Unlike in Arabidopsis, soybean laccase members in class IV classified into two subgroups: IVa subgroup (Gm18 and Gm19) and IVb subgroup (the six other genes). None of the soybean laccases clustered into class VI. The specific gene list of each category is presented in the Table S2a-c.

Duplication and Expansion Events in Soybean Laccase Gene Family
Gene duplication is considered one of the major driving forces of genome evolution [62]. Duplicated genes are the source for creating novel genetic variation. It is generally accepted that the soybean genome went through three rounds of whole-genome duplication (WGD) events and retained more than 50% of the repeated segments [63,64]. They were the γ WGD event about~130 to 240 million years ago (Mya), the legume WGD event about~58 Mya, and the Glycine genus WGD event about~13 Mya [65][66][67]. The general Ks rate range of soybean genome γ WGD events is bigger than 1.5, while that of the legume WGD event is 0.3 to 1.5, and that of the Glycine genus WGD event is 0 to 0.3 [40,44,47]. For an in-depth study the WGD events of laccase genes, we conducted collinear analysis using the base MCScanX software. In total, 85 orthologous gene pairs were identified within the pairwise syntenic blocks (Figure 3). We calculated the non-synonymous nucleotide substitution rate (Ks) and synonymous nucleotide substitution rate (Ka) between collinear gene pairs (Table S3). We calculated that about 25.9% of gene pairs (22 of 85) were separated during the Glycine genus WGD event, nearly 43.5% of gene pairs (37 of 85) were separated during the legume WGD event, and approximately 30.6% of gene pairs (26 of 85) were separated during the γ WGD event. Additionally, we also calculated the Ka/Ks ratios to study the gene divergence after the WGD events. We observed that all 85 synteny pairs underwent purifying selection with Ka/Ks < 1 [48].
The most common patterns of gene duplication include interspersed repeats, tandem repeats, and segmental duplication [68], while segmental and tandem duplications are usually considered as the main causes of large gene family expansion in plants [69]. To gain insight into the soybean laccase gene duplication pattern, we analyzed the tandem and segmental duplication events of laccase genes. To identify tandem duplication, we searched adjacent gene clusters on the same chromosome with no more than one interferential gene being inserted. We found that 43.0% (40 of 93) of genes in the soybean laccase gene family were tandem repeats ( Figure 1 and Table S4).

Duplication and Expansion Events in Soybean Laccase Gene Family
Gene duplication is considered one of the major driving forces of genome evolution [62]. Duplicated genes are the source for creating novel genetic variation. It is generally accepted that the soybean genome went through three rounds of whole-genome duplication (WGD) events and  The most common patterns of gene duplication include interspersed repeats, tandem repeats, and segmental duplication [68], while segmental and tandem duplications are usually considered as the main causes of large gene family expansion in plants [69]. To gain insight into the soybean laccase gene duplication pattern, we analyzed the tandem and segmental duplication events of laccase genes. To identify tandem duplication, we searched adjacent gene clusters on the same chromosome with no more than one interferential gene being inserted. We found that 43.0% (40 of 93) of genes in the soybean laccase gene family were tandem repeats ( Figure 1 and Table S4).
Segmental duplications are blocks of genomic sequence ranging from 1-200 kb that map to different loci in a genome and share high sequence identity (often >90%) [70][71][72]. In this study, we defined segmental duplication as genes located on the same homologous blocks [73]. Thus, collinearity gene pairs can be regarded as segmental duplication members in the laccase gen family. In total 45.1% (42 of 93) of laccase genes involved segmental duplication (Table S4).

Gene Structure and Motif Compositions of GmLacs
To reveal the structural diversity of soybean laccase genes, we constructed the exon/intron organizations and searched for conservative motifs based on the phylogenetic tree of all soybean laccase alignments ( Figure 4). Overall, nearly all the closest genes on the phylogenetic tree showed remarkably similar gene structures. However, there was still a small proportion of the gene clades that exhibited different intron/exon organizations. For instance, GmLac13 contained six exons, Segmental duplications are blocks of genomic sequence ranging from 1-200 kb that map to different loci in a genome and share high sequence identity (often >90%) [70][71][72]. In this study, we defined segmental duplication as genes located on the same homologous blocks [73]. Thus, collinearity gene pairs can be regarded as segmental duplication members in the laccase gen family. In total 45.1% (42 of 93) of laccase genes involved segmental duplication (Table S4).

Gene Structure and Motif Compositions of GmLacs
To reveal the structural diversity of soybean laccase genes, we constructed the exon/intron organizations and searched for conservative motifs based on the phylogenetic tree of all soybean laccase alignments (Figure 4). Overall, nearly all the closest genes on the phylogenetic tree showed remarkably similar gene structures. However, there was still a small proportion of the gene clades that exhibited different intron/exon organizations. For instance, GmLac13 contained six exons, whereas its nearby paralogous genes, GmLac68, GmLac67, GmLac40, and GmLac5, had eight exons even though their evolutionary relationships reached 98-100% bootstrap values, respectively.
To further reveal the conserved motifs of the GmLac proteins, we analyzed 93 GmLac proteins using the MEME program and identified 10 motifs (Figure 4; Figure S1). As expected, the motif compositions of peer groups had similar structures and organizations. Some classes of GmLac proteins with completely the same motif composition may indicate the possibility of functional redundancies among those members. Some motifs were found in all 93 putative members. These conserved motifs may play critical roles in mediating GmLac function. Despite an overall similarity among the 93 GmLac proteins, there were some differences in the number of motifs. For example, the number of motifs varied from 2-10 among those GmLac proteins. Varying numbers of motifs across the GmLac proteins may suggest functional divergence among some of the GmLac proteins. even though their evolutionary relationships reached 98-100% bootstrap values, respectively.
To further reveal the conserved motifs of the GmLac proteins, we analyzed 93 GmLac proteins using the MEME program and identified 10 motifs (Figure 4; Figure S1). As expected, the motif compositions of peer groups had similar structures and organizations. Some classes of GmLac proteins with completely the same motif composition may indicate the possibility of functional redundancies among those members. Some motifs were found in all 93 putative members. These conserved motifs may play critical roles in mediating GmLac function. Despite an overall similarity among the 93 GmLac proteins, there were some differences in the number of motifs. For example, the number of motifs varied from 2-10 among those GmLac proteins. Varying numbers of motifs across the GmLac proteins may suggest functional divergence among some of the GmLac proteins.

Expression Patterns of GmLac Genes
Gene expression patterns often reveal how, when, and where genes are active. To verify the expression profiles of GmLac genes, we analyzed publicly available RNA-Seq data in Phytozome. Various organs and tissues including root, root hair, nodules, pod, seed, stem, shoot apical meristem (SAM), leaves, and flower transcriptomes were downloaded (Table S5a,b). According to the fragments per kilobase million (FPKM) values of 90 GmLacs excluding the three members located in scaffolds, a differential transcriptional regulation pattern for the laccase genes was observed. GmLac25 has the highest expression abundance and is specifically expressed in flowers (Table S5a). The overall mean FPKM value of all 90 GmLac genes among the nine transcriptomes was about 10. We defined the genes with mean values three-fold higher than the overall mean FPKM value across all transcriptomes as highly expressed GmLac genes. With this criterion, we identified 11 genes that were highly expressed in multiple tissues ( Table 1). Most of the highly expressed genes showed transcript abundance in the stem, indicating that they may be involved in stem cell-wall maintenance to provide mechanical strength for the aboveground portion of the plant. Conversely, the expression of some genes was restricted to a specific organ or tissue and very poorly expressed in others. For example, GmLac1 and GmLac24 specifically and highly expressed in the roots, GmLac77 did so in pods, GmLac2, GmLac8, and GmLac9 did so in the stem, GmLac87 did so in the SAM, and GmLac25 did so in the flower ( Table 2). The promoters of these genes were active in specific organs or tissues. In addition, 11 genes were inactive in all tissues and organs studied (Table S5b). For the remaining members, there was either a mild expression pattern with moderate abundance or expression in several samples.  We also compared the expression pattern of collinear gene pairs (Table S6). In total, 64.7% (55 of 85) of collinear gene pair expression patterns were distinct. Only a fraction of gene pairs (23 of 85) maintained a similar expression pattern. These results additionally suggested a differentiation of the promoter elements of laccase genes and possibly their functions.

Expression Changes of GmLacs after Inoculation with Phytophthora sojae
It is well established that laccases are involved in lignin metabolism and plant defenses [34][35][36]74,75]. Phytophthora sojae is an oomycete pathogen that causes root and stem rot disease in soybean, which severely affects yield [76]. We were interested if any laccase genes were rapidly induced following P. sojae infection. An earlier transcriptomic dataset of soybean-P. sojae interaction is available on the NCBI database (https://www.ncbi.nlm.nih.gov/bioproject/PRJNA318321). We studied the responses of 93 soybean laccase genes following infection (Table S7a). Considering 0 h post infection (hpi) as a control, we calculated the expression change value of the GmLac genes at 0.5, 3, 6, and 12 hpi (Table S7b), and the results are depicted in a heatmap ( Figure 5). According to the clustering, 93 putative laccase genes can be roughly divided into three categories based their expression status. Most laccase genes (72 of 93) did not show obvious expression changes following P. sojae infection. Nine of the remaining were induced, including Gmlac28, Gmlac91, Gmlac92, Gmlac79, Gmlac23, Gmlac57, Gmlac47, Gmlac37, and Gmlac36. Although they were induced, their expression levels peaked at different stages. Expression abundances of the remaining 12 genes were suppressed. GmLac32 and GmLac75 were strongly repressed at 12 hpi, while the other nine (GmLac78, GmLac71, GmLac25, GmLac89, GmLac15, GmLac90, GmLac45, GmLac48, GmLac76, and GmLac18) were potentially suppressed at multiple stages. By examining the expression profile of these 21 genes (Table S5a), we noticed that these genes were barely expressed in normal conditions. On the other hand, after infection, the expression of these genes changed significantly, which strongly implies a close tie between these genes and oomycete pathogen resistance.
To further test the expression pattern of these genes, we performed qRT-PCR on a Williams 82 root sample with a soybean Phytophthora sojae strain (P7076) [77]. The results ( Figure 5) suggested that most (8/9) induced laccase genes were similarly induced with the exception of GmLac36. In particular, GmLac28, GmLac91, GmLac92, GmLac79, GmLac23, and GmLac57 expression abundances were significantly increased in both tests. Conversely, the quantitative results of repressed genes also displayed coincident RNA-seq data. GmLac78, GmLac75, GmLac32, GmLac32, GmLac45, GmLac89, and GmLac25 showed obvious repression in expression. GmLac48 showed no expression signal. Automatically, we suggest a putative defense function for the P. sojae-induced genes, while suppressed genes could equally be important, as P. sojae effectors may manage to suppress the expression of some defense-related GmLac genes to promote colonization.

Selection Pressure Analysis of Laccase Gene Family during Soybean Domestication
Cultivated soybean was domesticated from its relative Glycine soja in ancient China [78], which led to a series of "domestication syndromes" and subsequently resulted in markedly decreased nucleotide diversity (π) of the respective domestication genes [79]. Wild soybeans uniformly exhibit a softer and thinner stem, while cultivated soybeans uniformly show thicker, stronger, and erect stems, suggesting that these phenotypes are most likely domestication-related traits [80].
Lignin is abundant in the stem cortex and xylem cell wall, playing an important role in maintaining mechanical strength of stem [81]. Laccase is well known for its ability to facilitate lignin synthesis and promote plant tissue-specific lignification [82]. Therefore, we wanted to investigate whether the domestication of soybean was accompanied by changes in laccase gene nucleotide diversity. Thus, we performed selection pressure analysis of GmLac genes based on π ratio (π G.max /π G.soja ) between 62 wild species and 240 (130 landraces and 110 improved cultivars) cultivated soybean accessions [57]. A comparison of the π G.max /π G.soja values of a chromosomal region containing GmLac genes with the corresponding mean value of the respective chromosome (Table S9a) revealed that 57 of 90 (63.3%) GmLac genes exhibited above average π G.max /π G.soja values. This observation indicates that the nucleic acid diversity of GmLac genes significantly decreased following the domestication of soybean, and also suggests that most laccase genes were under selective pressure during the domestication of cultivated soybean (Table S9b). When judging whether a gene underwent selection pressure by comparing with the top 5% of π G.max /π G.soja values according to previous research [57,60], we identified 10 genes located in the selective sweep regions ( Figure 6; Table S9c). Of these 10 genes, GmLac8, GmLac12, GmLac66, and GmLac85 are highly expressed in the stem, and, more importantly, the expression of GmLac8 is highly specific to the stem (Table S9d). These genes are most likely involved in the transition from a soft and thin soybean stem to a thick and erect stem. To further support our conclusion that 10 genes have underwent strong selection pressure during the domestication of soybean, we conducted phylogenetic analysis using SNPs of the ten genes among 302 resequenced soybean accessions (Supplementary Materials). The 302 accessions can be divided into five groups (Figure 7). In Groups I through IV, all accessions except SRR1533282 are wild soybean species. In Group V, all but SRR1533169 and SRR1533198 are soybean cultivars. These results also imply that the nucleotide diversity of the 10 genes is much higher in the wild accessions than in the cultivated accessions.

Discussion
Laccase enzymes are multifunctional proteins. They were shown to contribute toward cell morphology [32,33], pigment formation [83], secondary cell-wall biosynthesis [84,85], and resistance to biotic and abiotic stresses in plant [74,75]. Functional studies of laccases were already carried out in many model and crop plants, including Arabidopsis [86], Brachypodium distachyon [87], sorghum [42], Oryza sativa [43], Populus trichocarpa [88], Gossypium raimondii, and Gossypium arboretum [74]. However, as one of the most important crops around the world, the functions of soybean laccase are rarely reported. In this study, gene structures, protein properties, phylogenetic relationships, replication characteristics, expression profiles, and selection pressures were investigated for members of the soybean laccase gene family.
The soybean Williams 82 genome encompasses 105 putative copper oxidase domain-containing genes that can be regarded as putative candidate laccase genes. Typical laccases usually have three copper oxidase domains; therefore, we investigated 93 GmLac genes that contained three copper

Discussion
Laccase enzymes are multifunctional proteins. They were shown to contribute toward cell morphology [32,33], pigment formation [83], secondary cell-wall biosynthesis [84,85], and resistance to biotic and abiotic stresses in plant [74,75]. Functional studies of laccases were already carried out in many model and crop plants, including Arabidopsis [86], Brachypodium distachyon [87], sorghum [42], Oryza sativa [43], Populus trichocarpa [88], Gossypium raimondii, and Gossypium arboretum [74]. However, as one of the most important crops around the world, the functions of soybean laccase are rarely reported. In this study, gene structures, protein properties, phylogenetic relationships, replication characteristics, expression profiles, and selection pressures were investigated for members of the soybean laccase gene family.
Through multiple sequence alignment of 17 Arabidopsis and 93 soybean laccases and construction of an N-J phylogenetic tree, we divided the soybean laccases into five classes according to the classification criteria applied in Arabidopsis thaliana. Class IV was divided into two different branches, IVa and IVb. Further genome-wide duplication (GWD), homologous, and historical duplication event analyses of the GmLac genes revealed that laccase genes are an ancient gene family that expanded as early as 321 Mya ago during the γ WGD event that was accompanied by the evolution of rosids. Laccase genes continued to expand in subsequent legume WGD and Glycine WGD events. According to the Ks value, purifying selection was the major force driving the evolution of the laccase gene family. Synteny analysis and the physical location of GmLac genes show that tandem repeats (41.93%) and segmental duplications (43.01%) were the main mechanisms of laccase gene duplication.
Gene structure and protein motif analyses showed that the most closely related members in the phylogenetic tree had similar exon/intron organization and common motif compositions. This may suggest the existence of a possible functional similarities among the closely related GmLac proteins. However, within the five basic classes, exon/intron organization and motif compositions were highly diverse, suggesting the occurrence of functional differentiation in each subgroup of the soybean laccase gene family.
Investigation of the expression profiles of GmLac genes among nine organs/tissues revealed the possible role of transcriptional regulation in mediating the functions of GmLac proteins. Some members have high expression abundance, while some genes are barely expressed. Both organ/tissue-specific and constitutive types of expressions were observed among members of the gene family.
In response to P. sojae infection, only a few GmLac genes were either induced or repressed. Although it is possible to assign a putative defense function for the P. sojae-induced genes, downregulation of some genes could equally be important in Phytophthora resistance. Exchange of the promoter of such a gene was shown to enhance resistance in transgenic soybean plants against this pathogen [96]. P. sojae effectors may manage to suppress the expression of some defense-related GmLac genes to avoid lignified cell barriers. Laccase genes are able to enhance plant resistance through accelerating the polymerization of monolignol to lignin. On one hand, increased lignin content helps create stronger physical barriers to prevent invasion; on the other hand, lignin synthesis reinforces cell biosynthesis, which is a branch of the phenylpropyl pathway, and phenylpropanoid metabolism is one of the most important secondary metabolic pathways involved in plant defense against biotic and abiotic stresses [75].
Our study on the natural variations in the GmLac gene family across a collection of both wild and cultivated accessions revealed that the laccase gene family may have been subjected to selection pressure following domestication ( Figure 6). Together with the differential expression patterns and selection pressure, four genes were found to be highly expressed in the stem. Based on the differences in stem characteristics between wild soybean and cultivated soybean and laccase's role in lignification, we speculate that these four genes may have made an important contribution to the transition from a soft and thin stem in wild species to a hard, thick, and erect stem in present day soybean cultivars.

Conclusions
A total of 93 putative GmLac genes were identified from the William 82 genome. Based on the phylogenetic relationships between these genes and Arabidopsis thaliana laccase genes, they were divided into five classes (Figures 1 and 2; Tables S1 and S2). Gene duplication and expansion analysis showed that, as expected, duplication of the members of the soybean laccase gene family occurred in all three WGD events in soybean ( Figure 3; Table S3). Purifying selection was the major driving force in laccase gene family evolution. The main gene duplication mechanisms for GmLac gene family were tandem and segmental duplications. The analyses of conserved domains suggested that, in general, adjacent members in the phylogenetic tree had common motif compositions. However, exceptions were observed within a clade for motif compositions. Differences were recorded for functional conservation, as well as in the motif compositions within clades (Figure 4; Table S4). Transcriptomic analysis revealed both tissue-and organ-specific specific expression patterns, as well as constitutive expression patterns, of the soybean laccase gene family (Tables 1 and 2; Table S5). This study also revealed that distinct GmLac gene sets are induced or repressed in response to invasions by oomycete pathogen P. sojae ( Figure 5; Table S7). Selective pressure analysis revealed a possible role of a few members of the gene family in transitioning from the soft and thin stem of wild soybean species to the strong, thick, and erect stem of the cultivated soybean accessions.
Supplementary Materials: The following are available online at http://www.mdpi.com/2073-4425/10/9/701/s1, Table S1: Basic information of soybean laccases; Table S2a: Gene numbers in each phylogenetic class in soybean and Arabidopsis thaliana, Table S2b: Specific gene list of each phylogenetic class of Arabidopsis thaliana, Table S2c: Specific gene list of each phylogenetic class of soybean; Table S3: Collinear gene pairs and their Ka, Ks value and estimated divergence time; Table S4: List of tandem repeat and segmental repeat laccase genes; Table S5a: Expression pattern of soybean laccase gene family; Table S5b: Genes that are barely expressed in 9 detecteted materials; Table S6: Comparison of expression patterns between collinearity pairs; Table S7a: Soybean laccase gene expression pattern following P. sojae infection; Table S7b: Change value of laccase gene inoculation with P. sojae ( FPKM); Table S8: Primers used for qRT-PCR; Table S9a: Sliding window analysis of nucleotide diversity of 93 laccase gene sequence; Table S9b: List of selection GmLacs above average value; Table S9c: List of GmLacs experience sellection pressure above top 5% genome-wide level; Table S9d: Expression pattern of the 10 highly selected GmLacs; Figure S1: Schematic diagram of the motifs of soybean laccase proteins; Data S1: SNPs (filtered) of all laccase gene of 302 soybean accessions.