The different Malus tissues used for the tissue-specific expression analyses were obtained from 7-year-old “Royal Gala” apple trees at the experimental station of the Shandong Agricultural University in 2020. The roots, stems, and leaves were collected, immediately frozen in liquid nitrogen, and stored at −80°C.

The tissue-cultured “Royal Gala” seedlings used for a low-S and hormone treatment were grown on a Murashige and Skoog (MS) agar medium supplemented with 30 g/L of sucrose and 0.3 mg/L of indole-3-acetic acid (IBA) for 30 days under 16 h light/8 h dark conditions at 26°C. Consistently growing and rooted tissue culture seedlings were selected and transferred into tissue culture glass bottles (250 ml) containing a half-strength Hoagland nutrient solution for treatment. The tissue culture glass bottles were wrapped with black plastic to simulate the dark environment of roots in the soil. The nutrient solution was continuously aerated with an air pump and the dissolved oxygen concentration was maintained at 8–8.5 mg/L. Additionally, the solution was replaced every 2 days. After 7 days of pre-cultivation, the seedlings were prepared for the following treatment: (1) the seedlings were treated in a half-strength Hoagland's nutrient solution (Control, CK) and an half-strength Hoagland's nutrient solution lacking S supplemented with 0.1 mM of MgSO₄ (low S treatment), and the samples of the roots and leaves were collected at 0, 1, 3, 6, 9, 12, and 15 days after the treatment; (2) for the hormone treatments, the seedlings were treated with a half-strength Hoagland's nutrient solution containing either 5 μM of ABA, 2 μM of indole-3-acetic acid (IAA), or 50 μM of methyl jasmonate (MeJA) for 0, 3, 6, 9, and 12 h and 1, 3, and 7 days before the roots were collected. All of the collected samples were quickly frozen in liquid nitrogen and stored at −80°C for a qRT-PCR analysis.

Malus hupehensis seedlings were used for gene cloning. Its seeds were soaked in water for 24 h and then mixed with fine sand at 4°C stratifications for 40 days. The germinated seeds were cultivated in plastic pots filled with soil that was mixed with turf, perlite, and vermiculite (3:1:1) in the greenhouse until they produced eight leaves.

The apple calli (“Orin” cultivar) used for genetic transformation were grown on an MS medium containing 0.5 mg/L of 2,4-dichlorophenoxyacetic acid (2,4-D) and 1.5 mg/L of 6-benzylaminopurine (6-BA) at 26°C in the dark and were sub-cultured in 16-day intervals.

To identify MdSultr genes in the Malus domestica genome, 12 Arabidopsis Sultrs protein sequences were downloaded from The Arabidopsis Information Resource (TAIR) (https://www.arabidopsis.org/) according to the gene accession numbers (Supplementary Table 1) reported previously (Takahashi, 2010) and used as queries to perform local Protein Basic Local Alignment Search Tool (BLASTP) searches in the apple protein database downloaded from Pfam (http://pfam.xfam.org/). Additionally, the Hidden Markov Model (HMM) profiles of the Sulfate_transp (PF00916) and STAS domain (PF01740) were downloaded from the Pfam database and utilized to search for the apple protein database using HMM 3.0 (Harvard University, USA). Afterward, the results were combined and duplicate gene sequences were deleted. Finally, the results were submitted to the Simple Modular Architecture Research Tool (SMART) database (http://smart.embl-heidelberg.de/) and the Conserved Domain Database (CDD) of National Center for Biotechnology Information (NCBI) (https://www.ncbi.nlm.nih.gov/cdd/?term=) to eliminate the genes that did not contain complete Sulfate_transp and STAS domains.

All candidate sequences were submitted to the Expert Protein Analysis System (ExPaSy) Proteomics Server online tools (http://web.expasy.org/protparam/) to predict molecular weight (MW), isoelectric point (PI), the number of amino acids (AAs), and the length of the open reading frame (ORF). The subcellular localization of MdSultrs was predicted using Plant-mPLoc (http://www.csbio.sjtu.edu.cn/bioinf/plant-multi/).

Multiple sequence alignments of Sultr protein sequences in apple and Arabidopsis were generated using the multiple sequence comparison by log-expectation (MUSCLE) program (Edgar, 2004) (GenBank accession numbers see in Supplementary Table 1), and the alignment data were used to construct a neighbor-joining (NJ) phylogenetic tree using MEGA 7.0 (King Abdulaziz University, Jeddah, Saudi Arabia) with the following parameters: 1,000 bootstrap replications, a passion model, and pairwise deletion (Kumar et al., 2016).

The chromosomal location information of the MdSultr genes was obtained from the Genome Database for Rosaceae (GDR) (http://www.rosaceae.org/) and submitted to MapGene2chrom (MG2C) (http://mg2c.iask.in/mg2c_v2.0/) to draw the chromosome location images. The information about the gene structure of MdSultrs was extracted from genome annotation files. The conserved motifs of MdSultrs were analyzed by the MEME program Version 5.3.3 (http://meme-suite.org/tools/meme). The region 2,000 bp upstream of the transcription start site of the MdSultrs was extracted and then submitted to the PlantCARE software (UGent, Gent, Belgium) (http://bioinformatics.psb.ugent.be/webtools/plantcare/html/) to analyze the cis-acting elements related to stress responsiveness and plant hormones on the promoter regions of MdSultrs. Additionally, TBtools (South China Agricultural University, Guangzhou, China) were used for visualization (Chen et al., 2020a).

To investigate the tissue-specific expression of nine MdSultrs, the gene expression profile data of apples were obtained from the National Center for Biotechnology Information (NCBI)-Gene Expression Omnibus (GEO) (https://www.ncbi.nlm.nih.gov/geo/) database (GSE42873). The expression matrix of nine MdSultrs from various tissues in different apple cultivars was extracted and visualized by TBtools (Chen et al., 2020a).

The coding sequence of MhSultr3;1a (GenBank accession No.: MZ634458) was amplified from the cDNA of M. hupehensis roots with specific primers, namely, MhSultr3;1a-F and MhSultr3;1a-R (Supplementary Table 2), using a Phanta Max Super-Fidelity DNA Polymerase Kit (Vazyme, Nanjing, China) according to the instructions of the manufacturer. Then, PCR amplification was performed as follows: 5 min at 95°C, 35 cycles at 94°C for 30 s, 56°C for 30 s, 72°C for 2 min, and a final extension at 72°C for 10 min. The PCR products were checked by sequencing (BGI, Shenzhen, China).

Multiple sequence alignments were performed using DNAman (Lynnon Biosoft, San Ramon, CA, United States). The full-length amino acid sequences of Sultr3;1s from M. domestica, Pyrus bretschneideri, Prunus avium, Prunus persica, Populus trichocarpa, and Arabidopsis (GenBank accession numbers see in Supplementary Table 1) were used for phylogenetic analysis with MEGA7.0 (Kumar et al., 2016).

The coding sequence of MhSultr3;1a was amplified with the specific primers: MhSultr3;1a-EF and MhSultr3;1a-ER (Supplementary Table 2). The amplified fragment was first cloned into the entry vector pDONR222 and subsequently into the binary vector pGWB405-GFP, which are controlled by the CaMV 35S promoter using the Gateway BP and LR recombination reactions, respectively (Invitrogen, Carlsbad, CA, USA). For subcellular localization, Agrobacterium GV3101 containing the 35S:MhSultr3;1a-GFP fusion vector or empty vector 35S: GFP was inserted into the leaf epidermal cells of Nicotiana benthamiana (Chen et al., 2018). The infiltrated tobacco was cultured in the dark for 48 h, and afterward, the infected leaves were cut and observed with a high-resolution laser confocal microscope (ZEISS LSM 880, Jena, Germany).

To construct the plasmid p112AINE-MhSultr3;1a, the coding sequence of MhSultr3;1a was amplified with specific primers, namely, MhSultr3;1aYF and MhSultr3;1aYR (Supplementary Table 2), and then inserted into the EcoRI and BamHI sites of the p112AINE yeast expression vector using the ClonExpress II One Step Cloning Kit (Vazyme, Nanjing, China) (Ding et al., 2016). The recombinant vector MhSultr3;1a-p112AINE and the empty vector p112AINE were transferred into the yeast double SO42- transporter mutant CP154-7A (MATa, his3, leu2, ura3, ade2, trp1, sul1::LEU2, and sul2::URA3) using the lithium acetate method and then selected on a synthetic defined (SD) yeast medium without Tryptophan at 28°C over 3 days (Cherest et al., 1997). For the complementation assays, the transformed yeast was cultured in a liquid SD/-Trp medium until OD₆₀₀ =1 and diluted to four sequential dilutions: 10⁻¹, 10⁻², 10⁻³, and 10⁻⁴; subsequently, 20 μl of each dilution were dropped onto plates with a yeast nitrogen base (YNB) medium (without ammonium SO42-), which was supplemented with essential amino acids and 0.1 mM of SO42- or 0.1 mM of homocysteine as the sole S source. All plates were incubated at 28°C for 3 days to observe growth.

To measure the growth rate of the yeast, the CP154-7A, p112AINE/CP154-7A, and MhSultr3;1a-p112AINE yeasts were incubated in a YNB liquid medium (without ammonium SO42-) that was supplemented with essential amino acids and 0.1 mM of sodium SO42-. The OD₆₀₀ values were measured every 6 h for 84 h using a microplate spectrophotometer (Fisher Scientific, Hampton, NH, USA) with three replicates.

The resulting vector, MhSultr3;1a-GFP, was transformed into Agrobacterium strain LBA4404 using the heat shock method. The transgenic apple calli were generated as described previously (An et al., 2017). Additionally, the transgenic apple calli were confirmed by PCR at DNA and RNA levels (Supplementary Figure 1). A Plant Genomic DNA Kit (Beijing Bioteke Corporation) was used to extract the DNA. The primers used for PCR are listed in Supplementary Table 3.

The 15-day-old wild type (WT) and the transgenic apple calli (OE1 and OE2) were transferred to an MS agar medium (CK) or MS agar medium (without S) containing 0.1 mm of MgSO₄ (low S) and then photographed after 15 days. The fresh weight, SO42- content, and Cys content were measured. The Cys content was measured using a Cys assay kit (Solarbio Science & Technology, Beijing, China). The SO42- content was determined as described previously (Lancilli et al., 2014). The apple calli grown on each plate were used as one biological replicate. A total of three biological replicates were analyzed. For each biological replicate, SO42- content, and Cys content were determined using 0.1 g of calli, respectively.

Total RNA was extracted from “Royal Gala,” M. hupehensis, and apple calli using the polysaccharides- and polyphenolics-rich RNAprep Pure Plant Kit (Tiangen Biotech, Beijing, China). Then, 1 μg of the total RNA was used to synthesize the cDNA using the PrimeScript™ RT reagent Kit and gDNA Eraser (perfect real-time) according to the instructions of the manufacturer (Takara Bio Inc., Shiga, Japan). Afterward, qRT-PCR was performed on a LightCycler^®96 (Roche) using the TB Green^® Premix Ex Taq™ (Tli RNaseH Plus) (TaKaRa, Japan) under the following conditions: 45 cycles of 95°C for 5 s and 60°C for 30 s. All primers used for the qRT-PCR are shown in Supplementary Table 4. Mdactin was used as the reference gene. Three biological and technical replications were assayed and the data were calculated using the 2^−ΔΔCt method.

Statistical analyses were performed with the Student's t-test and one-way ANOVA tests using the DPS 7.05 software (Shenzhen Chinese Technology Co., Ltd). All experiments were carried out in triplicates and expressed as the mean ± SD. The confidence level for statistical significance was p < 0.05. The graphs were made using Microsoft Excel 2010 (Microsoft Corp., Redmond, WA).

To identify MdSultrs members, the protein sequences of 12 AtSultrs, the conserved Sulfate_transp, and STAS domains were used as queries to search the M. domestica genome database using the BLASTP program. Afterward, the results were combined and the duplicate gene sequences were deleted. The genes that did not contain the complete Sulfate_transp and STAS domains were removed using SMART and CDD. Finally, nine putative MdSultrs were identified. As shown in Table 1, the ORF length of the MdSultrs ranged from 1,866 to 4,878 bp, encoding 622–1,625 AAs. MdSultr3;5 had the longest ORF and the most AAs. The predicted MW of the MdSultrs was between 68.2 and 181.1 kDa. The theoretical pI values ranged from 5.52 to 9.04. The subcellular localization prediction results show that MdSultr3;5 was predicted to localize to the nucleus, while MdSultr4;2 was predicted to localize to the chloroplast. The remaining seven MdSultrs were predicted to localize to the cellular membrane. According to the genomic location information obtained from the apple genome database, the nine MdSultrs were mapped onto apple chromosomes (Chrs) 1, 3, 11, 13, 15, 16, and 17. Among them, Chr13 contained the most MdSultr genes (three), whereas the other Chrs contained only one MdSultr gene each (Figure 1A).

To understand the phylogenetic relationship between M. domestica and Arabidopsis, an NJ-phylogenetic tree was constructed based on the protein sequence alignment of the full sequences of the 9 MdSultrs and 12 AtSultrs proteins. As shown in Figures 1B, 2A, the nine putative MdSultrs were divided into two groups, with Group 3 included eight members (MdSultr3;1a, MdSultr3;1b, MdSultr3;1c, MdSultr3;1d, MdSultr3;3a, MdSultr3;3b, MdSultr3;4, and MdSultr3;5), which had 66–76% similarity with Arabidopsis, and the remaining member (MdSultr4;2) was classified as Group 4, which had 66% similarity with Arabidopsis (Table 1).

As shown in Figure 2B, 18 putative conserved motifs with 8~50 amino acid residues were identified using the online MEME tool and subsequently annotated with Pfam and SMART (detailed in Supplementary Table 5). Motifs 5 and 9, which were STAS domains, and Motifs 1, 2, 4, 7, and 11, which were Sulfate_transp domains, were identified in all MdSultr proteins, and the remaining motifs were unknown functional elements. Notably, Motif 18 was only identified in the MdSultr3;3a and MdSultr3;3b proteins. Nine MdSultr protein sequences had highly conserved motifs and conserved motif orders. The similar motif distribution of MdSultrs indicated that they have a close evolutionary relationship and similar functions.

To identify the gene structure differences between MdSultrs, the gene and protein structure of MdSultrs were visualized using TBtools, and the results were consistent with the phylogenetic tree analysis. MdSultr3;5 and MdSultr4;2 had 39 and 19 exons, respectively, while the other MdSultrs had 12–13 each, and all members contained introns (Figure 2C).

To understand the transcriptional regulation mechanisms of the MdSultrs, the 2,000-bp promoter regions of the MdSultrs were isolated for the analysis of the potential cis-elements using the Plant-CARE database. As shown in Figure 2D, all MdSultrs possessed the anaerobic induction response element (ARE). Except for MdSultr3;1c, the other MdSultrs contained the MeJA response element (TGACG-motif/CGTCA-motif). Five MdSultrs had an ABA response element (ABRE) except for MdSultr3;4, MdSultr3;5, and MdSultr4;2. MdSultr3;1a, MdSultr3;1c, MdSultr3;5, and MdSultr4;2 had an auxin response element (TGA-element/AuxRR-core). Moreover, hormone-responsive elements, including gibberellin response (TATC-box/P-box) and salicylic acid (SA) response (TCA-element), and stress elements, such as defense and stress (TCA-rich repeats), low temperature (LTR), and drought (MBS) response elements, were found in the MdSultr promoter regions. These results suggest that MdSultrs may play important roles in the response to hormones and various stressors.

To explore the potential function of MdSultrs in M. domestica, the expression patterns of nine MdSultrs in different tissues, including flowers, fruits, leaves, roots, stems, seeds, and seedlings were analyzed using the GEO database (GSE42873). The results showed that nine MdSultrs was mainly expressed in the leaves, flowers, and fruits of different apple cultivars and hybrids. MdSultr3;1c showed higher expression levels compared with the other eight genes, in the roots of apple cultivars, namely, Golden Delicious (GD) and X8877 (Figure 3A).

In addition, the expression of nine MdSultrs in the roots, stems, and leaves of “Royal Gala” was determined by qRT-PCR. As shown in Figure 3B, MdSultr3;1a, MdSultr3;1d, MdSultr3;4, MdSultr3;5, and MdSultr4;2 had the highest expression levels in the roots;, while MdSultr3;1b, MdSultr3;1c, MdSultr3;3a, and MdSultr3;3b had the highest expression levels in leaves; and MdSultr4;2 was at relatively high levels in the roots and leaves. In summary, the different expression patterns of MdSultrs suggest that MdSultrs may play different roles in the growth and development of Malus.

To estimate the response of MdSultrs to low S, qRT-PCR was used to detect the expression levels of MdSultrs in the roots and leaves within 15 days of a low-S treatment. As shown in Figure 4, the expression patterns of MdSultrs in the roots and leaves were different under a low-S treatment. For example, in the roots, the expression of MdSultr3;1a, MdSultr3;1b, MdSultr3;1c, and MdSultr3;1d were significantly and rapidly increased, while MdSultr3;3a, MdSultr3;3b, MdSultr3;5, and MdSultr4;2 were upregulated after a long-term (12–15 d) low-S treatment (Figure 4A). In leaves, MdSultr3;1a, MdSultr3;1b, MdSultr3;1c, and MdSultr 3;1d were significantly downregulated, while the expression of MdSultr3;3a and MdSultr3;3b significantly increased after 3 days of a low-S treatment. MdSultr3;4, MdSultr3;5, and MdSultr4;2 were rapidly upregulated after 1–3 days of low-S conditions (Figure 4B).

To explore the potential roles of MdSultrs in responding to the hormones in roots, the expression patterns of MdSultrs were investigated under various plant hormones. As shown in Figure 5, the expression of the nine MdSultrs was found to be highly upregulated in response to ABA (Figure 5A). Except for MdSultr3;3a and MdSultr4;2, which were significantly upregulated by IAA, the other MdSultrs were downregulated (Figure 5B). Most MdSultrs were upregulated by the exogenous MeJA except for MdSultr3;4 (Figure 5C). These results indicate that MdSultrs may be involved in the crosstalk of hormone signaling with S metabolism.

Considering that MdSultr3;1a was especially expressed in roots and induced by low S (Figures 3B, 4), and that previous studies have shown that the genes cloned from M. hupehensis are highly homologous to the sequence of the M. domestica genome (Liu et al., 2012; Sun et al., 2017; Zhang et al., 2020), MhSultr3;1a was isolated from the M. hupehensis roots. Multiple alignment showed that the MhSultr3;1a protein was highly homologous to the MdSultr3;1a from apples (99.85% identity), PbSultr3;1 from P. bretschneideri (98.31%), PaSultr3;1 from P. avium (92.60%), PpSultr3;1 from Prunus persica (92.45%), PtSultr3;1 from P. trichocarpa (82.18%), and AtSultr3;1 from Arabidopsis (70.23%) (Supplementary Figure 2A). The phylogenetic tree showed that MhSultr3;1a was homologous to several Sultrs, but was most closely related to MdSultr3;1a (Supplementary Figure 2B). Subcellular localization revealed that MhSultr3;1a was localized to the cell membrane and the nuclear membrane (Figure 6).

To investigate whether MhSultr3;1a had an SO42- transport function, complementation tests were performed in the yeast mutant CP154-7A lacking both SUL1 and SUL2 SO42- transporters. The yeast mutant transformed with MhSultr3;1a and the empty vector p112A1NE grew well with 0.1 mM of homocysteine, but only the yeast-transformed MhSultr3;1a grew well on a medium containing 0.1 mM of sodium SO42- as the sole S source (Figure 7A). In contrast, the CP154-7A-transformed p112A1NE did not grow well, indicating that MhSultr3;1a can rescue the growth defect of CP154-7A. The growth curve was further analyzed to extend the results, as shown in Figure 7B. The CP154-7A-transformed MhSultr3;1a grew well on a liquid YNB medium containing 0.1 mM of SO42-, whereas the CP154-7A- and CP154-7A-transformed p112A1NE failed to grow. Thus, MhSultr3;1a complemented the yeast CP154-7A mutant and encoded a functional SO42- transporter.

To further characterize the function of MhSultr3;1a, two lines of transgenic apple calli (OE1 and OE2) were obtained. As shown in Supplementary Figure 1B, OE1 and OE2 showed relatively higher expression levels compared with the WT. No growth differences were observed between the OE1 and OE2 transgenic and WT apple calli under normal conditions, and their growth was inhibited by low-S treatment. Notably, the OE1 and OE2 transgenic apple calli had larger sizes than the WT when transferred to low S supply conditions (Figure 7C), and the fresh weights of OE1 and OE2 were higher than those of the WT (Figure 7D). The overexpression of MhSultr3;1a improved the growth of yeast and apple calli under low-S conditions.

As shown in Figure 8, the contents of SO42- (Figure 8A) and Cys (Figure 8B) under normal conditions showed no significant differences in the transgenic apple calli overexpressing MhSultr3;1a (OE1 and OE2) and the WT, however, under low-S conditions, the contents of SO42- and Cys in OE1 and OE2 were significantly higher than those in the WT. Especially in OE1, the Cys content was 34.2% higher than that in the WT under low-S conditions. Similarly, the Cys content in OE2 was also higher than that in the WT by 23.06% (Figure 8B). The overexpression of MhSultr3;1a accumulated more SO42- and Cys contents to meet the demands of S-containing compounds under S-limiting conditions.