Genome-wide identification and characterization of laccase family members in eggplant (Solanum melongena L.)

Identification of laccase genes in the eggplant genome

To identify members of the SmLAC gene family in eggplant, the genome assembly and annotation profile of S.melongena 67/3 line (Version 3.0) were downloaded from Sol Genomics Network (SGN, http://solgenomics.net/). A total of 17 laccase sequences from Arabidopsis were obtained according to published research (Turlapati et al., 2011). The reference sequences were used to search predicted proteins in eggplant based on genome files with an E-value cutoff of 1E−5. The candidate SmLAC proteins (SmLACs) were obtained after removing redundant and duplicate protein sequences. To further demonstrate the reliability of the candidate SmLAC proteins, a BLASTP search of the Swiss-Prot databases, was implemented, then all the putative proteins were confirmed by the Pfam (http://pfam.xfam.org/), NCBI CD-search (https://www.ncbi.nlm.nih.gov/Structure/cdd/wrpsb.cgi) and SMART database (http://smart.embl-heidelberg.de/). Finally, laccase members in eggplant with Cu-oxidase, Cu-oxidase_2, and Cu-oxidase_3 domains (PF00394, PF07731, and PF07732) were identified.

The annotation profile (Version 3.0) of eggplant SmLACs localization on chromosomes was obtained from Sol Genomics Network (https://solgenomics.net/). TBtools was used to map the identified SmLACs to the chromosomes based on the annotation information and gene density of the eggplant genome. In order to better identify biological function of SmLACs in eggplant, the physicochemical properties of SmLACs, including MW (molecular weight), pI (theoretical isoelectric point), instability index, aliphatic index and grand average of hydropathicity, were analyzed by online software ExPASy (https://www.expasy.org). CELLO (v. 2.5, subcellular localization predictor, http://cello.life.nctu.edu.tw) was used to predict the subcellular localization of SmLACs.

Phylogenetic analysis of SmLACs

In order to further clarify the classification, evolution and function of laccase members between eggplant and Arabidopsis, multiple sequence alignments were performed between SmLACs and 17 AtLACs from Arabidopsis containing typical characteristics using MEGA-X with the default parameters. An unrooted phylogenetic tree was constructed using the maximum-likelihood (ML) method with the following parameters: poisson correction, partial deletion, and 1,000 bootstrap replicates.

Gene structure, conserved motifs and conserved domains analysis of SmLACs

In order to better understand relationship between SmLACs function and evolution, conserved motifs, conserved domains and gene structure of eggplant laccase members were analyzed. In order to identify conserved motifs, among all the SmLACs, their protein sequences were subjected to online software MEME (http://meme-suite.org/tools/meme) using default parameters, except that the number of motifs was set to 10. Then, conserved domains were searched with SMART (http://smart.embl-heidelberg.de/#). All annotation information of these SmLACs was retrieved from the sol genomics network. Finally, conserved motifs, conserved domains and gene structure of eggplant laccase members were displayed using gene structure viewer in TBtools (Chen et al., 2020).

Cis-elements, gene replication and synteny analysis of SmLACs

For identifying the cis-elements at the promoter regions of SmLAC genes, the 2.0 kb upstream sequence of the start codon (ATG) of each SmLAC gene was extracted from the eggplant genome sequence. A search was performed by using the PlantCARE server (http://bioinformatics.psb.ugent.be/webtools/plantcare/html/). Simultaneously, gene duplications of SmLACs were analyzed and illustrated with TBtools (Chen et al., 2020). Next, according to calculation method reported in a previous study (Liu et al., 2019b), the values of nonsynonymous (Ka) and synonymous (Ks) substitution rates of duplicated genes were calculated. MCScanX in TBtools (Chen et al., 2020) was later used to detect the synteny of SmLAC genes between eggplant and Arabidopsis, eggplant and tomato, with other parameters (filter collinearity in up genome, filter collinearity in down genome, filter gene in small collinearity in block:30).

Expression analysis of SmLACs in different tissues

On the basis of RNA-seq data with SRA accession number SRP078398 reported in a previous study (Barchi et al., 2019), Expression values (FPKM) of SmLACs from 18 different tissues were extracted. It included seven different developmental stages of fruits, four different vegetative tissues, two different reproductive tissues, two different embryo tissues, two different developmental stages of buds and verticillium-inoc. roots 6 hpi. Finally, the expression values of laccase genes were transformed from log₂ (FPKM+1), then the gene-wise normalized and hierarchically clustered heatmap was generated by TBtools software (Chen et al., 2020).

Plant materials, growth conditions and cold stress

Plant materials used in this study were eggplant (Solanum melongena L.) variety Sanyueqie. Seedlings of eggplant were grown in pots containing organic substrate and vermiculite (3:1) mixtures, and then cultured at 25 °C under a photoperiod of 16 h (day) and 8 h (night) in a phytotron for two months. Then all the seedlings were treated at 4 °C for 12 h in a low temperature-programmable incubator, the second and third fully expanded leaves from the top were collected, immediately frozen in liquid nitrogen and then stored at −80 °C until they were used RNA extraction.

RNA extraction, RNA-seq and quantification of gene expression levels under cold stress

The expression of the corresponding genes was also investigated using unpublished RNA-Seq datasets (Sequence Read Archive accessions SRR15036047, SRR15036048, SRR15036049 and SRR15036050). Specifically, total RNA was extracted with Invitrogen TRIzol® Reagent, then the transcriptome libraries were sequenced on the Illumina technology sequencing platform using the paired-end sequencing method (Novogene Co., Ltd, Tianjin, China). Two independent biological replicates for each treatment were performed for the transcriptome sequencing. Based on the subsequent published reference genome sequence from eggplant (Barchi et al., 2019), we used the TBtools (Chen et al., 2020) software to complete the analysis from the raw data to the gene expression level. In RNA-seq analysis, the gene expression levels can be estimated by counting the sequencing sequences (reads) located in genomic regions or gene exon regions. In addition to being directly proportional to the gene expression levels, the reads counts are also positively related to the length of the gene and the sequencing depth. In order to make the estimated gene expression levels of different genes and different experiments comparable, fragments per kb (kilo-base) per million fragments (FPKM) considers the effect of sequencing depth and gene length on fragment counts and is currently the most commonly used method for quantifying expression level of genes (Trapnell et al., 2010). The expression values of laccase genes were transformed from log₂ (FPKM+1), then the gene-wise normalized and a hierarchically clustered heatmap was generated by TBtools software (Chen et al., 2020).

Plant materials, growth conditions and stress treatments

In this study, we used eggplant variety “Sanyueqie” for plant materials. Seedlings of eggplant were grown in pots containing organic substrate and vermiculite (3:1) mixtures, and then cultured at 25 °C under a photoperiod of 16 h (day) and 8 h (night) in a phytotron for two months. Then all the seedlings were treated according to the following methods: (1) for cold stress, the plants were maintained at 4 °C in a low temperature-programmable incubator for 0, 1, 6, 12, 24 and 36 h; (2) for drought stress, the plants were irrigated with 20% PEG (cultured at 25 °C under a 16 h light/8 h dark cycle) for 0, 1, 6, 12, 24 and 36 h; and (3) for salinity stress, the plants were dipped in 200 mM NaCl and maintained at 25 °C with a 16 h light/8 h dark cycle for 0, 1, 6, 12, 24 and 36 h. At different time points, we collected the second and third fully expanded leaves from the top, immediately put them in liquid nitrogen, and then stored these samples at −80 °C for RNA extraction. In all abiotic stress treatments, we performed three biological replications independently.

RNA extraction, and expression analysis

We extracted total RNA from leaves of treated and control group eggplant seedlings using RNAprep Pure Plant Kit (TIANGEN) according to the manufacturer’s specifications, and then performed quantification with a two-step reaction process: reverse transcription (RT) and qPCR. Reverse transcription was performed in 2,720 Thermal Cycler (Applied Biosystems™). qPCR reaction was performed on a CFX96 Real-Time PCR (Bio-Rad, Hercules, CA, USA). At the end of the PCR cycles, we performed analysis of the melting curve to validate the specific generation of the expected PCR product. Adenine phosphoribosyl transferase (APRT, accession JX448345, (Gantasala et al., 2013)) from eggplant was used as the internal control. The relevant primers are listed in Table S1. All reactions were run in triplicate. The relative expression values of SmLAC genes were calculated using the 2^−ΔΔCt method.

Statistical analysis

Data analyses were conducted using the SPSS version 25.0 (SPSS, Chicago, IL, USA) statistical software. For all analysis, the level of significance was set at p < 0.05. Column charts were drawn using Excel 2019 software. Sample variability was given as standard deviation (SD).

Identification of laccase gene members in eggplant genome

To identify the laccase family members in eggplant, we downloaded the genome assembly and annotation profile of eggplant (Version 3.0) from Sol Genomics Network. Subsequently, seventeen laccase sequences from Arabidopsis as reference sequences were used to search predicted proteins in eggplant based on genome files with an E-value cutoff of 1E−5. After removing redundant peptide sequences with only one or two typical domains, 48 laccase family members in eggplant were identified, for specific coding region sequences and peptide sequences, please see Tables S2 and S3. These members were then numbered according to their location on the chromosomes and designated SmLAC1 through SmLAC48 (Table S4). Moreover, 42 SmLAC gene members were widely distributed on 12 chromosomes and six SmLAC genes were associated with scaffolds (Fig. 1). Notably, chromosome 5 contained the largest number of laccase genes, 11 genes (SmLAC14-SmLAC24) were spread out along chromosome 5. On the contrary, chromosomes 8, 10 and 12 contained only one gene respectively. In addition, six genes (SmLAC43-SmLAC48) were not mapped onto any chromosome. Then physico-chemical properties of 48 SmLAC family members were further analyzed. The length of SmLAC proteins ranged from 494 to 835 amino acids (581 amino acids on average), with molecular weight ranging from 55.32 to 93.41 kDa (64.91 kDa on average), and isoelectric points (pI) varying from 5.38 to 9.60 (8.50 on average). In addition, most of SmLAC proteins instability index values were below 40.0, whereas instability index values of only four SmLAC proteins were above 40, which indicated that most of the eggplant laccase proteins belonged to stable proteins. Furthermore, the grand average of hydropathicity of all eggplant laccase proteins was below 0.0, it implied that all eggplant laccase proteins were hydrophilic proteins.

Analysis of the phylogeny of eggplant laccase family

To further study the phylogenetic relationships of laccases between eggplant and Arabidopsis, an unrooted phylogenetic tree was constructed using the maximum-likelihood (ML) method. On the basis of six different groups previously reported for Arabidopsis, 25 eggplant laccase members were divided into six groups, but the laccase members from eggplant were not equally distributed in each group, eight family members were observed belonging to group III, seven members belonging to group II, four members belonging to group I, and one member each belonging to group V and group VI (Fig. 2). Moreover, eggplant laccase members in group IV were classified into two subgroups: subgroup IVa (SmLAC14) and subgroup IVb (SmLAC10, SmLAC11 and SmLAC21). In addition, the remaining 23 SmLAC members found no hits among the AtLACs, they were temporarily classified as group VII and group VIII (Fig. 2).

Analysis of gene structure, motif compositions and conserved domains of eggplant laccase family members

To show the structural diversity of eggplant laccase family members, exon/intron organizations, conserved motifs and domains were analyzed based on the phylogenetic tree of all eggplant laccase members alignments (Fig. 3). Our results showed that almost all the closest genes encoding those SmLACs on the phylogenetic tree showed remarkably analogous gene structures. Moreover, nearly 40% of SmLAC gene members included six exons, and the number of exons varied from 1 to 13 (0 to 12 introns) (Fig. 3). However, there was still a small proportion of the gene clades that exhibited different intron/exon organizations. For instance, SmLAC44 had eight exons, whereas its nearby paralogous genes, SmLAC5, SmLAC28, SmLAC38, SmLAC43 and SmLAC45, typically had six exons even though their evolutionary relationships reached 87–100% bootstrap values, respectively. Conserved motifs of the 48 SmLACs were analyzed through the MEME program. It was discovered that six conserved motifs including motif 1, motif 2, motif 3, motif 4, motif 5 and motif 10 were found in all 48 putative SmLAC members. Moreover, the motif compositions in the same group possessed similar structures and organizations. SmLACs from group I, group II (except for SmLAC44 and SmLAC48), group III (except for SmLAC30), group IV, group V and group VI contained completely the same motif composition, which suggested the possibility of functional redundancies among these members (Fig. 3). The conserved motifs might play critical roles in mediating functions of SmLACs. Furthermore, the 48 eggplant laccase members contained three conserved domains of Cu-oxidase, Cu-oxidase_2 and Cu-oxidase_3 identified by SMART (Fig. 3).

Analysis of the cis-regulatory elements predicted in the promoters of SmLACs

To find possible regulatory patterns of the SmLACs, we extracted the 2.0 kb upstream sequence of the start codon (ATG) of each SmLAC gene (Table S5), and examined the putative cis-regulatory elements in the promoter regions of eggplant laccase genes by PlantCARE. All SmLAC promoter regions were predicted to contain many light responsive elements (Table 1), such as G-box, GT1-motif and so on, which indicated important roles in eggplant morphogenesis. In addition, abscisic acid (ABA) and anaerobic induction (AAI) relative elements were more frequently identified in the promoter regions of laccases (Table 1). Simultaneously, other important representative cis-regulatory elements were also predicted. Among them, in the promoters of the SmLACs, a great number of cis-elements associated with hormone response were identified in many laccase genes (Table 1). For example, 39 out of 48 laccases included ABA-responsive cis-regulatory elements, 27 laccases contained gibberellins(GA) response elements, 26 laccases had salicylic acid (SA) responsive elements, and 21 laccases contained cis-regulatory elements in relation with auxin and methyl jasmonate (Me-JA) respectively (Table 1).

In the same way, a great number of stress-responsive cis-elements, in relation with anaerobic induction, defense and stress, drought, low temperature and wound, were also detected in the promoters of the SmLACs (Table 1). For instance, stress-responsive cis-regulatory elements, responding to AAI, were detected most in the promoters of 35 SmLACs, while stress-responsive cis-elements, related to wound, were predicted least in the promoters of two SmLACs. At the same time, tissue specific cis-regulatory elements, related to the development of meristem, endosperm, seed and zein metabolism, were also predicted within the promoters of many laccase genes from eggplant (Table 1). Among them, cis-regulatory elements, related to development of meristem, were found most in the promoters of 17 SmLACs, but cis-regulatory elements, related to the development of seed, were discovered only in the promoter of SmLAC14 (Table 1). Furthermore, MYB binding sites, involved in flavonoid biosynthetic genes regulation, were also found in the promoters of five SmLACs, including SmLAC19, SmLAC32, SmLAC35, SmLAC40 and SmLAC48, which indicated that these laccases may be related to the process of flavonoid biosynthesis in eggplant.

Analysis of gene duplication and collinearity in eggplant laccase family

Gene duplication, including tandem repeats, segmental repeats and interspersed repeats, is considered to be the main driving force under the process of genome evolution. Interestingly, we discovered that two pairs of eggplant laccase genes (SmLAC15/SmLAC16 and SmLAC36/SmLAC37) contained the same protein coding sequences and physico-chemical properties, so they were temporarily identified as interspersed repeats. In order to gain insight into the eggplant laccase gene duplication pattern, we further analyzed the tandem and segmental duplication events of eggplant laccase genes by TBtools. In this study, we discovered that eight pairs of 48 laccase genes belonged to tandem repeats (Fig. 1, Table S6). However, there were no segmental repeats in the eggplant laccase gene family. Then we obtained Ka and Ks values of each duplicated SmLAC gene pair, using simple Ka/Ks calculator in Tbtools. It was discovered that the Ks values of tandem repeats were between 0.043 and 2.038. Moreover, 75% (6/8) of Ks values were less than 1. In addition, we calculated Ka/Ks values of all SmLAC gene pairs, and discovered that all the Ka/Ks values were less than 1, which revealed these genes had evolved in the selection of strong purification. Furthermore, we calculated the approximate dates of SmLAC gene duplication events, in accordance with the method reported in a previous study(Liu et al., 2019b). The duplication events of SmLAC genes happened from 8.27 Mya (Ks = 0.043) to 391.86 Mya (Ks = 2.04), with an average of 113.05 Mya (Ks = 0.59). Detailed information in regard to duplication gene pairs, duplication type, Ka, Ks, Ka/Ks and approximate duplication date (Mya) of all identified duplicated SmLAC genes were listed in Table S6.

In order to detect the homology of eggplant laccase genes, we analyzed the collinearity, between eggplant and tomato, eggplant and Arabidopsis thaliana, eggplant and rice, with Dual Systeny Plot for MCscanX in TBtools. As shown in Fig. 4 and Table S7, there were twenty-four pairs of collinearity between twenty-one SmLAC genes and twenty-four SlLAC genes, six pairs of collinearity between six SmLAC genes and six AtLAC genes, but only two pairs of collinearity between one SmLAC gene and two OsLAC genes. Moreover, twenty-one SmLAC genes were collinear with LAC genes from Arabidopsis, tomato or rice. Unfortunately, none of SmLAC genes were collinear with LAC genes in other three species.

Analysis of expression patterns of SmLACs in different tissues

To determine the function of SmLAC genes in eggplant in different developmental stages, we analyzed the expression patterns of forty-eight SmLAC genes in various vegetative and reproductive stages of eggplant, based on RNA-seq data reported in previous study (Barchi et al., 2019). Overall, according to the data, we classified the expression patterns of SmLAC genes into three types (no-expression, constitutive expression and specific expression) (Fig. 5 and Table S8). Among them, gene SmLAC19 was not expressed in all of the examined developmental stages. By contrast, seventeen SmLAC genes were ubiquitously expressed in all eggplant organs/tissues, but differentially expressed in different tissues. Among them, SmLAC6, SmLAC8, SmLAC34 and SmLAC47 genes showed high expression levels in all tissues/organs, which indicated that these genes may involve in metabolism, plant growth and development. But thirteen SmLAC genes were discovered at high level in many tissues/organs, but low FPKM values were found in some organs. For instance, gene SmLAC13 was accumulated at a relatively higher level in roots, fruits and flowers, yet at the lowest level in senescent leaves. In addition, the other thirty-one SmLAC genes showed relatively high level in some organs, but a relatively low level or no-expression in many tissues/organs. For example, SmLAC2, SmLAC24, SmLAC39 and SmLAC40 were detected at relatively high levels in flowers and open buds, while at a lower to no-expression level in other tissues (Fig. 5 and Table S8).

Expression analysis of SmLAC genes under cold stress

According to our unpublished RNA-seq data under cold stress, a heat map was displayed on the basis of log₂ transformed (FPKM+1) values. Compared with normal growth conditions, twenty-one SmLAC genes were up-regulated, which indicated that 43.75% of laccase members were induced by cold stress (Fig. 6 and Table S9). On the contrary, twenty-two SmLAC genes were down-regulated under cold stress for 12 h, which suggested that 43.75% of laccase members exhibited suppressed expression. In addition, five SmLAC genes (SmLAC2, SmLAC11, SmLAC16, SmLAC19 and SmLAC24) were detected with no expression values in both conditions, which was basically consistent with the above published RNA-seq data regarding the expression of the laccase family gene in the leaf tissue. Unfortunately, no significant differentially expressed genes have been found based on RNA-seq data under cold stress. From this result, it was suggested that cold stress could not significantly affect the expression of laccase family genes.

Expression patterns of SmLAC genes in response to cold, drought and salt tress

On the basis of the potential cis-regulatory elements of SmLACs promoters, it was indicated that SmLAC genes may respond to various stresses. Therefore, six SmLAC genes (SmLAC12, SmLAC15, SmLAC17, SmLAC23 and SmLAC41) were selected to detect their expression patterns under various stress conditions. As shown in Fig. 7, for cold stress, the expression levels of SmLAC15, SmLAC17, SmLAC23 and SmLAC41 were significantly up-regulated on the whole. Among them, the expression levels of SmLAC15, SmLAC23 and SmLAC41 were upregulated most significantly at 12 h (27.7, 2.3 and 32.8 folds as much as that of the control respectively), while the relative expression value of gene SmLAC17 showed the most significant level at 24 h. On the contrary, gene SmLAC12 was significantly down-regulated at different time points, except for 24 h. Whereas, gene SmLAC26 showed significantly down-regulated expression at 1, 6, and 36 h, but increased sharply to a significant up-regulated expression level at 12 h. Under drought treatment, the expression values of SmLAC15, SmLAC17 and SmLAC26 exhibited up-regulated change within 6 h, and showed peak expression at 6 h, then declined. With the increase of processing time, SmLAC26 was still significantly up-regulated at 12, 24 and 36 h, but SmLAC17 reached significantly up-regulated expression at 24 and 36 h again, and SmLAC15 only at the time point of 36 h. In addition, the expression values of SmLAC12, SmLAC23 and SmLAC41 exhibited similar change tendencies (first declined, increased and then declined). Overall, compared with the control, gene SmLAC23 emerged significantly down-regulated levels at different time points, but SmLAC12 and SmLAC41 possessed relatively complex expression changes, first displayed continuously down-regulated expression within 12 h, then SmLAC12 attained a significant down-regulation level again at 36 h, while SmLAC41 showed a significant up-regulation level at 24 h. As for salt treatment, both of SmLAC12 and SmLAC17 displayed significantly down-regulated levels at different time points, but SmLAC26 showed significantly inhibited expression level only at 6, 12 and 24 h. On the contrary, SmLAC15 and SmLAC41 acquired significantly induced expression levels at 12 and 36 h, while SmLAC23 obtained significantly induced expression levels at 6 and 36 h, but none of the three genes reached the levels of significant difference at other time points. From the above results, we can learn that the selected SmLAC genes except for SmLAC12 were significantly induced by one or multiple stresses.