Genome-wide characterization of the xyloglucan endotransglucosylase/hydrolase gene family in Solanum lycopersicum L. and gene expression analysis in response to arbuscular mycorrhizal symbiosis

Identification of XTH family members in Solanum lycopersicum

All gene and protein sequence information was retrieved by searching the Phytozome v13 database (http://www.phytozome.net) and the Solanaceae crops genome database (https://solgenomics.net/). To identify all SlXTH proteins, the BLASTP algorithm using the SlXTH14 amino acid sequence was employed to search all potential XTH proteins in the Solanum lycopersicum genome. SlXTH14 was selected based on transcriptomic analysis as previously reported by Cervantes-Gámez et al. (2016), in which several cell wall biogenesis-related genes were differentially expressed in tomato leaves in response to AM symbiosis. The Hidden Markov Model (HMM, https://www.ebi.ac.uk/Tools/hmmer/) was used to search the profiles of the SlXTH protein domains PF00722 and PF06955, as previously reported by Wang et al. (2018). The online program SMART (http://smart.embl-heidelberg.de/) was used to identify the conserved domain of candidate SlXTHs, and only the proteins containing both domains PF00722 and PF06955 were kept for further analysis. The chromosome coordinates of each SlXTH genomic sequence (File S1), as well as their coding (File S2), transcript (File S3), and protein (File S3) sequences were obtained from the Phytozome v13 database. Physicochemical parameters for each protein, including predicted molecular weight and isoelectric point (PI), were obtained using tools available at the ExPASy bioinformatics resource portal (https://www.expasy.org/). The subcellular localizations were predicted with ProtComp 9.0 (http://linux1.softberry.com). The SignalP 5.0 server (https://services.healthtech.dtu.dk/service.php?SignalP-5.0) was used to predict the presence of signal peptides. Finally, SlXTH genes were nominated as previously reported (Saladié et al., 2006), and new sequences were named according to the following numbers.

Gene structure and motif analysis

Genomic and complete coding DNA (CDS) corresponding to each identified SlXTH gene were analyzed for exon-intron distribution. The Gene Structure Display Server (GSDS 2.0) (http://gsds.gao-lab.org/) was employed to obtain graphical representation of the exon-intron organization by comparing the CDS sequences of the SlXTH genes to the corresponding genomic DNA sequences (Hu et al., 2015). Protein structural motif analysis was performed using the MEME program (https://meme-suite.org/meme/) to predict conserved motifs (10 maximum motifs) in SlXTH proteins as previously reported (Bailey et al., 2009). The consensus sequence was analyzed to identify the conserved catalytic motif (DEIDFEFLG) of SlXTH proteins, and the web logo was illustrated using the MEME tool.

Structurally based sequence alignment and structural prediction of SlXTH proteins

The bioinformatics online tool ESPript (https://espript.ibcp.fr/ESPript/ESPript/) was used to predict the secondary structures in the SlXTH protein sequences and the secondary elements. SlXTH sequences were aligned using ClustalW with default settings to identify shared structural features of SlXTHs, and the PDB databank (https://www.rcsb.org/) was used to locate the XTH crystal protein structure (PDB id: 2UWA; PDB id:1UN1) as previously reported (Johansson et al., 2004; Baumann et al., 2007). Three-dimensional (3D) structures predicting models of SlXTH proteins were constructed based on the oligomeric state, the maximized percentage identity, ligands, the model quality estimation (QMEAN) and the global quality estimation score (GMQE), using the SWISS-MODEL template library (https://swissmodel.expasy.org/) (Biasini et al., 2014).

In silico chromosomal mapping, gene duplication and Ka/Ks estimation

The chromosomal location of each SlXTH was obtained from the Phytozome v13 database. The physical location and relative distances of SlXTH genes were schematically represented on their respective tomato chromosome using the online server MG2C (http://mg2c.iask.in/mg2c_v2.0/). To analyze gene duplication events, tandem and segmental duplications were considered. A gene pair on the same chromosome located five or fewer gene loci apart and showing more than 90% sequence similarity was considered a tandem duplication, whereas sister gene pairs located on different chromosomes were considered segmental duplication events. To estimate the selective pressure and divergence time of SlXTH genes, amino acid and coding sequences from segmental gene pair duplications were analyzed using the Toolkit for Biologists Tools (TBtools) software (https://github.com/CJ-Chen/TBtools) to determine the Ka (non-synonymous), Ks (synonymous), and Ka/Ks ratio parameters (Chen et al., 2020). The approximate time (T) duplication event was estimated using the (T) = Ks/ 2 λ ×10⁻⁶ million years ago (Mya) for each gene pair, where λ = 1. 5 × 10⁻⁸ substitutions per site per year for dicot plants (Koch, Haubold & Mitchell-Olds, 2000).

Gene ontology (GO) annotation

Gene ontology annotation analysis of SlXTH genes was conducted using the Blast2GO software (https://www.blast2go.com/) (Conesa & Götz, 2008). Amino acid sequences of each SlXTH gene were uploaded to the program and the biological process (BP), cellular compartments (CC), and molecular functions (MF) were determined. In addition, a Blast2GO analysis was performed to BLASTp search, InterPro Scan, mapping and annotation with default settings.

Analysis of cis-acting regulatory elements from the SlXTH genes

To predict the cis-acting elements in the promoter of SlXTH genes, 2.0 kb upstream of the initiation codon (ATG) of each SlXTH gene were extracted from the Phytozome v13 database. The upstream sequences were submitted to the online PlantCare (http://bioinformatics.psb.ugent.be/webtools/plantcare/html/) database for the prediction (Lescot, 2002) and visualized using GraphPad Prism 6 software.

Phylogenetic analysis of SlXTH proteins

For insight into the evolutionary relationship among different XTH gene family members, we performed a multiple sequence alignment of the full-length XTH protein sequences from other solanaceous plants such as N. tabacum, S. tuberosum, Petunia axillaris and the model plant A. thaliana using ClustalW with default parameters. We analyzed the results with MEGA X (http://www.megasoftware.net) (Kumar et al., 2018). The phylogenetic tree was constructed based on the neighbor-joining algorithm with 1000 bootstrap replications, and was visualized using the iTOL online tool (https://itol.embl.de/) (Letunic & Bork, 2016).

Gene expression analysis

RNA-seq expression data for S. lycopersicum tissues including leaves, roots, buds, and flowers were downloaded from the GEO database at NCBI (http://www.ncbi.nlm. nih.gov/geo/) and the Solanaceae crops genome database (SRA049915: accession numbers SRX118613, SRX118614, SRX118615, SRX118616) (FDR less than 3%; q-value threshold <0.03) (The Tomato Genome Consortium, 2012). The XTH expression data were estimated using the expressed ‘reads per kilobase of exon per million fragments mapped’ (RPKM) value. RPKM values for different tissue were subjected to hierarchical clustering analysis with TBtools software (https://github.com/CJ-Chen/TBtools). Finally, the data were normalized to examine differences in the expression of the same gene in different samples and represented as a heatmap with TBtools.

Plant material and growth conditions

The arbuscular mycorrhizal fungi Rhizophagus irregularis was provided by the CIIDIR-SINALOA at the Instituto Politécnico Nacional in Sinaloa, Mexico. The inoculum was grown according to previously reported methods (Bécard & Fortin, 1988).

S. lycopersicum (var. Missouri) seeds were surface-sterilized. Tomato seeds were planted in germination trays with a mixture of sterilized vermiculite and sand (3:1 v/v) and maintained at 25 °C. Four-week-old tomato plants were transplanted individually to pots (1 L) with the same substrate. At this time, tomato plants were inoculated with 500 spores of R. irregularis (M+ treatment). AM spores were prepared from an axenic carrot root culture colonized with R. irregularis and extracted as previously described by Cervantes-Gámez et al. (2016). Control samples consisted of mock-inoculated plants (i.e., non-colonized, M- treatment) with the last rinse of the spore inoculum wash. All plants were watered once per week with distilled water and twice per week with 30 mL of half-strength Hoagland nutrient solution with 50 µM KH₂PO₄ as the final phosphate concentration to favor the mycorrhizal colonization (Hoagland & Arnon, 1950). One-half of the root system and whole leaves from each M+ and M- tomato plant were harvested four weeks after R. irregularis inoculation. The collected plant material was immediately frozen in liquid nitrogen and stored at −80 °C for subsequent RNA extraction. Five biological replicates were employed per treatment and two independent experiments were performed.

The other half of each plant root system was fixed in 50% ethanol, clarified in 20% KOH, neutralized in 0.1 M HCl, and stained in 0.05% trypan blue in lactoglycerol (Phillips & Hayman, 1970). Roots were maintained in lactoglycerol 1:1:1 (water/lactic acid/glycerol) and observed by light microscopy (BOECO Germany, BM-180). Mycorrhizal colonization was confirmed as previously reported by Mendoza-Soto et al. (2022).

RNA extraction, primer design, and qRT-PCR analysis

Nine out of the 37 identified SlXTH genes were selected to experimentally determine their expression in leaves and root tissues of mycorrhizal colonized and non-colonized plants based on the occurrence of defense-related regulatory elements within their promoter sequences.

Total RNA was isolated from leaves and roots of non-colonized (M-) and colonized (M+) plants using TRIzol reagent (Invitrogen, Carlsbad, CA), following the manufacturer’s protocol. The complementary DNA synthesis was performed as previously reported (Cervantes-Gámez et al., 2016). The 3′  untranslated regions (UTR) of each gene were used to design qPCR primers for gene specificity. The primers used are listed in Table S1. Melting temperature (Tm) and GC content were calculated using Oligo Calc (http://biotools.nubic.northwestern.edu/OligoCalc.html). qRT-PCR was performed using SYBR Green (QIAGEN, USA) and quantified on a Rotor-Gene Q (QIAGEN, USA) real-time PCR thermal cycler. qRT-PCR was programmed for 40 cycles, denaturing at 95 °C for 15 s, annealing at 58 °C for 30 s, and extension at 72 °C for 30 s. Amplification of a single PCR product was verified by thermal gradient PCR and melting curve qRT-PCR analysis. The elongation factor 1-α (SlEF1- α) gene was used for normalization. The relative expression of SlXTH genes was calculated by the 2^−ΔCT method (Livak & Schmittgen, 2001). Five biological replicates for each condition (non-colonized and colonized plants) were evaluated, and two independent experiments were performed with similar results. Data from one of the experiments are shown.

Data analysis

For the relative expression of each SlXTH gene, the paired Student’s t-test was used to evaluate the significance of differences between non-colonized (M-) and colonized (M+) tomato plants. All data were checked for normal distributions (Shapiro–Wilk’s test) before statistical analyses, which were performed using the scientific data analysis and graphing software SigmaPlot for Windows, version 11.0.

Identification and characterization of the SlXTH gene family

A comprehensive genome-wide screening of the tomato database was executed to identify all SlXTH genes. As a result, 37 SlXTH genes were identified, including some novel members of the family. All SlXTH genes identified within the tomato genome showed the conserved PF00722 (glycosyl hydrolases family 16) and PF06955 (xyloglucan endotransglucosylase C-terminus) domains by using Pfam analysis (Fig. S1), which verifies and validates the sequence search results. The 37 SlXTH genes were named SlXTH1 to SlXTH37 based on a previously reported work in which SlXTH1 to SlXTH25 were already identified in S. lycopersicum (Saladié et al., 2006). The characteristics of each sequence, including gene ID number and length of the genomic, transcript, coding DNA (CDS) and amino acid sequences, as well as molecular weight (MW), isoelectric point (PI), and chromosome coordinates, are summarized in Table 1 and Table S2. The length of SlXTH proteins ranged from 274 (SlXTH31) to 372 (SlXTH26) amino acids, with the predicted CDS ranging from 825 to 1,119 bp and the calculated MW varying between 0.69 and 1.63 kDa. SlXTH26 was the most significant XTH protein (Table 1). The theoretical PI values of SlXTH ranged from 4.85 (SlXTH27) to 9.51 (SlXTH14) due to the differences in ionic strength and pH in the amino acids present in these proteins (Table 1).

Subcellular localization prediction revealed that most of the SlXTH proteins (32 out of 37 SlXTHs) were located on the plasma membrane. In contrast, SlXTH5, SlXTH6, SlXTH14, SlXTH26, and SlXTH36 were predicted to localize in the extracellular region (Table S2). In addition, the signal peptide prediction indicated that all SlXTH proteins contain signal peptide sequences exceptSlXTH13, SlXTH18 and SlXTH22 (Table S2).

Gene structure and conserved motif analysis of SlXTH proteins

To investigate the structural diversity of the SlXTH genes, exon and intron structures of the 37 SlXTH genes were determined by aligning their CDS and genomic sequences using the GSDS server. We also constructed a phylogenetic tree using full-length deduced amino acid sequences of the SlXTH genes, presented with the exon and intron distribution in Fig. 1. The phylogenetic tree shows that the SlXTH genes are divided into two major subfamilies: subfamily I/II and III. Subfamily I/II has 30 gene members, while subfamily III has seven. Structural analysis of SlXTH genes showed that each subfamily’s most closely related genes share similar exon and intron numbers. For example, members from subfamily I/II mostly contain three introns and four exons distributions, exceptSlXTH13, which presents five exons in its coding region. All members of subfamily III also contain three introns and four exons in their coding region. On the other hand, the presence or absence of 5′  and 3′  UTR was not exclusively associated with either of the two subfamilies. For example, SlXTH6 and SlXTH8 from subfamily III, and SlXTH30, SlXTH31, SlXTH33, SlXTH22, SlXTH18, SlXTH29, SlXTH11, SlXTH9, and SlXTH13 from subfamily I/II do not have any 5′  or 3′  UTRs (Fig. 1).

To further characterize the SlXTH family, MEME motif detection software was used to predict potentially conserved motifs. A total of ten conserved motifs with lengths of ten amino acids were identified. Motif compositions differed in members from the two subfamilies (Fig. 1 and Fig. S2). For example, all SlXTH members from subfamilies I/II and III showed the presence of the ten highly conserved motifs, except for SlXTH36 and SlXTH31 from subfamily I/II, which only have seven and six conserved motifs, respectively (Fig. S2A). Furthermore, multiple sequence alignment of all SlXTH proteins revealed the conserved amino acid motif DEIDFEFLG, which is responsible for the catalytic activity as well as being the most characteristic motif of this family (Fig. S2B), indicating that this conserved core motif is an essential for XTH proteins, and suggesting that all of these proteins have a similar function. The majority of the SlXTH proteins within the same subfamily showed identical gene structure and motif compositions, consistent with the phylogenetic analysis of the whole XTH gene family.

Structural prediction of SlXTH proteins

The alignments of SlXTHs with a xyloglucan endotransglycosylase crystal protein structure (PDB id: 2UWA and PDB id:1UN1) were used to predict the secondary structures of the SlXTH proteins with ESPript (Figs. S3 and S4). All SlXTH protein members in subfamily I/II and subfamily III had similar structures to the reference crystal protein structure. Twenty-eight subfamily I/II members showed a conserved position of the N-glycosylation site at amino acid 99. However, this was not found in two SlXTH proteins (SlXTH31 and SlXTH36) (Fig. S3, see label *). Amino acid 116 was also conserved in all members of subfamily III (Fig. S4, see label *). The active site (ExDxE) containing the residues responsible for catalytic activity was highly conserved in all SlXTH family members (Fig. S3 and S4, see label AS). In addition, all members possess the XET/XEH C-terminal extension, a characteristic fingerprint among XTHs from other plant species.

A tertiary (3D) protein model of SlXTHs might help to understand XTH enzyme’s structure and possible mode of action (Table S3). Most SlXTHs from the same subfamily (I/II and III) showed similar 3D structures with percentage identities between 36.54 and 77.32%, indicating a reliable structure prediction. In addition, essential ligands were predicted based on their chemical identity. For example, ligands for β-D-glucopyranose, α-D-xylopyranose and β-D-galactopyranose were identified in 34 SlXTH proteins but not in SlXTH8, SlXTH14, and SlXTH21, in which no ligands sites were detected (Table S3).

Chromosome mapping and gene duplication analysis of SlXTH genes

The chromosome coordinates of all SlXTH genes were obtained from the Phytozome v13 database (Table S1), and their chromosomal locations were mapped using the online server MG2C (http://mg2c.iask.in/mg2c_v2.0/) (Fig. 2). SlXTH genes were heterogeneously distributed among all chromosomes across the tomato genome. The most significant number of SlXTH genes were located on chromosomes 12, 3, and 7, with five, six, and seven SlXTH genes, respectively. In contrast, chromosomes 4, 6, 8, and 10 had only one SlXTH gene. The other chromosomes contained between two and four SlXTH genes (Fig. 2). Finally, no SIXTH gene was found on chromosome 0.

Tandem and segmental duplications reveal information about the expansion of new gene family members and evolutionary functions in plants (Ganko, Meyers & Vision, 2007). Tandem duplications during SlXTH evolution were investigated using the Smith-Waterman algorithm alignment. Two SlXTH gene pairs (SlXTH3/SlXTH37 and SlXTH24/SlXTH35) were confirmed to be tandem duplicated since sequence similarity was higher than 90% (Table S4). Both SlXTH gene pairs are on chromosome 3 (Fig. 2, see label *). A total of 12 segmental duplication events were identified based on phylogenetic analysis, which includes nine sister pairs (SlXTH36/SlXTH3, SlXTH15/SlXTH27, SlXTH4/SlXTH1, SlXTH16/SlXTH28, SlXTH18/SlXTH29, SlXTH10/SlXTH11, SlXTH9/SlXTH17, SlXTH2/SlXTH19, SlXTH24/SlXTH37) from subfamily I/II, and three sister pairs (SlXTH14/SlXTH6, SlXTH21/SlXTH8, SlXTH26/SlXTH5) from subfamily III (Fig. 1 and 2, see names in blue). The Ka/Ks parameters were evaluated to determine the divergence after duplication. Interestingly, 10 of the 12 sister pairs had Ka/Ks <0.5, which indicates purification selection during evolution. Furthermore, divergence times were estimated to have occurred between 7.4 and 233.33 million years ago (Table S5).

Gene ontology (GO) analysis of SlXTH genes

GO analysis was performed on the entire SlXTH gene family using Blast2GO software (Fig. S5). SlXTH genes are involved in biological processes such as cell wall organization, cell wall biogenesis, and xyloglucan metabolic processes (Fig. S5A). Molecular function and cellular compartment results revealed that all members of the SlXTH family were located in the cell wall and apoplastic region. However, some SlXTH gene members were found to be integral membrane component (Fig. S5B) and had hydrolase and transferase activities (Fig. S5C). The biological processes, molecular functions, and cellular features of each SlXTH protein are specified in Table S6.

Phylogenetic analysis of the SlXTH proteins

One hundred fifty-one full-length XTH protein sequences from S. lycopersicum, S. tuberosum, P. axillaris, N. tabacum and A. thaliana were used to construct a phylogenetic tree based on the neighbor-joining method (Fig. 3; all sequences are provided in File S5). According to this analysis, XTH members are divided into three major subfamilies: an ancestral subfamily (purple branch), subfamily I/II (black branch), and subfamily III, which is divided into subfamilies III-A (pink branch) and III-B (brown branch) (Fig. 3). Nine SlXTHs were clustered in the ancestral subfamily, which includes three SlXTH proteins (SlXTH30, SlXTH31, and SlXTH36). In subfamily III-A, nine XTHs were grouped, two of which were SlXTH proteins. Subfamily III-B contained 21 proteins, six of which were SlXTH proteins. The remaining SlXTHs belonged to subfamily I/II, which includes most of the XTH members from P. axillaris, A. thaliana, S. tuberosum, N. tabacum, and S. lycopersicum (Fig. 3).

Analysis of cis-acting regulatory elements from the SlXTH genes

To further study the potential regulatory elements in the promoters of each member of the SlXTH gene family, 2.0 kb of the promoter sequence of each gene was extracted from the tomato genome database and a cis-acting regulatory element analysis was conducted (File S6). This analysis revealed that SlXTH promoters many regulatory elements, including some involved in cell development, stress-related elements, and hormone regulation (Fig. 4). Methyl jasmonate (MeJa)-responsive regulatory elements were identified in 21 SlXTH promoters (Fig. 4, see brown rectangles with a dot). Defense- and stress-responsive cis-elements were found in SlXTH7, SlXTH15, SlXTH29, SlXTH33, and SlXTH36 (Fig. 4, see purple rectangles with a minus sign). Wound-responsive elements were found in the promoter sequences of SlXTH7, SlXTH8, SlXTH13, SlXTH19, and SlXTH21 genes (Fig. 4, see pink rectangles with an asterisk). All 37 SlXTH genes contain many light-responsive elements (Fig. 4, see a brownish rectangle with a plus sign). A drought-inducible MYB binding site (MBS) was found in almost all SlXTH gene promoters, except SlXTH10, SlXTH11, SlXTH24, SlXTH28, and SlXTH32 (Fig. 4, see the light pink rectangles). The other SlXTH members showed different regulatory elements involved in cell development, with roles in meristem expression, endosperm expression, and the cell cycle (Fig. 4). This result indicates that the SlXTH gene family members are involved in different biological processes and can respond to various biotic and abiotic stresses.

Expression profile of SlXTH genes in selected tomato tissues

The expression patterns of SlXTH genes in different tissues were analyzed using the temporal and spatial expression information from public RNA-seq projects (SGN database) in RPKM values (Fig. 5). Nineteen SlXTH genes were expressed in at least one tissue, while 18 were either not expressed in any of the tested tissues or their expression was relatively low (Fig. 5). SlXTH5, SlXTH7, and SlXTH16 were highly expressed in leaves, whereas SlXTH1, SlXTH2, SlXTH8, and SlXTH11 showed less expression in this tissue, and expression was almost undetected in the other members (Fig. 5). SlXTH16 was highly expressed in leaves, roots, and buds, whereas flowers showed low expression. SlXTH1 and SlXTH21 were mainly expressed in flowers, whereas SlXTH14 was expressed in roots and buds. Some SlXTH members, such as SlXTH6 and SlXTH9, were explicitly expressed in roots. Expression in the other SlXTH genes was either low or undetected (Fig. 5). No RNA-seq study of AM-colonized shoots in tomatoes is available in the SGN database.

Expression profile of SlXTH genes in response to arbuscular mycorrhizal symbiosis

To study the possible role of SlXTH family members in the response of tomato plants to AM colonization, we experimentally evaluated the expression levels of some SlXTH genes in roots and leaves of colonized (M+) and non-colonized (M-) plants by qRT-PCR. Eight SlXTH genes (SlXTH2, SlXTH3, SlXTH6, SlXTH7, SlXTH9, SlXTH14, SlXTH21, and SlXTH35) were selected based on the fact that they presented at least one cis-regulatory element responsive to defense, stress, wounds, and MeJa, which are all known or postulated to be involved in the modulation of plant defense and priming (Fig. S6). In addition, SlXTH17, which does not contain any defense-responsive regulatory elements, was also included in the analysis (Fig. 4 and S6). Microscopy observations of tomato roots revealed that symbiotic structures such as intraradical hyphae, vesicles, and arbuscules were observed in colonized (M+) plants (Fig. 6A, see labels ih, V, and *). As expected, no symbiotic structures were observed in non-colonized (mock, M-) tomato plants (Fig. 6B). In addition, RT-PCR was performed using tomato mycorrhiza-specific phosphate transporters (SlPT4) as a molecular marker of mycorrhiza colonization in tomato roots, and SlPT4 transcript accumulation was only detected in the roots of colonized (M+) plants (Fig. 6C).

Differential expression of SlXTHs was observed in leaves and roots in response to AM symbiosis. In leaves, only SlXTH2 showed higher relative expression in M+ plants as compared to M- plants, whereas SlXTH3, SlXTH6, SlXTH7, SlXTH9, SlXTH14, SlXTH21 and SlXTH35 showed downregulation (Fig. 7). The expression of SlXTH17 was unchanged regardless of the symbiotic status of tomato plants in leaves. In roots, most SlXTH genes exhibit no differential change in expression profile in M+ plants compared to M- plants (Fig. 8). Only SlXTH7 and SlXTH35 were upregulated in response to AM colonization (M+) compared to the control plants (M-), whereas SlXTH3 and SlXTH21 were downregulated. SlXTH17 was not expressed in tomato roots, regardless of the plant’s symbiotic status (Fig. 8). These results indicate that several SlXTH genes likely play critical roles in the tomato response to AM symbiosis, suggesting that these genes, which contain regulatory elements involved in plant defense, could participate in the defense priming process and that they are regulated by AM symbiosis.