The protein sequences of Arabidopsis thaliana laccase family genes (17 members) were retrieved from the Arabidopsis information resource TAIR10 (TAIRv10, http://www.arabidopsis.org/). The sequences were cross-checked with the Pfam database for the presence of three domains of laccases including CuRO_1_LCC_plant (cd13849), CuRO_2_LCC_plant (cd13875), and CuRO_3_LCC_plant (cd13897). Arabidopsis laccase protein sequences were aligned as queries against the Aeluropus Proteome and genome version 1 (Hashemi-Petroudi et al., 2022) using the local BLASTP and TBLASTN program, respectively. The resulting peptide sequences were verified with the blastP tool in NCBI and the Pfam database for the presence of conserved laccase (TIGR03389). In order to reduce the redundancy, % identity of all sequences was computed with the Decrease Redundancy program (https://web.expasy.org/decrease_redundancy), and any pairs of sequences with more than 99% identity were removed from the analysis. The use of the Decrease Redundancy program resulted in reduction from 21 putative AlLACs to 15 defined protein sequences. The naming of each laccase gene in Aeluropus (AlLAC) was based on the closest known orthologue of Arabidopsis.

To construct a phylogenetic tree, AlLACs with their orthologous in Arabidopsis (TAIR10) and rice (IRGSP-1.0) were analyzed using the multiple sequence alignment Muscle tool of MEGA v. 11 software (Tamura et al., 2021) with the default parameters. A phylogenetic tree was constructed using the maximum-likelihood (ML) method with 1000 bootstrap replicates.

The MEME database was used to predict the conserved motifs of AlLAC protein sequences. The number of conserved motifs was adjusted to 10, and other parameters were set as the default. The gene structure of AlLAC genes was constructed using Tbtools (Chen et al., 2020) based on GFF format file for exon and intron location information known from Arabidopsis orthologs. The physical and chemical properties, including number of amino acids, molecular weight (kDa), theoretical pI, grand average of hydropathicity (GRAVY), total number of negatively charged residues (Asp + Glu), total number of positively charged residues (Arg + Lys), and instability index of AlLACs, were calculated using the online ExPASy-ProtParam tool (Gasteiger et al., 2005). The subcellular localization of AlLACs were predicted using WoLF PSORT (Horton et al., 2007) based on default settings.

The Phyre2 server was applied to predict the 3D structure of AlLAC proteins (Kelley et al., 2015). Similar structures were analyzed using the Phyre investigator tool to recognize the pocket site related to the binding region.

PlantCare (Lescot et al., 2002) was used to study cis-regulatory elements in the 1000 bp promoter region. The sequence information was retrieved from Aeluropus genome version 1 (Hashemi-Petroudi et al., 2022). The identified cis-regulatory elements were classified based on their function, and drawn as a graph.

The STRING v11.5 database was used to identify the interactions of AlLAC proteins based on their orthologues in Arabidopsis. The first shell of the network was adjusted to ≤ 20 and the second shell was fixed to no more than 10. Gene ontology (GO) enrichment analysis was used to identify the significant (FDR ≤ 0.05) molecular function, biological process, and cellular component terms presented in the LAC-interaction network using the STRING. Cytoscape v3.8.2 (Shannon et al., 2003) was used to construct the LAC-interaction network.

Seeds of Aeluropus were grown at 25 ± 2°C under 16 hours of light and 8 hours of darkness in the greenhouse at Sari Agricultural Sciences and Natural Resources University. After three weeks, the samples were transferred in groups of three to plastic containers, each containing five liters of Hoagland’s nutrient solution under hydroponic culture (Hoagland and Arnon, 1950). After two months, plants of similar size were selected for exposure to stress conditions. For salinity stress, plants were treated gradually by adding 100 mM salt (NaCl) every 48 hours to a final concentration of 600 mM NaCl. For osmotic stress, plants were treated with 20% PEG 6000 of -0.80 MPa. PEG was added to a plastic container and samples (roots and leaves) were collected at different exposure periods of 0 hour (h), as control sample, 6 hours, 48 hours and one week. For cold stress, the plants were exposed to a temperature of 4°C. Sampling of roots and leaves tissues was done 6 h, 48 h and one week after stress exposure, in leaves three biological replicates. The abscisic acid (ABA) treatment was also done by spraying 100 micromolar hormones on the leaves; two leaf and root tissues were sampled at 3, 6, 24 and 48 hours after applying the treatment, in three biological replicates. Three untreated vases were used as controls.

Total RNA was extracted from leaf and root tissue of all three biological replicates using the Threezol kit (Threezol, Riragene, Iran). The quantity and quality of the RNA samples were measured by spectrophotometry and 1.5% agarose gel electrophoresis, respectively. DNase I treatment (DNase I RNase-free, Thermo Scientific) was used to remove genomic DNA. After combining the RNA of the biological replicates, cDNA synthesis was performed using a kit (Thermo Scientific) according to the company’s instructions, and diluted five times.

The expression levels of the target genes were measured with a Bio-Rad CFX96 machine, using The Maxima SYBR Green/ROX qPCR Master Mix kit (Thermo Scientific) in three technical replicates. The cycling profile started at 95°C for 15 seconds, then 40 cycles of 95°C for 15 seconds and 60°C for 60 seconds. At least one negative control (NTC) was considered for each primer. Melting curve analysis of samples and threshold cycle were calculated with CFX software (Bio-Rad). Normalization of gene expression was done with the geometric mean method using specific reference genes of each tissue. Selecting the appropriate internal control gene is very important when determining gene expression, opting for genes that show the least amount of gene expression changes in different stress conditions (either different tissues or sampled at different times). Specific reference genes for leaf tissue were GTF and U2SnRNP, and specific reference genes for root tissue were PRS3 and EF1a (Hashemi et al., 2016). Five AlLAC genes, AlLAC5, AlLAC12.2, AlLAC14, AlLAC16.1, and AlLAC17.1, were selected based on bioinformatics analyses and primers were designed using AlleleID software (version 7.5, Table S1). Quantitative analysis of the data related to the relative expression level of the studied genes was done using the 2^-△△CT method (Livak and Schmittgen, 2001).

The phylogenetic tree for members of the LAC gene family in Arabidopsis (AtLACs), rice (OsLACs), and A. littoralis (AlLACs) separated LACs into five main groups ( Figure 1 ). Furthermore, AlLACs were more similar to their orthologs in rice, suggesting that a close evolutionary process may have occurred in this gene family in rice and A. littoralis. Moreover, it seems that the diversity of LAC family has more occurred after derivation of monocots and dicots. AlLACs from group I, including AlLAC11.1, AlLAC6, AlLAC16.2, and AlLAC11.2 showed a greater genetic distance than other members, suggesting that these members may have a higher potential for further molecular functional investigation. In addition, the gene structure, conserved motifs, and encoded domain distribution of AlLACs were analyzed ( Figure 2 and Table S3 ). Ten conserved motifs and three encoded domains were identified in AlLACs. Motif 4 was in the Cu-oxidase 1 domain region, while motifs 6, 1, and 3 were predicted in Cu-oxidase 2 domain, and motifs 5 and 2 were located in Cu-oxidase 3 domain ( Figure 2 ). Some of the conserved motifs (motif 9 and motif 10) were identified in the extra-domain region that can be used to identify the AlLACs. Among the various gene structures observed in AlLACs, the maximum exon/intron number was observed in AlLAC11.1.

The predicted 3D structure and binding sites of AlLACs revealed diverse structures among members ( Figure 3 ). Moreover, leucine (L), proline (P), valine (V), phenylalanine (F), glycine (G), and alanine (A) were frequently predicted in the binding site region of AlLACs ( Figure 4 ). The residues detected in the pocket sites of all AlLACs demonstrate the key positions of the important amino acids likely related to their function and interaction points.

The interaction network for AlLAC proteins was drawn based on their orthologs in the model plant, Arabidopsis. Sks11, sks13, and SKU5 proteins, which have oxidoreductase activity, were the most similar to AlLAC11.2, AlLAC16.2, and AlLAC6 proteins, respectively ( Figure 5 ). There were no direct interactions between the laccases themselves. Three interaction groups were observed: the first group contained Sks11 (AT3G13390), sks13 (AT3G13400), and SKU5 (AT4G12420); the second group contained LAC7 and LAC17; and the third group contained LAC5 and LAC6. LAC5 showed the most and strongest interactions compared to the other laccases. Gene ontology (GO) enrichment analysis revealed several biological process terms; iron incorporation into metallo-sulfur cluster, lignin catabolism, regulation of symbiotic processes, plant-type primary cell wall biogenesis, and L-ascorbic acid biosynthesis were significantly linked with the LAC-interaction network ( Table S4 ). In addition, molecular function terms including cysteine desulfurase activity, glucose-6-phosphate isomerase activity, ferrochelatase activity, copper ion binding, oxidoreductase activity, cellulose synthase (udp-forming) activity, and antioxidant activity were significantly associated with the LAC-interaction network ( Table S4 ).

The upstream region of AlLAC genes was analyzed to detect putative cis-regulatory elements. The cis-regulatory elements that are involved with stress, growth and development, phytohormones, and transcription factor (TF) site were observed in the promoter region ( Figure 6 ). The common and unknown function of cis-regulatory elements were frequently observed upstream of AlLACs ( Figure 6A ). Cis-regulatory elements to methyl jasmonate (MeJA) and ABA responsiveness (ABRE) were frequently detected in the promoter region ( Figure 6B ). Moreover, acting elements related to ethylene (ERE), salicylic acid (SA), gibberellic acid (GA), and auxin responsiveness were also observed in the promoter region ( Figure 6B ), as were stress-responsive MYB elements such as as-1, an acting element involved in oxidative stress-responsive, and STRE (stress-controlled transcription factors) elements ( Figure 6C ). Dehydration-responsive element (DRE) and MBS (MYB binding site involved in drought-inducibility) elements were also observed upstream of AlLACs ( Figure 6C ). These results suggest that laccases probably cooperate with factors involved in stress response and are probably present in stress-dependent signaling pathways.

The expression levels of five candidate AlLAC genes (AlLAC5, AlLAC14, AlLAC16.1, AlLAC17.1, and AlLAC12.2) were investigated under ABA application in shoot and root tissues of A. littoralis ( Figure 7 ). The AlLAC genes for real time RT-PCR analysis were selected on their mRNA abundance and salt stress inducibility reported in (Younesi-Melerdi et al., 2020). Based on their expression profile, all selected genes illustrated tissue-specific expression, and they were highly expressed in shoot tissues. Interestingly, all studied AlLACs were sharply upregulated after 48 hours of ABA application in shoot tissues, while AlLAC14 was upregulated at all-time points of ABA application in shoot tissues. Besides, all AlLACs showed a down-regulation after 3 h of ABA application in root tissues. Overall, ABA may indirectly affect the expression of AlLACs in the shoot.

Expression patterns of five selected AlLAC genes in response to abiotic stresses – cold, osmotic (using PEG application), and salt stress – were evaluated in shoot and root tissues ( Figures 8 - 10 ). In response to cold stress, AlLAC5 was upregulated after seven days in shoot tissues, while AlLAC14 was highly induced after 48 hours and seven days in both shoot and root tissues ( Figure 8 ). AlLAC16.1 and AlLAC17.1 were less induced in response to cold stress, and AlLAC12.2 showed upregulation at all studied time points, especially after seven days, in root tissues. In response to salt stress, AlLAC genes showed diverse expression patterns. For instance, AlLAC5, AlLAC17.1, and AlLAC12.2 were upregulated in roots, while AlLAC14 showed a high upregulation in the shoot ( Figure 9 ). Furthermore, AlLAC5, and AlLAC17.1 were upregulated in both root and shoot, suggesting that these genes are directly associated with response to salt stress. In addition, AlLAC16.1 showed a down-regulation in root tissues at all-time points in response to salinity ( Figure 9 ). The expression profile of AlLACs was analyzed under PEG application for induction of osmotic stress. Interestingly, all studied genes showed an upregulation in root tissues, and high expression was recorded after 48 h ( Figure 10 ). In shoot tissues, AlLAC14 was sharply induced after 48 hours in response to drought stress. Overall, AlLAC14 appears to be more expressed in shoot tissues, and AlLACs are more involved in the late response of A. littoralis to abiotic stress.