The genome of Populus trichocarpa v3.0 and Arabidopsis were downloaded from Phytozome (https://phytozome.jgi.doe.gov) (Goodstein et al., 2012). According to published research, Arabidopsis thaliana has a total of 17 laccase genes (Mccaig et al., 2005), this number is 5 in Zea mays (Caparrós-Ruiz et al., 2006), and 53 in Populus trichocarpa (Bryan et al., 2016). Each of the 75 sequences has three domains, which named Cu-oxidase, Cu-oxidase_2, and Cu-oxidase_3 (PF00394, PF07731, and PF07732). An alignment of 17 AtLAC peptides, 5 ZmLAC peptides,53 PtrLAC peptides, 46 P. deltoides laccase peptides and 12 P. tomentosa laccase peptides were performed using ClustalW in MEGA 10.0 with the default parameters. Full-length sequences of all proteins were used to construct phylogenetic tree with the neighbor-joining(NJ) method (1000 replicates) (Kumar et al., 2018). Finally, the tree was annotated and decorated by ITOL online tools (http://itol.embl.de/) (Letunic and Bork, 2016).

The evolutionary tree construction method of 53 genes was the same as shown in 2.1. The exon-intron organization of laccase genes was identified by the online program Gene Structure Display Server 2.0 (http://gsds.cbi.pku.edu.cn) (Hu et al., 2015). The MEME online program (http://meme.nbcr.net/meme/intro.html) was used to identify conserved motifs in Populus trichocarpa laccases (Bailey et al., 2009), and the maximum number of motif was set to 15. The domain information of each Populus laccase protein was obtained from Pfam database (https://pfam.xfam.org/search). Visualization of the motif and domain position was accomplished by TBtools software (https://github.com/CJ-Chen/TBtools) (Chen et al., 2020).

Chromosome length and the distribution of 53 LAC genes of P. trichocarpa were obtained by using the Phytozome database (https://phytozome.jgi.doe.gov). Multiple Collinearity Scan toolkit (MCScanX) (Wang et al., 2012) was adopted to analyze the gene duplication events with the default parameters. The Gene Location Visualize tool of TBtools software (https://github.com/CJ-Chen/TBtools) was used to Visualize the results of tandem repeat analysis, the Circle Gene View function in TBtools was used to exhibit the collinearity of laccase gene family in P. trichocarpa, the Dual Systeny Plot in TBtools was used to show the synteny relationships between P. trichocarpa and other 5 species. Ka/Ks Calculator 2.0 was used to calculate the Ka/Ks ratio of PtrLACs to reveal the natural selection pressure experienced by laccase genes during the evolution (Wang et al., 2010).

The RNA-Seq data and FPKM (Fragments per kilobase of transcript per million fragments mapped) values of 40 samples including different tissues and developmental stages were obtained from phytozome database. The expression level of PtrLAC genes were presented in the form of heatmap by TBtools software.

The upstream sequences (2000 bp) of 53 PtrLAC genes were extracted by PGSC. Then, PlantCare (http://bioinformatics.psb.ugent.be/webtools/plantcare/html/) was used for predictive analysis (Lescot et al., 2002), and 9 cis-acting elements were obtained after screening. They are abscisic acid responsiveness, auxin responsiveness, defense and stress responsiveness, drought-inducibility, gibberellin-responsiveness, light responsiveness, low-temperature responsiveness, MeJA-responsiveness, MYB binding site involved in light responsiveness.

The number of amino acids per protein was obtained using the phytozome database, the molecular weight (Da) and isoelectric point (pl) of the protein were predicted by using the online ExPASy tool (https://web.expasy.org/protparam/) (Gasteiger et al., 2003). Subcellular locations of PtrLAC members were calculated by the online software CELLO (http://cello.life.nctu.edu.tw/).The signal peptides of PtrLACs were predicted by using the SignalP online tool (http://www.cbs.dtu.dk/services/SignalP/) (Almagro Armenteros, 2019).

Protein structure was predicted by AlphaFold2. The structure prediction process was roughly as described in AlphaFold Paper 2, consisting of five steps:MSA construction, template search, five model reasoning, average PLDD-based model ordering, and constraint relaxation of the predicted structure. Structure factors of ZmLAC3, which coded as 6KLG (ZmLAC3 native), 6KLI (ZmLAC3–SinA complex) and 6KLJ (ZmLAC3–ConA complex), were downloaded from Protein Data Bank (http://www.rcsb.org/) (Andreaeopsida, 2013). AutoDock Vina software were used for molecular docking of PtrLAC with SinA and ConA (Trott and Olson, 2009). Structural comparison between ZmLAC3 and PtrLACs as well as Polar interaction between monolignols and amino acid were analyzed by PyMol 2.0 (Seeliger and De Groot, 2010).

45-day-old 84K poplars (P.alba×P.glandulosa) were used as materials for qRT-PCR. Primers gaattcgagctcGGTACCattcaaacctgcgttgatcc and aagcttgcatgcCTGCAGacacttgggaagatcagatggtgg were used to clone a 4.7 kb genomic region into pCAMBIA1305.1-native-eGFP for generating proLAC23::LAC23-eGFP. 2-mo-old pLAC23-LAC23-EGFP line was used for fluorescence observation.

The stem transverse section preparation was followed as described by Li et al. (Li et al., 2011). Images were observed and captured by using the LEICA DM2500 fluorescence microscope. The autofluorescence of lignin was observed using the third gear (UV excitation light) of the microscope, and EGFP fluorescence was observed using the fifth gear (blue excitation light) of the microscope.

Total RNA was isolated from young leaves, young stems and young roots. cDNA of these three samples were obtained by reverse transcription. Gene-specific forward and reverse primers ( Table S1 ) were used for quantitative RT-PCR to estimate the abundance of gene-specific laccase transcripts in different tissues and stages. The SuperReal PreMix Plus (SYBR Green) kit used was purchased from Tiangen Biotech (Beijing) Co., LTD. The instrument used in the qRT-PCR experiment was Bio-Red.

It has been reported that there are 53 laccase genes in Populus trichocarpa, which are located on 16 chromosomes ( Figure 1 ). Most of the laccase genes are concentrated on five chromosomes (chr 01, chr 06, chr 09, chr 16 and chr 19), and some chromosomes have only one gene (chr 04, chr 07, chr 12, chr 13, chr 14 and scaffold_87).

Segmental and tandem duplications are usually considered to be the main causes of large gene family expansion in plants (Cannon et al., 2004), thus we analyzed the tandem and segmental duplication events of PtrLACs family. The results showed that 22 genes (41.5%) were considered as tandem duplicates, among which two pairs of independent tandem duplication genes were located in chr 10, chr 11 and chr 16. Three groups of three tandem duplicated genes located on chr 05, chr 06 and chr 15. Five tandem duplicated genes located on chr 19. Segmental duplications are regions of 1000 or more base pairs with a DNA sequence identity of 90% or greater that are present in more than one copy (Giannuzzi et al., 2011)^- (Van de Peer et al., 2009). Therefore, the collinearity gene pair can be considered as the fragment duplication genes in Populus laccase gene family. In our study, a total of 22 genes (41.5%) were involved in fragment replication ( Figure 2 ). Based on above results, it could be inferred that tandem duplication and fragment replication play an equal important role in the expansion of Populus laccase gene family.

To better understand the evolutionary constraints acting on Populus laccase gene family, we calculated the non-synonymous nucleotide substitution rate (Ka) and synonymous nucleotide substitution rate (Ks) between both segmental and tandem duplicated PtrLAC gene pairs (So, 2002). The Ka/Ks value of almost all of the gene pairs was less than one (except the tandem repeat pair PtrLAC45/PtrLAC46), which indicated that Populus laccase gene family underwent a strong purifying selective pressure during the evolution.

Existing molecular evolution analysis combined with fossil evidence shows that the molecular evolution rate (E value) of Populus is about 1/6 of that of Arabidopsis (Tuskan and Torr, 2012), while the E value of Arabidopsis is 1.5×10^-8 replacement/synonymous replacement site/year (Koch et al., 1998). Therefore, the E value of Populus in this paper is 1/6 of that of Arabidopsis, that is, 2.5×10^-9 replacement/synonymous replacement site/year. Salicoid duplication event occurs between about 60 and 65 MYA (Million Years Ago) (Tuskan and Torr, 2012). Estimation of the repetition time of collinearity gene pairs in PtrLAC gene family show that about 13.3% (4/30) gene pairs were separated during the WGD event, 20.0% (6/30) of the repeats were replicated earlier than 150 million years ago, and the rest collinear gene pairs were replicated after the WGD events. In addition, tandem repeats were often replicated later than segmental repeats, suggesting that fragment duplication and tandem duplication were the main driving forces of gene family expansion in the early and late stages, respectively ( Table S2, S3 ). These analyses helped us screen out some ancient laccase genes, such as PtrLAC4, PtrLAC41, PtrLAC19, PtrLAC25 and PtrLAC18.

To study the evolutionary relationships among the 75 laccase genes in Populus Arabidopsis, and maize, we built a phylogenetic tree using full-length amino acid sequences. In total, 53 sequences fromP. trichocarpa, 17 sequences from Arabidopsis, 5 sequences from maize, 46 sequences from P. deltoides, and 12 sequences from P. tomentosa were assessed in the phylogenetic tree ( Figure 3 ). The phylogenetic analysis indicated that PtrLACs could be divided into six large groups corresponding to group 1 to 6 in Arabidopsis as defined by Bonnie C (Mccaig et al., 2005). The group 1, 2 and 4 each contains 12 PtrLAC members, group 3 and group 5 each possesses 6 members, and group 6 contains 5 members. In addition, we found that members of group 4 involved in tandem duplication accounted for 66.7% (8/12), all the genes belong to group 5 and group 6 were involved in gene duplication events. These phenomena suggest that gene duplication is responsible for expansion of these 3 subfamilies.

The subcellular prediction shows that most of the laccase proteins are likely to be extracellular, but some have been shown to be localized on the cell membrane or in lysosomes and peroxisomes. A signal peptide analysis shows that most laccases are secretory proteins and have a signal peptide of about 27 amino acids, but the signal peptide of PtrLAC51 was not predicted, indicating that it may be intracellularly localized ( Table S4 ). We found that most laccases in poplar are alkaline, which is similar to the results in most plants, such as A. thaliana and maize.

To better understand the structural features of PtrLACs and gain more insight into the evolution of the laccase family in P. trichocarpa, we analyzed the intron-exons composition of all sequences as shown in Figure 4 . All the sequences in the group 1, group 3 and group 4 contains 5 introns and 6 exons, most sequences in group 2 also have such structure, except PtrLAC8, which contains 4 introns and 5 exons. The sequences in group 5 and most of the sequences in group 6 have four introns and five exons. Notably, in group 6, PtrLAC51 contained only 3 exons and 2 introns, while PtrLAC36 contained 7 exons and 6 introns. Overall, the family members with closer genetic relationships have more similar gene structures ( Figure 4B ).

To further clarify the characteristics of the Populus laccase family, we detected 15 conserved motifs in 53 sequences and analyzed their distribution in combination with the evolutionary tree. As expected, most of the laccase protein sequences in Populus trichocarpa had all 15 motifs. However, in group 6, none of the sequences had motif 15. In addition, PtrLAC51 also lost motif 2 and motif 6, PtrLAC36 doesn’t include motif 6, 9 and 11. In group 4, PtrLAC3 doesn’t have motif 15, and PtrLAC16 just contains six motifs (motif 6, 2, 3, 8, 5 and 1). These differences suggest that family members of the two branches may have lost part of their motifs during evolutionary process, resulting in a new function. ( Figure 4C ).

In order to define the function of each member of the Populus laccase gene family, we use the Pfam database (https://pfam.xfam.org/search) to analyze the domain of each sequence. Results showed that all the members have Cu-oxidase, Cu-oxidase_2, and Cu-oxidase_3 (PF00394, PF07731, and PF07732) these three conserved domains, which associated with the redox reaction of substrates. Beyond that, we also detected 9 domains and some of them have specific functions: YAF2_RYBP was a C-terminal binding motif, which is first found in YAF2 and RYBP proteins and usually forms the RYBP -/YAF2-PRC1 complex (Wang et al., 2015b); ResB domain was related to the synthesis of cytochrome C and is essential for plant growth (Le Brun et al., 2000); NAR2 is a plant protein with a C-terminal transmembrane region, which often works with NRT2 and can transport nitrate at low concentration, this NAR2-NRT2 system plays an important role in the regulation of lateral root growth (Yong et al., 2010); CBiM is a substrate specific component of the cobalt transport complex CBiMNQO, which is related to the synthesis of coenzyme B12 (Santos et al., 2008) ( Figure 5 ). Laccases with these special domains may enable to perform other specific functions in plants, and should be the focus of future research.

To further reveal the phylogenetic mechanisms of Populus laccase gene family, we constructed five comparative syntenic maps of Populus associated with five representative species, including three monocots (Setaria viridis, Oryza sativa, Sorghum bicolor) and two dicots (Arabidopsis thaliana and Glycine max). Results show that, 16 PtrLACs had a collinear relationship with Arabidopsis thaliana, and the number was 5, 7, 6 and 24 for sorghum, rice, Setaria and soybean. As shown in Figure 6 , the number of orthologous gene pairs between Populus and other species was 22 gene pairs (with Arabidopsis), 10 gene pairs (with Setaria), 13 gene pairs (with rice), 7 gene pairs (with sorghum) and 59 gene pairs (with soybean), respectively. Some PtrLAC genes were found to be associated with at least three syntenic gene pairs, such as PtrLAC19 (between Populus and Arabidopsis), PtrLAC25 (between Populus and Setaria), PtrLAC25 and PtrLAC32 (between Populus and Rice). This phenomenon is particularly prominent between Populus and soybean, there are 9 collinearity gene pairs of this type in between (PtrLAC4, PtrLAC50, PtrLAC3, PtrLAC19, PtrLAC21, PtrLAC25, PtrLAC32, PtrLAC37, PtrLAC48).

We also found that, there were 59 collinear gene pairs of laccase genes between Populus and soybean, but only 7 between Populus and Sorghum, which may relate to the phylogenetic relationships between Populus and other five plant species. What is noteworthy is that, many PtrLACs (PtrLAC4, PtrLAC21, PtrLAC33, PtrLAC48) that are collinearly related to Arabidopsis and soybean have not been found in other three monocotyledons, which may indicate that these orthologous pairs formed after the divergence of dicotyledonous and monocotyledonous plants. Some other PtrLACs, however, were identified have collinearity with at least four species (PtrLAC11, PtrLAC25, PtrLAC20, PtrLAC19), suggesting that these genes maybe present before the ancestors differentiated.

In order to understand the expression of PtrLACs in different tissues and stages, we use the RNA-seq data published in Phytozome to draw a heat map and analyzed its FPKM (Fragments Per Kilobase Million) values. we compared the expression patterns of PtrLACs in different developmental processes such as bud set and bud flush, male/female catkin development, leaf expansion and root/stem response to different nitrogen nutrition as shown in Figure 7 .

At different stages of apical bud growth, the expression of laccase genes was more active in mid-spring. During the formation of dormant buds, the expression level of PtrLAC43 is higher in the pre-dormant buds and early dormant buds, which may be related to the protection of meristem in winter. In general, most PtrLAC genes were expressed at low levels in the apical meristem of Populus, which was consistent with the results in Arabidopsis. High expression level of laccase genes was not detected at different stages of stem development, but it was observed in some members of the group 4 (PtrLAC44, PtrLAC45, PtrLAC43, PtrLAC47) that the expression level was higher in late fall f2 and early winter f1, indicating a dynamic response to seasonal variation. It is worth noting that the expression levels of most PtrLACs are high in both nodes and internodes, which may be related to the involvement of laccase in the lignification and thickening of the plant secondary wall.

In male catkin development, PtrLAC5 is consistently highly expressed, most laccase genes had high expression in the early stage of male catkin development and 3 genes (PtrLAC13, 19, 25) showed high expression level in the middle development of male catkin, showed dynamic changes of transcription abundance. In female catkin development, compared with other laccase genes, PtrLAC30, 11, 21, 19 showed consistently higher expression, most of the PtrLACs showed certain dynamic changes with the developmental stage.

In addition, we also analyzed the expression of PtrLAC genes under different nitrogen sources. The results showed that PtrLAC gene was not significantly differentially expressed in stems due to different nitrogen sources, while in roots, it was found that the expression levels of some laccase genes were significantly increased under nitrate nitrogen treatment compared with other nitrogen sources, such as PtrLAC11, PtrLAC13, PtrLAC14, PtrLAC18, PtrLAC19, PtrLAC20, PtrLAC21, PtrLAC26, and PtrLAC41. Accordingly, we speculated that most PtrLAC genes may be involved in the assimilation of nitrate nitrogen in roots ( Figure 7 ).

In order to understand the relative expression level of PtrLACs in different tissues at the early stage of secondary growth, we selected 84K poplar (P.alba×P.glandulosa) trees which were 45-day-old, and analyzed the relative expression levels of PtrLACs in root, stem and leaf by Quantitative Real-time PCR. It turns out that only part of PtrLACs were detected in the early stage of secondary growth, including 20 genes in roots, 17 genes in stems and 12 genes in leaves. We also found some genes with spatial expression specificity in roots and stems, such as PtrLAC5, PtrLAC7, PtrLAC23, PtrLAC12, PtrLAC18 and PtrLAC22, which are currently only expressed in stems. During this period, PtrLAC13, PtrLAC19, PtrLAC21, PtrLAC26, PtrLAC28, PtrLAC31, PtrLAC36, PtrLAC37, PtrLAC42 and PtrLAC52 were only expressed in roots. In addition, the expression levels of PtrLAC genes were relatively close in the same sample, no significant difference was detected. By comparing the relative expression of one PtrLAC in different samples, we found that the expression level of most PtrLAC genes was higher in roots than in leaves and stems. Expression levels of PtrLAC9, PtrLAC10, PtrLAC16, PtrLAC25, PtrLAC33, PtrLAC40 and PtrLAC43 in roots were significantly higher than those in other tissues. ( Figure 8 ).

To further investigate the potential regulatory mechanisms of PtrLACs, we used PlantCARE to predict cis-acting elements within 2.0-kb upstream of each gene sequences, nine cis-acting elements were analyzed and displayed in Figure 8 . We found that the number of cis-acting elements involved in photoresponse was significantly higher than others and some of them have MYB protein binding sites (Baldoni et al., 2015). In view of this interesting finding, we conjectured that the regulation of laccase expression is closely related to photosynthesis. Plastocyanin (Pc) is an important part of the electron transport chain of photosynthesis, which transfer electrons through the redox changes of copper ions in proteins (Höhner et al., 2020). The activity of laccase is also dependent on copper ion (Sun et al., 2019), which leads to the competition effect of copper ions between laccase and Pc. It has been reported that under copper ion stress, the transcriptional regulator SPL7 binds to the copper responsive cis-element (GTAC) in miR397 promoter to promote miR397 synthesis, thereby inhibiting the expression of laccase (Yamasaki et al., 2009). Therefore, we speculated that under copper ion stress, light may inhibit the expression of laccase genes, so that copper ion in plants can meet pc requirements and ensure the progress of photosynthesis. In addition, many PtrLAC photoresponsive cis-elements contain MYB protein binding sites, suggesting that MYB may be participate in regulation of this physiological process as a trans-factor. Besides, PtrLACs also possess cis-acting elements involved in plant hormone response, especially to jasmonic acid and abscisic acid, indicating that laccase may be related to inducing plant chemical defense and bud dormancy. We also detected the cis-elements related to drought induction, low temperature response, defense and stress response in PtrLACs, suggesting a potential stress response under certain conditions. Anyhow, the cis-element analysis illustrated that PtrLACs can respond to different abiotic factors, among which light and copper ion concentration may play a crucial role in the regulation of laccase expression. ( Figure 9 )

In order to explore the relationship between the structure and function of laccase proteins, we selected some representative laccases and simulated its 3D structure by using AlphaFold2 (Cramer, 2021). Older genes may be more closely related to the basic physiological functions in plants. Therefore, we selected PtrLAC25, a gene that had a collinear relationship with all five species in Figure 6 and participated in fragment repetition 200 million years ago, and PtrLAC41, a gene that had a collinear relationship with multiple PtrLACs both in WGD events and 200 million years ago. AtLAC2 has been suggested an involvement in dehydration response (Cai et al., 2006). PtrLAC12 as a member of group 2 is likely to be involved in Populus developmental lignification. In phylogenetic analysis, PtrLAC12 and AtLAC2 had a strong genetic relationship and close genetic distance. Combined with the results of cis-acting elements prediction, we speculated that PtrLAC12 may also played a role in drought stress. Therefore, as a potentially functionally complex laccase, PtrLAC12 was used as a member of the structural simulation. The Cu-oxidase_3 domain of PtrLAC51 was severely lost, which was the most different from other laccase genes in sequence structure, so it was also within our selection range. PtrLAC23 is also a member of group 2 and it was important to note that it did not participate in any repetition events of PtrLAC gene family and its sequence was complete, without motif loss. Meanwhile, it has been proved that miRNA novel-m1190-3p targets PtrLAC23 in the secondary stem of P. trichocarpa (Wang et al., 2021), so we speculated that this gene may be related to wood formation, thus we also predict its structure.

ZmLAC3 is the only plant laccase with known structure and has been confirmed to be implicated in the polymerisation of lignin in maize. Structural comparison with ZmLAC3 protein show that five laccases (PtrLAC12, PtrLAC23, PtrLAC25, PtrLAC41, PtrLAC51) all have low RMSD values (< 3 Å) with ZmLAC3, which means their protein structures are almost the same with ZmLAC3 (Reva et al., 1998) (Pitera, 2014). Populus lignin is mainly composed of G and S type monomers, while H type monomer content is very low (Humphreys and Chapple, 2002). Therefore, coniferyl alcohol (ConA) and sinapyl alcohol (SinA) were selected for molecular docking with selected laccases. Results show that the selected five laccases generally prefer ConA to SinA, and PtrLAC23 has the highest affinity with ConA among them.

The results of Figure 10 indicated that PtrLAC23 (Potri.009G156600) may be involved in the lignification process of Populus. To further prove this prediction, we constructed the fusion protein expression vector of PtrLAC23 and enhanced Green Fluorescent protein (EGFP) initiated by the endogenous promoter: pLAC23-LAC23-EGFP, which was introduced into 84K poplar by Agrobacterium tumefaciens to obtain stable genetic proLAC23::LAC23-EGFP transgenic material. Stem sections of the transgenic plants (growing for 2 months) were observed by LEICA DM2500 fluorescence microscope. Basic principles and flow of the above operations are shown in Figure S1 .

We use the green fluorescence of EGFP and spontaneous fluorescence of lignin to locate PtrLAC23 and lignin, respectively. It turns out that PtrLAC23 is mainly located in xylem, and a small amount is distributed in phloem region. Its distribution position is almost the same as that of lignin, which can be inferred that PtrLAC23 has a certain correlation with the synthesis of lignin. ( Figure 11 ). To further verify the function of PtrLAC23, we determined the relative expression level of PtrLAC23 in poplar stems at different growth stages. Results showed that compared with the 33 days, 56 days and 69 days samples, the expression level of PtrLAC23 was significantly increased in stems of 80 day-old Populus ( Figure 12 ).