2023 Volume 3
Article Contents
About this article
ARTICLE   Open Access    

SWEET transporters and their potential roles in response to abiotic and biotic stresses in mulberry

  • # These authors contributed equally: Xiaoru Kang, Shuai Huang, Yuming Feng

More Information
  • Mulberry (Morus spp., Moraceae) is a traditional economic crop plant and is also being gradually utilized as a beverage plant. SWEETs (Sugars Will Eventually be Exported Transporter) are important sugar transporters involved in various biological processes and responses to various stresses. However, SWEETs in mulberry are still poorly studied without a comprehensive functional analysis of SWEETs. In the present study, a total of 24 SWEETs were identified using the Morus alba (Ma) genome. Phylogenetic analysis showed that these 24 MaSWEETs were clustered with SWEETs from Arabidopsis, Populus and Oryza and fell into four clades. MaSWEETs in the same clade are likely to pose similar intron/exon patterns. These MaSWEETs distributed on 12 chromosomes and tandem duplication and block duplication were responsible for the expansion of SWEETs in mulberry. Transmembrane domains and conserved active sites of Tyr and Asp were observed in MaSWEETs. Cis-elements in promoter regions of MaSWEETs indicated the possible function of MaSWEETs in response to hormones and environment stimulus. MaSWEETs showed quite different expression preference in tissues and organs indicating the possible function divergence. In addition, most MaSWEETs showed a disturbed expression levels in response to various abiotic stresses and Ciboria shiraiana infection. MaSWEET1a was functionally characterized as a negative regulator of resistance to C. shiraiana infection based on in vivo transient overexpression of MaSWEET1a in tobacco and down-regulation of MaSWEET1a/b in mulberry. Our results provided foundation for further functional dissection of SWEETs in mulberry and a potential regulator for genetic modification.
  • Columnar cacti are plants of the Cactaceae family distributed across arid and semi-arid regions of America, with ecological, economic, and cultural value[1]. One trait that makes it possible for the columnar cactus to survive in the desert ecosystem is its thick epidermis covered by a hydrophobic cuticle, which limits water loss in dry conditions[1]. The cuticle is the external layer that covers the non-woody aerial organs of land plants. The careful control of cuticle biosynthesis could produce drought stress tolerance in relevant crop plants[2]. In fleshy fruits, the cuticle maintains adequate water content during fruit development on the plant and reduces water loss in fruit during postharvest[3]. Efforts to elucidate the molecular pathway of cuticle biosynthesis have been carried out for fleshy fruits such as tomato (Solanum lycopersicum)[4], apple (Malus domestica)[5], sweet cherry (Prunus avium)[6], mango (Mangifera indica)[7], and pear (Pyrus 'Yuluxiang')[8].

    The plant cuticle is formed by the two main layers cutin and cuticular waxes[3]. Cutin is composed mainly of oxygenated long-chain (LC) fatty acids (FA), which are synthesized by cytochrome p450 (CYP) enzymes. CYP family 86 subfamily A (CYP86A) enzymes carry out the terminal (ω) oxidation of LC-FA[9]. Then, CYP77A carries out the mid-chain oxidation to synthesize the main cutin monomers. In Arabidopsis, AtCYP77A4 and AtCYP77A6 carry out the synthesis of mid-chain epoxy and mid-chain dihydroxy LC-FA, respectively[10,11]. AtCYP77A6 is required for the cutin biosynthesis and the correct formation of floral surfaces[10]. The expression of CYP77A19 (KF410855) and CYP77A20 (KF410856) from potato (Solanum tuberosum) restored the petal cuticular impermeability in Arabidopsis null mutant cyp77a6-1, tentatively by the synthesis of cutin monomers[12]. In eggplant (Solanum torvum), the over-expression of StoCYP77A2 leads to resistance to Verticillium dahlia infection in tobacco plants[13]. Although the function of CYP77A2 in cutin biosynthesis has not yet been tested, gene expression analysis suggests that CaCYP77A2 (A0A1U8GYB0) could play a role in cutin biosynthesis during pepper fruit development[14].

    It has been hypothesized that the export of cuticle precursors is carried out by ATP binding cassette subfamily G (ABCG) transporters. ABCG11/WBC11, ABCG12, and ABCG13 are required for the load of cuticle lipids in Arabidopsis[1517], but ABCG13 function appears to be specific to the flower epidermis[18]. The overexpression of TsABCG11 (JQ389853) from Thellungiella salsugineum increases cuticle amounts and promotes tolerance to different abiotic stresses in Arabidopsis[19].

    Once exported, the cutin monomers are polymerized on the surface of epidermal cells. CD1 code for a Gly-Asp-Ser-Leu motif lipase/esterase (GDSL) from tomato required for the cutin formation through 2-mono(10,16-dihydroxyhexadecanoyl)glycerol esterification[20]. GDSL1 from tomato carries out the ester bond cross-links of cutin monomers located at the cuticle layers and is required for cuticle deposition in tomato fruits[21]. It has been shown that the transcription factor MIXTA-like reduces water loss in tomato fruits through the positive regulation of the expression of CYP77A2, ABCG11, and GDSL1[22]. Despite the relevant role of cuticles in maintaining cactus homeostasis in desert environments[1], the molecular mechanism of cuticle biosynthesis has yet to be described for cactus fruits.

    Stenocereus thurberi is a columnar cactus endemic from the Sonoran desert (Mexico), which produces an ovoid-globose fleshy fruit named sweet pitaya[23]. In its mature state, the pulp of sweet pitaya contains around 86% water with a high content of antioxidants and natural pigments such as betalains and phenolic compounds, which have nutraceutical and industrial relevance[23]. Due to the arid environment in which pitaya fruit grows, studying its molecular mechanism of cuticle biosynthesis can generate new insights into understanding species' adaptation mechanisms to arid environments. Nevertheless, sequences of transcripts from S. thurberi in public databases are scarce.

    RNA-sequencing technology (RNA-seq) allows the massive generation of almost all the transcripts from non-model plants, even if no complete assembled genome is available[24]. Recent advances in bioinformatic tools has improved our capacity to identify long non-coding RNA (lncRNA), which have been showed to play regulatory roles in relevant biological processes, such as the regulation of drought stress tolerance in plants[25], fruit development, and ripening[2629].

    In this study, RNA-seq data were obtained for the de novo assembly and characterization of the S. thurberi fruit peel transcriptome. As a first approach, three transcripts, StCYP77A, StABCG11, and StGDSL1, tentatively involved in cuticle biosynthesis, were identified and quantified during sweet pitaya fruit development. Due to no gene expression analysis having been carried out yet for S. thurberi, stably expressed constitutive genes were identified for the first time.

    Sweet pitaya fruits (S. thurberi) without physical damage were hand harvested from plants in a native conditions field located at Carbó, Sonora, México. They were collocated in a cooler containing dry ice and transported immediately to the laboratory. The superficial part of the peels (~1 mm deep) was removed carefully from the fruits using a scalpel. Peel samples from three fruits were pooled according to their tentative stage of development defined by their visual characteristics, frozen in liquid nitrogen, and pulverized to create a single biological replicate. Four samples belonging to four different plants were analyzed. All fruits harvested were close to the ripening stage. Samples named M1 and M2 were turning from green to ripe [~35−40 Days After Flowering (DAF)], whereas samples M3 and M4 were turning from ripe to overripe (~40−45 DAF).

    Total RNA was isolated from the peels through the Hot Borate method[30]. The concentration and purity of RNA were determined in a spectrophotometer Nanodrop 2000 (Thermo Fisher) by measuring the 260/280 and 260/230 absorbance ratios. RNA integrity was evaluated through electrophoresis in agarose gel 1% and a Bioanalyzer 2100 (Agilent). Pure RNA was sequenced in the paired-end mode in an Illumina NextSeq 500 platform at the University of Arizona Genetics Core Facility. Four RNA-seq libraries, each of them from each sample, were obtained, which include a total of 288,199,704 short reads with a length of 150 base pairs (bp). The resulting sequence data can be accessed at the Sequence Read Archive (SRA) repository of the NCBI through the BioProject ID PRJNA1030439. Libraries are named corresponding to the names of samples M1, M2, M3, and M4.

    FastQC software (www.bioinformatics.babraham.ac.uk/projects/fastqc) was used for short reads quality analysis. Short reads with poor quality were trimmed or eliminated by Trimmomatic (www.usadellab.org/cms/?page=trimmomatic) with a trailing and leading of 25, a sliding window of 4:25, and a minimum read length of 80 bp. A total of 243,194,888 reads with at least a 25 quality score on the Phred scale were used to carry out the de novo assembly by Trinity (https://github.com/trinityrnaseq/trinityrnaseq/wiki) with the following parameters: minimal k-mer coverage of 1, normalization of 50, and minimal transcript length of 200 bp.

    Removal of contaminating sequences and ribosomal RNA (rRNA) was carried out through SeqClean. To remove redundancy, transcripts with equal or more than 90% of identity were merged through CD-hit (www.bioinformatics.org/cd-hit/). Alignment and quantification in terms of transcripts per million (TPM) were carried out through Bowtie (https://bowtie-bio.sourceforge.net/index.shtml) and RSEM (https://github.com/deweylab/RSEM), respectively. Transcripts showing a low expression (TPM < 0.01) were discarded. Assembly quality was evaluated by calculating the parameters N50 value, mean transcript length, TransRate score, and completeness. The statistics of the transcriptome were determined by TrinityStats and TransRate (https://hibberdlab.com/transrate/). The transcriptome completeness was determined through a BLASTn alignment (E value < 1 × 10−3) by BUSCO (https://busco.ezlab.org/) against the database of conserved orthologous genes from Embryophyte.

    To predict the proteins tentatively coded in the S. thurberi transcriptome, the best homology match of the assembled transcripts was found by alignment to the Swiss-Prot, RefSeq, nr-NCBI, PlantTFDB, iTAK, TAIR, and ITAG databases using the BLAST algorithm with an E value threshold of 1 × 10−10 for the nr-NCBI database and of 1 × 10−5 for the others[3134]. An additional alignment was carried out to the protein databases of commercial fruits Persea americana, Prunus persica, Fragaria vesca, Citrus cinensis, and Vitis vinifera to proteins of the cactus Opuntia streptacantha, and the transcriptomes of the cactus Hylocereus polyrhizus, Pachycereus pringlei, and Selenicereus undatus. The list of all databases and the database websites of commercial fruits and cactus are provided in Supplementary Tables S1 & S2. The open reading frame (ORF) of the transcripts and the protein sequences tentative coded from the sweet pitaya transcriptome was predicted by TransDecoder (https://github.com/TransDecoder/TransDecoder/wiki), considering a minimal ORF length of 75 amino acids (aa). The search for protein domains was carried out by the InterPro database (www.ebi.ac.uk/interpro). Functional categorization was carried out by Blast2GO based on GO terms and KEGG metabolic pathways[35].

    LncRNA were identified based on the methods reported in previous studies[25,29,36]. Transcripts without homology to any protein from Swiss-Prot, RefSeq, nr-NCBI, PlantTFDB, iTAK, TAIR, ITAG, P. americana, P. persica, F. vesca, C. cinensis, V. vinifera, and O. streptacantha databases, without a predicted ORF longer than 75 aa, and without protein domains in the InterPro database were selected to identify tentative lncRNA.

    Transcripts coding for signal peptide or transmembrane helices were identified by SignalP (https://services.healthtech.dtu.dk/services/SignalP-6.0/) and TMHMM (https://services.healthtech.dtu.dk/services/TMHMM-2.0/), respectively, and discarded. Further, transcripts corresponding to other non-coding RNAs (ribosomal RNA and transfer RNA) were identified through Infernal by using the Rfam database[37] and discarded. The remaining transcripts were analyzed by CPC[38], and CPC2[39] to calculate their coding potential. Transcripts with a coding potential score lower than −1 for CPC and a coding probability lower than 0.1 for CPC2 were considered lncRNA. To characterize the identified lncRNA, the length and abundance of coding and lncRNA were calculated. Bowtie and RSEM were used to align and quantify raw counts, respectively. The edgeR package[40] was used to normalize raw count data in terms of counts per million (CPM) for both coding and lncRNA.

    To obtain the transcript's expression, the aligning of short reads and quantifying of transcripts were carried out through Bowtie and RSEM software, respectively. A differential expression analysis was carried out between the four libraries by edgeR package in R Studio. Only the transcripts with a count equal to or higher than 0.5 in at least one sample were retained for the analysis. Transcripts with log2 Fold Change (log2FC) between +1 and −1 and with a False Discovery Rate (FDR) lower than 0.05 were taken as not differentially expressed (NDE).

    For the identification of the tentative reference genes two strategies were carried out as described below: i) The NDE transcripts were aligned by BLASTn (E value < 1 × 10−5) to 43 constitutive genes previously reported in fruits from the cactus H. polyrhizus, S. monacanthus, and S. undatus[4143] to identify possible homologous constitutive genes in S. thurberi. Then, the homologous transcripts with the minimal coefficient of variation (CV) were selected; ii) For all the NDE transcripts, the percentile 95 value of the mean CPM and the percentile 5 value of the CV were used as filters to recover the most stably expressed transcripts, based on previous studies[44]. Finally, transcripts to be tested by quantitative reverse transcription polymerase chain reaction (qRT-PCR) were selected based on their homology and tentative biological function.

    The fruit harvesting was carried out as described above. Sweet pitaya fruit takes about 43 d to ripen, therefore, open flowers were tagged, and fruits with 10, 20, 30, 35, and 40 DAF were collected to cover the pitaya fruit development process (Supplementary Fig. S1). The superficial part of the peels (~1 mm deep) was removed carefully from the fruits using a scalpel. Peel samples from three fruits were pooled according to their stage of development defined by their DAF, frozen in liquid nitrogen, and pulverized to create a single biological replicate. One biological replicate consisted of peels from three fruits belonging to the same plant. Two to three biological replicates were evaluated for each developmental stage. Two technical replicates were analyzed for each biological replicate. RNA extraction, quantification, RNA purity, and RNA integrity analysis were carried out as described above.

    cDNA was synthesized from 100 ng of RNA by QuantiTect Reverse Transcription Kit (QIAGEN). Primers were designed using the PrimerQuest™, UNAFold, and OligoAnalyzer™ tools from Integrated DNA Technologies (www.idtdna.com/pages) and following the method proposed by Thornton & Basu[45]. Transcripts quantification was carried out in a QIAquant 96 5 plex according to the PowerUp™ SYBR™ Green Master Mix protocol (Applied Biosystems), with a first denaturation step for 2 min at 95 °C, followed by 40 cycles of denaturation step at 95 °C for 15 s, annealing and extension steps for 30 s at 60 °C.

    The Cycle threshold (Ct) values obtained from the qRT-PCR were analyzed through the algorithms BestKeeper, geNorm, NormFinder, and the delta Ct method[46]. RefFinder (www.ciidirsinaloa.com.mx/RefFinder-master/) was used to integrate the stability results and to find the most stable expressed transcripts in sweet pitaya fruit peel during development. The pairwise variation value (Vn/Vn + 1) was calculated through the geNorm algorithm in R Studio software[47].

    An alignment of 17 reported cuticle biosynthesis genes from model plants were carried out by BLASTx against the predicted proteins from sweet pitaya. Two additional alignments of 17 charaterized cuticle biosynthesis proteins from model plants against the transcripts and predicted proteins of sweet pitaya were carried out by tBLASTn and BLASTp, respectively. An E value threshold of 1 × 10−5 was used, and the unique best hits were recovered for all three alignments. The sequences of the 17 characterized cuticle biosynthesis genes and proteins from model plants are showed in Supplementary Table S3. The specific parameters and the unique best hits for all the alignments carried out are shown in Supplementary Tables S4S8.

    Cuticle biosynthesis-related transcripts tentatively coding for a cytochrome p450 family 77 subfamily A (CYP77A), a Gly-Asp-Ser-Leu motif lipase/esterase 1 (GDSL1), and an ATP binding cassette transporter subfamily G member 11 (ABCG11) were identified by best bi-directional hit according to the functional annotation described above. Protein-conserved domains, signal peptide, and transmembrane helix were predicted through InterProScan, SignalP 6.0, and TMHMM, respectively. Alignment of the protein sequences to tentative orthologous of other plant species was carried out by the MUSCLE algorithm[48]. A neighbor-joining (NJ) phylogenetic tree with a bootstrap of 1,000 replications was constructed by MEGA11[49].

    Fruit sampling, primer design, RNA extraction, cDNA synthesis, and transcript quantification were performed as described above. Relative expression was calculated according to the 2−ΔΔCᴛ method[50]. The sample corresponding to 10 DAF was used as the calibrator. The transcripts StEF1a, StTUA, StUBQ3, and StEF1a + StTUA were used as normalizer genes.

    Normality was assessed according to the Shapiro-Wilk test. Significant differences in the expression of the cuticle biosynthesis-related transcripts between fruit developmental stages were determined by one-way ANOVA based on a completely randomized sampling design and a Tukey honestly significant difference (HSD) test, considering a p-value < 0.05 as significant. Statistical analysis was carried out through the stats package in R Studio.

    RNA was extracted from the peels of ripe sweet pitaya fruits (S. thurberi) from plants located in the Sonoran Desert, Mexico. Four cDNA libraries were sequenced in an Illumina NextSeq 500 platform at the University of Arizona Genetics Core Facility. A total of 288,199,704 reads with 150 base pairs (bp) in length were sequenced in paired-end mode. After trimming, 243,194,888 (84.38%) cleaned short reads with at least 29 mean quality scores per read in the Phred scale and between 80 to 150 bp in length were obtained to carry out the assembly. After removing contaminating sequences, redundancy, and low-expressed transcripts, the assembly included 174,449 transcripts with an N50 value of 2,110 bp. Table 1 shows the different quality variables of the S. thurberi fruit peel transcriptome. BUSCO score showed that 85.4% are completed transcripts, although out of these, 37.2% were found to be duplicated. The resulting sequence data can be accessed at the SRA repository of the NCBI through the BioProject ID PRJNA1030439.

    Table 1.  Quality metrics of the Stenocereus thurberi fruit peel transcriptome.
    Metric Data
    Total transcripts 174,449
    N50 2,110
    Smallest transcript length (bp) 200
    Largest transcript length (bp) 19,114
    Mean transcript length (bp) 1,198.69
    GC (%) 41.33
    Total assembled bases 209,110,524
    TransRate score 0.05
    BUSCO score (%) C: 85.38 (S:48.22, D:37.16),
    F: 10.69, M: 3.93.
    Values were calculated through the TrinityStats function of Trinity and TransRate software. Completeness analysis was carried out through BUSCO by aligning the transcriptome to the Embryophyte database through BLAST with an E value threshold of 1 × 10−3. Complete (C), single (S), duplicated (D), fragmented (F), missing (M).
     | Show Table
    DownLoad: CSV

    A summary of the homology search in the main public protein database for the S. thurberi transcriptome is shown in Supplementary Table S1. From these databases, the higher homologous transcripts were found in RefSeq with 93,993 (53.87 %). Based on the E value distribution, for 41,685 (44%) and 68,853 (49%) of the hits, it was found a strong homology (E value lower than 1 × 10−50) to proteins in the Swiss-Prot and RefSeq databases, respectively (Supplementary Fig. S2a & b). On the other hand, 56,539 (52.34%) and 99,599 (71.11%) of the matches showed a percentage of identity higher than 60% in the Swiss-Prot and RefSeq databases, respectively (Supplementary Fig. S2c & d).

    Figure 1 shows the homology between transcripts from S. thurberi and proteins of commercial fruits, as well as proteins and transcripts of cacti. Transcripts from S. thurberi homologous to proteins from fruits of commercial interest avocado (P. americana), peach (P. persica), strawberry (F. vesca), orange (C. sinensis), and grapefruit (V. vinifera) ranged from 77,285 (44.30%) to 85,421 (48.96%), with 70,802 transcripts homologous to all the five fruit protein databases (Fig. 1a).

    Transcripts homologous to transcripts or proteins from the cactus dragon fruit (H. polyrhizus), prickly pear cactus (O. streptacantha), Mexican giant cardon (P. pringlei), and pitahaya (S. undatus) ranged from 76,238 (43.70%) to 114,933 (65.88%), with 64,009 transcripts homologous to all the four cactus databases (Fig. 1b). Further, out of the total of transcripts, 44,040 transcripts (25.25%) showed homology only to sequences from cactus, but not for model plants Arabidopsis, tomato, or the commercial fruits included in this study (Fig. 1c).

    Figure 1.  Venn diagram of the homology search results against model plants databases, commercial fruits, and cactus. The number in the diagram corresponds to the number of transcripts from S. thurberi homologous to sequences from that plant species. (a) Homologous to sequences from Fragaria vesca (Fa), Persea americana (Pa), Prunus persica (Pp), Vitis vinifera (Vv), and Citrus sinensis (Cs). (b) Homologous to sequences from Opuntia streptacantha (Of), Selenicereus undatus (Su), Hylocereus polyrhizus (Hp), and Pachycereus pringlei (Pap). (c) Homologous to sequences from Solanum lycopersicum (Sl), Arabidopsis thaliana (At), from the commercial fruits (Fa, Pa, Pp, Vv, and Cs), or the cactus included in this study (Of, Su, Hp, and Pap). Homologous searching was carried out by BLAST alignment (E value < 1 × 10−5). The Venn diagrams were drawn by ggVennDiagram in R Studio.

    A total of 45,970 (26.35%), 58,704 (33.65%), and 48,186 (27.65%) transcripts showed homology to transcription factors, transcriptional regulators, and protein kinases in the PlantTFDB, iTAK-TR, and iTAK-PK databases, respectively (Supplementary Tables S1, S9S11). For the PlantTFDB, the homologous transcripts belong to 57 transcriptional factors (TF) families (Fig. 2 & Supplementary Table S9), from which, the most frequent were the basic-helix-loop-helix (bHLH), myeloblastosis-related (MYB-related), NAM, ATAF, and CUC (NAC), ethylene responsive factor (ERF), and the WRKY domain families (WRKY) (Fig. 2).

    Figure 2.  Transcription factor (TF) families distribution of S. thurberi fruit peel transcriptome. The X-axis indicates the number of transcripts with hits to each TF family. Alignment to the PlantTFDB database by BLASTx was carried out with an E value threshold of 1 × 10−5. The bar graph was drawn by ggplot2 in R Studio.

    Based on the homology found and the functional domain searches, gene ontology terms (GO) were assigned to 68,559 transcripts (Supplementary Table S12). Figure 3 shows the top 20 GO terms assigned to the S. thurberi transcriptome, corresponding to the Biological Processes (BP) and Molecular Function (MF) categories. For BP, organic substance metabolic processes, primary metabolic processes, and cellular metabolic processes showed a higher number of transcripts (Supplementary Table S13). Further, for MF, organic cyclic compound binding, heterocyclic compound binding, and ion binding were the processes with the higher number of transcripts. S. thurberi transcripts were classified into 142 metabolic pathways from the KEGG database (Supplementary Table S14). The pathways with the higher number of transcripts recorded were pyruvate metabolism, glycerophospholipid metabolism, glycolysis/gluconeogenesis, and citrate cycle. Further, among the top 20 KEEG pathways, the cutin, suberin, and wax biosynthesis include more than 30 transcripts (Fig. 4).

    Figure 3.  Top 20 Gene Ontology (GO) terms assigned to the S. thurberi fruit peel transcriptome. Bars indicate the number of transcripts assigned to each GO term. Assignment of GO terms was carried out by Blast2GO with default parameters. BP and MF mean Biological Processes and Molecular Functions GO categories, respectively. The graph was drawn by ggplot2 in R Studio.
    Figure 4.  Top 20 KEGG metabolic pathways distribution in the S. thurberi fruit peel transcriptome. Bars indicate the number of transcripts assigned to each KEGG pathway. Assignment of KEGG pathways was carried out in the Blast2GO suite. The bar graph was drawn by ggplot2 in R Studio.

    Out of the total of transcripts, 43,391 (24.87%) were classified as lncRNA (Supplementary Tables S15 & S16). Figure 5 shows a comparison of the length (Fig. 5a) and expression (Fig. 5b) of lncRNA and coding RNA. Both length and expression values were higher in coding RNA than in lncRNA. In general, coding RNA ranged from 201 to 18,629 bp with a mean length of 1,507.18, whereas lncRNA ranged from 200 to 5,198 bp with a mean length of 481.51 (Fig. 5a). The higher expression values recorded from coding RNA and lncRNA were 12.83 and 9.45 log2(CPM), respectively (Fig. 5b).

    Figure 5.  Comparison of coding RNA and long non-coding RNA (lncRNA) from S. thurberi transcriptome. (a) Box plot of transcript length distribution. The Y-axis indicates the length of each transcript in base pairs. (b) Box plot of expression levels. The Y-axis indicates the log2 of the count per million of reads (log2(CPM)) recorded for each transcript. Expression levels were calculated by the edgeR package in R studio. (a), (b) The lines inside the boxes indicate the median. The higher and lower box limits represent the 75th and 25th percentiles, respectively. The box plots were drawn by ggplot2 in R Studio.

    To identify the transcripts without significant changes in expression between the four RNA-seq libraries, a differential expression analysis was carried out. Of the total of transcripts, 4,980 were not differentially expressed (NDE) at least in one paired comparison between the libraries (Supplementary Tables S17S20). Mean counts per million of reads (CPM) and coefficient of variation (CV)[44] were calculated for these NDE transcripts. Transcripts with a CV value lower than 0.113, corresponding with the percentile 5 of the CV, and a mean CPM higher than 1,138.06, corresponding with the percentile 95 of the mean CPM were used as filters to identify the most stably expressed transcripts (Supplementary Table S21). Based on its homology and its tentative biological function, five transcripts were selected to be tested as tentative reference genes. Besides, three NDE transcripts homologous to previously identified stable expressed reference genes in other species of cactus fruit[4143] were selected (Supplementary Table S22). Homology metrics for the eight tentative reference genes selected are shown in Supplementary Table S23. The primer sequences used to amplify the transcripts by qRT-PCR and their nucleotide sequence are shown in Supplementary Tables S24 & S25, respectively.

    The amplification specificity of the eight candidate reference genes determined by melting curves analysis is shown in Supplementary Fig. S3. For the eight tentative reference transcripts selected, the cycle threshold (Ct) values were recorded during sweet pitaya fruit development by qRT-PCR (Supplementary Table S26). The Ct values obtained ranged from 16.85 to 30.26 (Fig. 6a). Plastidic ATP/ADP-transporter (StTLC1) showed the highest Ct values with a mean of 27.34 (Supplementary Table S26). Polyubiquitin 3 (StUBQ3) showed the lowest Ct values in all five sweet pitaya fruit developmental stages (Fig. 6a).

    Figure 6.  Expression stability analysis of tentative reference genes. (a) Box plot of cycle threshold (Ct) distribution of candidate reference genes during sweet pitaya fruit development (10, 20, 30, 35, and 40 d after flowering). The black line inside the box indicates the median. The higher and lower box limits represent the 75th and 25th percentiles, respectively. (b) Bar chart of the geometric mean (geomean) of ranking values calculated by RefFinder for each tentative reference gene (X-axis). The lowest values indicate the best reference genes. (c) Bar chart of the pairwise variation analysis and determination of the optimal number of reference genes by the geNorm algorithm. A pairwise variation value lower than 0.15 indicates that the use of Vn/Vn + 1 reference genes is reliable for the accurate normalization of qRT-PCR data. The Ct data used in the analysis were calculated by qRT-PCR in a QIAquant 96 5 plex (QIAGEN) according to the manufacturer's protocol. The box plot and the bar graphs were drawn by ggplot2 and Excel programs, respectively. Abbreviations: Actin 7 (StACT7), alpha-tubulin (StTUA), elongation factor 1-alpha (StEF1a), COP1-interactive protein 1 (StCIP1), plasma membrane ATPase 4 (StPMA4), BEL1-like homeodomain protein 1 (StBLH1), polyubiquitin 3 (StUBQ3), and plastidic ATP/ADP-transporter (StTLC1).

    The best stability values calculated by NormFinder were 0.45, 0.51, 0.97, and 0.99, corresponding to the transcripts elongation factor 1-alpha (StEF1a), alpha-tubulin (StTUA), plastidic ATP/ADP-transporter (StTLC1), and actin 7 (StACT7), respectively (Supplementary Table S27). For BestKeeper, the most stable expressed transcripts were StUBQ3, StTUA, and StEF1a, with values of 0.72, 0.75, and 0.87, respectively. In the case of the delta Ct method[51], the transcripts StEF1a, StTUA, and StTLC1 showed the best stability.

    According to geNorm analysis, the most stable expressed transcripts were StTUA, StEF1a, StUBQ3, and StACT7, with values of 0.74, 0.74, 0.82, and 0.96, respectively. All the pairwise variation values (Vn/Vn + 1) were lower than 0.15, ranging from 0.019 for V2/V3 to 0.01 for V6/V7 (Fig. 6c). The V value of 0.019 obtained for V2/V3 indicates that the use of the best two reference genes StTUA and StEF1a is reliable enough for the accurate normalization of qRT-PCR data, therefore no third reference gene is required[47]. Except for BestKeeper analysis, StEF1a and StTUA were the most stable transcripts for all of the methods carried out in this study (Supplementary Table S27). The comprehensive ranking analysis indicates that StEF1a, followed by StTUA and StUBQ3, are the most stable expressed genes and are stable enough to be used as reference genes in qRT-PCR analysis during sweet pitaya fruit development (Fig. 6b).

    Three cuticle biosynthesis-related transcripts TRINITY_DN17030_c0_g1_i2, TRINITY_DN15394_c0_g1_i1, and TRINITY_DN23528_c1_g1_i1 tentatively coding for the enzymes cytochrome p450 family 77 subfamily A (CYP77A), Gly-Asp-Ser-Leu motif lipase/esterase 1 (GDSL1), and an ATP binding cassette transporter subfamily G member 11 (ABCG11/WBC11), respectively, were identified and quantified. The nucleotide sequence and predicted amino acid sequences of the three transcripts are shown in Supplementary File 1. The best homology match for StCYP77A (TRINITY_DN17030_c0_g1_i2) was for AtCYP77A4 (AT5G04660) from Arabidopsis and SmCYP77A2 (P37124) from eggplant (Solanum melongena) in the TAIR and Swiss-Prot databases, respectively (Supplementary Table S23).

    TransDecoder, InterPro, and TMHMM analysis showed that StCYP77A codes a polypeptide of 518 amino acids (aa) in length that comprises a cytochrome P450 E-class domain (IPR002401) and a transmembrane region (residues 10 to 32). The phylogenetic tree constructed showed that StCYP77A is grouped in a cluster with all the CYP77A2 proteins included in this analysis, being closer to CYP77A2 (XP_010694692) from B. vulgaris and Cgig2_012892 (KAJ8441854) from Carnegiea gigantean (Supplementary Fig. S4).

    StGDSL1 (TRINITY_DN15394_c0_g1_i1) alignment showed that it is homologous to a GDSL esterase/lipase from Arabidopsis (Q9LU14) and tomato (Solyc03g121180) (Supplementary Table S23). TransDecoder, InterPro, and SignalP analysis showed that StGDSL1 codes a polypeptide of 354 aa in length that comprises a GDSL lipase/esterase domain IPR001087 and a signal peptide with a cleavage site between position 25 and 26 (Supplementary Fig. S5).

    Supplementary Figure S6 shows the analysis carried out on the predicted amino acid sequence of StABCG11 (TRINITY_DN23528_c1_g1_i1). The phylogenetic tree constructed shows three clades corresponding to the ABCG13, ABCG12, and ABCG11 protein classes with bootstrap support ranging from 40% to 100% (Supplementary Fig. S6a). StABCG11 is grouped with all the ABCG11 transporters included in this study in a well-separated clade, being closely related to its tentative ortholog from C. gigantean Cgig2_004465 (KAJ8441854). InterPro and TMHMM results showed that the StABCG11 sequence contains an ABC-2 type transporter transmembrane domain (IPR013525; PF01061.27) with six transmembrane helices (Supplementary Fig. S6b).

    The predicted protein sequence of StABCG11 is 710 aa in length, holding the ATP binding domain (IPR003439; PF00005.30) and the P-loop containing nucleoside triphosphate hydrolase domain (IPR043926; PF19055.3) of the ABC transporters of the G family. Multiple sequence alignment shows that the Walker A and B motif sequence and the ABC signature[15] are conserved between the ABCG11 transporters from Arabidopsis, tomato, S. thurberi, and C. gigantean (Supplementary Fig. S6c).

    According to the results of the expression stability analysis (Fig. 6), four normalization strategies were tested to quantify the three cuticle biosynthesis-related transcripts during sweet pitaya fruit development. The four strategies consist of normalizing by StEF1a, StTUA, StUBQ3, or StEF1a+StTUA. Primer sequences used to quantify the transcripts StCYP77A (TRINITY_DN17030_c0_g1_i2), StGDSL1 (TRINITY_DN15394_c0_g1_i1), and StABCG11 (TRINITY_DN23528_c1_g1_i1) by qRT-PCR during sweet pitaya fruit development are shown in Supplementary Table S24.

    The three cuticle biosynthesis-related transcripts showed differences in expression during sweet pitaya fruit development (Supplementary Table S28). The same expression pattern was recorded for the three cuticle biosynthesis transcripts when normalization was carried out by StEF1a, StTUA, StUBQ3, or StEF1a + StTUA (Fig. 7). A higher expression of StCYP77A and StGDSL1 are shown at the 10 and 20 DAF, showing a decrease at 30, 35, and 40 DAF. StABCG11 showed a similar behavior, with a higher expression at 10 and 20 DAF and a reduction at 30 and 35 DAF. Nevertheless, unlike StCYP77A and StGDSL1, a significant increase at 40 DAF, reaching the same expression as compared with 10 DAF, is shown for StABCG11 (Fig. 7).

    Figure 7.  Expression analysis of cuticle biosynthesis-related transcripts StCYP77A, StGDSL1, and StABCG11 during sweet pitaya (Stenocereus thurberi) fruit development. Relative expression was calculated through the 2−ΔΔCᴛ method using elongation factor 1-alpha (StEF1a), alpha-tubulin (StTUA), polyubiquitin 3 (StUBQ3), or StEF1a + StTUA as normalizing genes at 10, 20, 30, 35, and 40 d after flowering (DAF). The Y-axis and error bars represent the mean of the relative expression ± standard error (n = 4−6) for each developmental stage in DAF. The Ct data for the analysis was recorded by qRT-PCR in a QIAquant 96 5 plex (QIAGEN) according to the manufacturer's protocol. The graph line was drawn by ggplot2 in R Studio. Abbreviations: cytochrome p450 family 77 subfamily A (StCYP77A), Gly-Asp-Ser-Leu motif lipase/esterase 1 (StGDSL1), and ATP binding cassette transporter subfamily G member 11 (StABCG11).

    Characteristics of a well-assembled transcriptome include an N50 value closer to 2,000 bp, a high percentage of conserved transcripts completely assembled (> 80%), and a high proportion of reads mapping back to the assembled transcripts[52]. In the present study, the first collection of 174,449 transcripts from S. thurberi fruit peel are reported. The generated transcriptome showed an N50 value of 2,110 bp, a TransRate score of 0.05, and a GC percentage of 41.33 (Table 1), similar to that reported for other de novo plant transcriptome assemblies[53]. According to BUSCO, 85.4% of the orthologous genes from the Embryophyta databases completely matched the S. thurberi transcriptome, and only 3.9% were missing (Table 1). These results show that the S. thurberi transcriptome generated is not fragmented, and it is helpful in predicting the sequence of almost all the transcripts expressed in sweet pitaya fruit peel[24].

    The percentage of transcripts homologous found, E values, and identity distribution (Supplementary Tables S1 & S2; Supplementary Fig. S2) were similar to that reported in the de novo transcriptome assembly for non-model plants and other cactus fruits[4143,54] and further suggests that the transcriptome assembled of S. thurberi peel is robust[52]. Of the total of transcripts, 70,802 were common to all the five commercial fruit protein databases included in this study, which is helpful for the search for conserved orthologous involved in fruit development and ripening (Fig. 2a). A total of 34,513 transcripts (20%) show homology only to sequences in the cactus's databases, but not in the others (Supplementary Tables S1 & S2; Fig. 1c). This could suggest that a significant conservation of sequences among plants of the Cactaceae family exists which most likely are to have a function in this species adaptation to desert ecosystems.

    To infer the biological functionality represented by the S. thurberi fruit peel transcriptome, gene ontology (GO) terms and KEGG pathways were assigned. Of the main metabolic pathways assigned, 'glycerolipid metabolism' and 'cutin, suberine, and wax biosynthesis' suggests an active cuticle biosynthesis in pitaya fruit peel (Fig. 4). In agreement with the above, the main GO terms assigned for the molecular function (MF) category were 'organic cyclic compound binding', 'transmembrane transporter activity', and 'lipid binding' (Fig. 3). For the biological processes (BP) category, the critical GO terms for the present research are 'cellular response to stimulus', 'response to stress', 'anatomical structure development', and 'transmembrane transport', which could suggest the active development of the fruit epidermis and cuticle biosynthesis for protection to stress.

    The most frequent transcription factors (TF) families found in S. thurberi transcriptome were NAC, WRKY, bHLH, ERF, and MYB-related (Fig. 2), which had been reported to play a function in the tolerance to abiotic stress in plants[55,56]. Although the role of NAC, WRKY, bHLH, ERF, and MYB TF in improving drought tolerance in relevant crop plants has been widely documented[57,58], their contribution to the adaptation of cactus to arid ecosystems has not yet been elucidated and further experimental pieces of evidence are needed.

    It has been reported that the heterologous expression of ERF TF from Medicago truncatula induces drought tolerance and cuticle wax biosynthesis in Arabidopsis leaf[59]. In tomato fruits, the gene SlMIXTA-like which encodes a MYB transcription factor avoids water loss through the positive regulation of genes related to the biosynthesis and transport of cuticle compounds[22]. Despite the relevant role of cuticles in maintaining cactus physiology in desert environments, experimental evidence showing the role of the different TF-inducing cuticle biosynthesis has yet to be reported for cactus fruits.

    Out of the transcripts, 43,391 were classified as lncRNA (Supplementary Tables S15 & S16). This is the first report of lncRNA identification for the species S. thurberi. In fruits, 3,679 lncRNA has been identified from tomato[26], 3,330 from peach (P. persica)[29], 3,857 from melon (Cucumis melo)[28], 2,505 from hot pepper (Capsicum annuum)[27], and 3,194 from pomegranate (Punica granatum)[36]. Despite the stringent criteria to classify the lncRNA of sweet pitaya fruit (S. thurberi), a higher number of lncRNAs are shown when compared with previous reports. This finding is most likely due to the higher level of redundancy found during the transcriptome analysis. To reduce this redundancy, further efforts to achieve the complete genome assembly of S. thurberi are needed.

    Previous studies showed that lncRNA is shorter and has lower expression levels than coding RNA[6062]. In agreement with those findings, both the length and expression values of lncRNA from S. thurberi were lower than coding RNA (Fig. 5). It has been suggested that lncRNA could be involved in the biosynthesis of cuticle components in cabbage[61] and pomegranate[36] and that they could be involved in the tolerance to water deficit through the regulation of cuticle biosynthesis in wild banana[60]. Nevertheless, the molecular mechanism by which lncRNA may regulate the cuticle biosynthesis in S. thurberi fruits has not yet been elucidated.

    A relatively constant level of expression characterizes housekeeping genes because they are involved in essential cellular functions. These genes are not induced under specific conditions such as biotic or abiotic stress. Because of this, they are very useful as internal reference genes for qRT-PCR data normalization[63]. Nevertheless, their expression could change depending on plant species, developmental stages, and experimental conditions[64]. Reliable reference genes for a specific experiment in a given species must be identified to carry out an accurate qRT-PCR data normalization[63]. An initial screening of the transcript expression pattern through RNA-seq improves the identification of stably expressed transcripts by qRT-PCR[44,64].

    Identification of stable expressed reference transcripts during fruit development has been carried out in blueberry (Vaccinium bracteatum)[65], kiwifruit (Actinidia chinensis)[66], peach (P. persica)[67], apple (Malus domestica)[68], and soursop (Annona muricata)[69]. These studies include the expression stability analysis through geNorm, NormFinder, and BestKeeper algorithms[68,69], some of which are supported in RNA-seq data[65,66]. Improvement of expression stability analysis by RNA-seq had led to the identification of non-previously reported reference genes with a more stable expression during fruit development than commonly known housekeeping genes in grapevine (V. vinifera)[44], pear (Pyrus pyrifolia and P. calleryana)[64], and pepper (C. annuum)[70].

    For fruits of the Cactaceae family, only a few studies identifying reliable reference genes have been reported[4143]. Mainly because gene expression analysis has not been carried out previously for sweet pitaya (S. thurberi), the RNA-seq data generated in this work along with geNorm, NormFinder, BestKeeper, and RefFinder algorithms were used to identify reliable reference genes. The comprehensive ranking analysis showed that out of the eight candidate genes tested, StEF1a followed by StTUA and StUBQ3 were the most stable (Fig. 6b). All the pairwise variation values (Vn/Vn + 1) were lower than 0.15 (Fig. 6c), which indicates that StEF1a, StTUA, and StUBQ3 alone or the use of StEF1a and StTUA together are reliable enough to normalize the gene expression data generated by qRT-PCR.

    The genes StEF1a, StTUA, and StUBQ3 are homologous to transcripts found in the cactus species known as dragonfruit (Hylocereus monacanthus and H. undatus)[41], which have been tested as tentative reference genes during fruit development. EF1a has been proposed as a reliable reference gene in the analysis of changes in gene expression of dragon fruit (H. monacanthus and H. undatus)[41], peach (P. persica)[67], apple (M. domestica)[68], and soursop (A. muricata)[69]. According to the expression stability analysis carried out in the present study (Fig. 6) four normalization strategies were designed. The same gene expression pattern was recorded for the three target transcripts evaluated when normalization was carried out by the genes StEF1a, StTUA, StUBQ3, or StEF1a + StTUA (Fig. 7). Further, these data indicates that these reference genes are reliable enough to be used in qRT-PCR experiments during fruit development of S. thurberi.

    The plant cuticle is formed by two main layers: the cutin, composed mainly of mid-chain oxygenated LC fatty acids, and the cuticular wax, composed mainly of very long-chain (VLC) fatty acids, and their derivates VLC alkanes, VLC primary alcohols, VLC ketones, VLC aldehydes, and VLC esters[3]. In Arabidopsis CYP77A4 and CYP77A6 catalyze the synthesis of midchain epoxy and hydroxy ω-OH long-chain fatty acids, respectively[10,11], which are the main components of fleshy fruit cuticle[3].

    The functional domain search carried out in the present study showed that StCYP77A comprises a cytochrome P450 E-class domain (IPR002401) and a membrane-spanning region from residues 10 to 32 (Supplementary Fig. S4). This membrane-spanning region has been previously characterized in CYP77A enzymes from A. thaliana and Brassica napus[11,71]. It suggests that the protein coded by StCYP77A could catalyze the oxidation of fatty acids embedded in the endoplasmic reticulum membrane of the epidermal cells of S. thurberi fruit. Phylogenetic analysis showed that StCYP77A was closer to proteins from its phylogenetic-related species B. vulgaris (BvCYP772; XP_010694692) and C. gigantea (Cgig2_012892) (Supplementary Fig. S4). StCYP77A, BvCYP77A2, and Cgig2_012892 were closer to SlCYP77A2 and SmCYP77A2 than to CYP77A4 and CYP77A6 proteins, suggesting that StCYP77A (TRINITY_DN17030_c0_g1_i2) could correspond to a CYP77A2 protein.

    Five CYP77A are present in the Arabidopsis genome, named CYP77A4, CYP77A5, CYP77A6, CYP77A7, and CYP77A9, but their role in cuticle biosynthesis has only been reported for CYP77A4 and CYP77A6[72]. It has been suggested that CYP77A2 from eggplant (S. torvum) could contribute to the defense against fungal phytopathogen infection by the synthesis of specific compounds[13]. In pepper fruit (C. annuum), the expression pattern of CYP77A2 (A0A1U8GYB0) and ABCG11 (LOC107862760) suggests a role of CYP77A2 and ABCG11 in cutin biosynthesis at the early stages of pepper fruit development[14].

    In the case of the protein encoded by StGDSL1 (354 aa), the length found in this work is similar to the length of its homologous from Arabidopsis (AT3G16370) and tomato (Solyc03g121180) (Supplementary Fig. S5). A GDSL1 protein named CD1 polymerizes midchain oxygenated ω-OH long-chain fatty acids to form the cutin polyester in the extracellular space of tomato fruit peel[20,21]. It has been suggested that the 25-amino acid N-signal peptide found in StGDSL1 (Supplementary Fig. S5), previously reported in GDSL1 from Arabidopsis, B. napus, and tomato, plays a role during the protein exportation to the extracellular space[21,73].

    A higher expression of StCYP77A, StGDSL1, and StABCG11 is shown at the 10 and 20 DAF of sweet pitaya fruit development (Fig. 7), suggesting the active cuticle biosynthesis at the early stages of sweet pitaya fruit development. In agreement with that, two genes coding for GDSL lipase/hydrolases from tomato named SGN-U583101 and SGN-U579520 are highly expressed in the early stages and during the expansion stages of tomato fruit development, but their expression decreases in later stages[74]. It has been shown that the expression of GDSL genes, like CD1 from tomato, is higher in growing fruit[20,21]. Like tomato, the increase in expression of StCYP77A and StGDSL1 shown in pitaya fruit development could be due to an increase in cuticle deposition during the expansion of the fruit epidermis[20].

    The phylogenetic analysis, the functional domains, and the six transmembrane helices found in the StABCG11 predicted protein (Supplementary Fig. S6), suggests that it is an ABCG plasma membrane transporter of sweet pitaya[15]. Indeed, an increased expression of StABCG11 at 40 DAF was recorded in the present study (Fig. 7). Further, this data strongly suggests that it could be playing a relevant role in the transport of cuticle components at the beginning and during sweet pitaya fruit ripening.

    In Arabidopsis, ABCG11 (WBC11) exports cuticular wax and cutin compounds from the plasma membrane[15,75]. It has been reported that a high expression of the ABC plasma membrane transporter from mango MiWBC11 correlates with a higher cuticle deposition during fruit development[7]. The expression pattern for StABCG11, StCYP77A, and StGDSL1 suggests a role of StABCG11 as a cutin compound transporter in the earlier stages of sweet pitaya fruit development (Fig. 7). Further, its increase at 40 DAF suggests that it could be transporting cuticle compounds other than oxygenated long-chain fatty acids, or long-chain fatty acids that are not synthesized by StCYP77A and StGDSL1 in the later stages of fruit development.

    Like sweet pitaya, during sweet cherry fruit (Prunus avium) development, the expression of PaWCB11, homologous to AtABCG11 (AT1G17840), increases at the earlier stages of fruit development decreases at the intermediate stages, and increases again at the later stages[76]. PaWCB11 expression correlated with cuticle membrane deposition at the earlier and intermediate stages of sweet cherry fruit development but not at the later[76]. The increased expression of StABCG11 found in the present study could be due to the increased transport of cuticular wax compounds, such as VLC fatty acids and their derivates, in the later stages of sweet pitaya development[15,75].

    Cuticular waxes make up the smallest amount of the fruit cuticle. Even so, they mainly contribute to the impermeability of the fruit's epidermis[3]. An increase in the transport of cuticular waxes at the beginning of the ripening stage carried out by ABCG transporters could be due to a greater need to avoid water loss and to maintain an adequate amount of water during the ripening of the sweet pitaya fruit. Nevertheless, further expression analysis of cuticular wax biosynthesis-related genes, complemented with chemical composition analysis of cuticles could contribute to elucidating the molecular mechanism of cuticle biosynthesis in cacti and their physiological contribution during fruit development.

    In this study, the transcriptome of the sweet pitaya (S. thurberi) fruit peel was assembled for the first time. The reference genes found here are a helpful tool for further gene expression analysis in sweet pitaya fruit. Transcripts tentatively involved in cuticle compound biosynthesis and transport are reported for the first time in sweet pitaya. The results suggest a relevant role of cuticle compound biosynthesis and transport at the early and later stages of fruit development. The information generated will help to improve the elucidation of the molecular mechanism of cuticle biosynthesis in S. thurberi and other cactus species in the future. Understanding the cuticle's physiological function in the adaptation of the Cactaceae family to harsh environmental conditions could help design strategies to increase the resistance of other species to face the increase in water scarcity for agricultural production predicted for the following years.

    The authors confirm contribution to the paper as follows: study conception and design: Tiznado-Hernández ME, Tafolla-Arellano JC, García-Coronado H, Hernández-Oñate MÁ; data collection: Tiznado-Hernández ME, Tafolla-Arellano JC, García-Coronado H, Hernández-Oñate MÁ; analysis and interpretation of results: Tiznado-Hernández ME, García-Coronado H, Hernández-Oñate MÁ, Burgara-Estrella AJ; draft manuscript preparation: Tiznado-Hernández ME, García-Coronado H. All authors reviewed the results and approved the final version of the manuscript.

    All data generated or analyzed during this study are included in this published article and its supplementary information files. The sequence data can be accessed at the Sequence Read Archive (SRA) repository of the NCBI through the BioProject ID PRJNA1030439.

    The authors wish to acknowledge the financial support of Consejo Nacional de Humanidades, Ciencias y Tecnologías de México (CONAHCYT) through project number 579: Elucidación del Mecanismo Molecular de Biosíntesis de Cutícula Utilizando como Modelo Frutas Tropicales. We appreciate the University of Arizona Genetics Core and Illumina for providing reagents and equipment for library sequencing. The author, Heriberto García-Coronado (CVU 490952), thanks the CONAHCYT (acronym in Spanish) for the Ph.D. scholarship assigned (749341). The author, Heriberto García-Coronado, thanks Dr. Edmundo Domínguez-Rosas for the technical support in bioinformatics for identifying long non-coding RNA.

  • The authors declare that they have no conflict of interest.

  • Supplemental Table S1 Primers used in this study.
    Supplemental Table S2 Putative cis-elements of MaSWEETs promoters.
    Supplemental Fig. S1 Prediction of MaSWEETs transmbrane region.
  • [1]

    Liu Y, Song Y, Ruan Y. 2022. Sugar conundrum in plant - pathogen interactions: roles of invertase and sugar transporters depend on pathosystems. Journal of Experimental Botany 73:1910−25

    doi: 10.1093/jxb/erab562

    CrossRef   Google Scholar

    [2]

    Chen LQ, Hou BH, Lalonde S, Takanaga H, Hartung ML, et al. 2010. Sugar transporters for intercellular exchange and nutrition of pathogens. Nature 468:527−32

    doi: 10.1038/nature09606

    CrossRef   Google Scholar

    [3]

    Chen L, Qu X, Hou BH, Sosso D, Osorio S, et al. 2012. Sucrose efflux mediated by SWEET proteins as a key step for phloem transport. Science 335:207−11

    doi: 10.1126/science.1213351

    CrossRef   Google Scholar

    [4]

    Eom JS, Chen L, Sosso D, Julius BT, Lin IW, et al. 2015. SWEETs, transporters for intracellular and intercellular sugar translocation. Current Opinion in Plant Biology 25:53−62

    doi: 10.1016/j.pbi.2015.04.005

    CrossRef   Google Scholar

    [5]

    Misra VA, Wafula EK, Wang Y, DePamphilis CW, Timko MP. 2019. Genome-wide identification of MST, SUT and SWEET family sugar transporters in root parasitic angiosperms and analysis of their expression during host parasitism. BMC Plant Biology 19:196

    doi: 10.1186/s12870-019-1786-y

    CrossRef   Google Scholar

    [6]

    Lin IW, Sosso D, Chen L, Gase K, Kim SG, et al. 2014. Nectar secretion requires sucrose phosphate synthases and the sugar transporter SWEET9. Nature 508:546−49

    doi: 10.1038/nature13082

    CrossRef   Google Scholar

    [7]

    Feng C, Han J, Han X, Jiang J. 2015. Genome-wide identification, phylogeny, and expression analysis of the SWEET gene family in tomato. Gene 573:261−72

    doi: 10.1016/j.gene.2015.07.055

    CrossRef   Google Scholar

    [8]

    Gao Y, Wang Z, Kumar V, Xu X, Yuan D, et al. 2018. Genome-wide identification of the SWEET gene family in wheat. Gene 642:284−92

    doi: 10.1016/j.gene.2017.11.044

    CrossRef   Google Scholar

    [9]

    Hu B, Wu H, Huang W, Song J, Zhou Y, et al. 2019. SWEET gene family in Medicago truncatula: genome-wide identification, expression and substrate specificity analysis. Plants 8:338

    doi: 10.3390/plants8090338

    CrossRef   Google Scholar

    [10]

    Wei Y, Xiao D, Zhang C, Hou X. 2019. The expanded SWEET gene family following whole genome triplication in Brassica rapa. Genes 10:722

    doi: 10.3390/genes10090722

    CrossRef   Google Scholar

    [11]

    Huang D, Chen Y, Liu X, Ni D, Bai L, et al. 2022. Genome-wide identification and expression analysis of the SWEET gene family in daylily (Hemerocallis fulva) and functional analysis of HfSWEET17 in response to cold stress. BMC Plant Biology 22:211

    doi: 10.1186/s12870-022-03609-6

    CrossRef   Google Scholar

    [12]

    Zhang L, Wang L, Zhang J, Song C, Li Y, et al. 2020. Expression and localization of SWEETs in Populus and the effect of SWEET7 overexpression in secondary growth. Tree Physiology 41:882−99

    doi: 10.1093/treephys/tpaa145

    CrossRef   Google Scholar

    [13]

    Jiang L, Song C, Zhu X, Yang J. 2021. SWEET transporters and the potential functions of these sequences in tea (Camellia sinensis). Frontiers in Genetics 12:655843

    doi: 10.3389/fgene.2021.655843

    CrossRef   Google Scholar

    [14]

    Xu Y, Tao Y, Cheung LS, Fan C, Chen L, et al. 2014. Structures of bacterial homologues of SWEET transporters in two distinct conformations. Nature 515:448−52

    doi: 10.1038/nature13670

    CrossRef   Google Scholar

    [15]

    Breia R, Conde A, Badim H, Fortes AM, Gerós H, et al. 2021. Plant SWEETs: from sugar transport to plant - pathogen interaction and more unexpected physiological roles. Plant Physiology 186:836−52

    doi: 10.1093/plphys/kiab127

    CrossRef   Google Scholar

    [16]

    Sosso D, Luo D, Li Q, Sasse J, Yang J, et al. 2015. Seed filling in domesticated maize and rice depends on SWEET-mediated hexose transport. Nature Genetics 47:1489−93

    doi: 10.1038/ng.3422

    CrossRef   Google Scholar

    [17]

    Miao L, Lv Y, Kong L, Chen Q, Chen C, et al. 2018. Genome-wide identification, phylogeny, evolution, and expression patterns of MtN3/saliva/SWEET genes and functional analysis of BcNS in Brassica rapa. BMC Genomics 19:174

    doi: 10.1186/s12864-018-4554-8

    CrossRef   Google Scholar

    [18]

    Yang J, Luo D, Yang B, Frommer WB, Eom JS. 2018. SWEET11 and 15 as key players in seed filling in rice. New Phytologist 218:604−15

    doi: 10.1111/nph.15004

    CrossRef   Google Scholar

    [19]

    Klemens PAW, Patzke K, Deitmer J, Spinner L, le Hir R, et al. 2013. Overexpression of the vacuolar sugar carrier AtSWEET16 modifies germination, growth, and stress tolerance in Arabidopsis. Plant Physiology 163:1338−52

    doi: 10.1104/pp.113.224972

    CrossRef   Google Scholar

    [20]

    Le Hir R, Spinner L, Klemens PA, Chakraborti D, de Marco F, et al. 2015. Disruption of the sugar transporters AtSWEET11 and AtSWEET12 affects vascular development and freezing tolerance in Arabidopsis. Molecular Plant 8:1687−90

    doi: 10.1016/j.molp.2015.08.007

    CrossRef   Google Scholar

    [21]

    Seo PJ, Park JM, Kang SK, Kim SG, Park CM. 2011. An Arabidopsis senescence-associated protein SAG29 regulates cell viability under high salinity. Planta 233:189−200

    doi: 10.1007/s00425-010-1293-8

    CrossRef   Google Scholar

    [22]

    Wang L, Yao L, Hao X, Li N, Qian W, et al. 2018. Tea plant SWEET transporters: expression profiling, sugar transport, and the involvement of CsSWEET16 in modifying cold tolerance in Arabidopsis. Plant Molecular Biology 96:577−92

    doi: 10.1007/s11103-018-0716-y

    CrossRef   Google Scholar

    [23]

    Yao L, Ding C, Hao X, Zeng J, Yang Y, et al. 2020. CsSWEET1a and CsSWEET17 mediate growth and freezing tolerance by promoting sugar transport across the plasma membrane. Plant and Cell Physiology 61:1669−82

    doi: 10.1093/pcp/pcaa091

    CrossRef   Google Scholar

    [24]

    Cox KL, Meng F, Wilkins KE, Li F, Wang P, et al. 2017. TAL effector driven induction of a SWEET gene confers susceptibility to bacterial blight of cotton. Nature Communications 8:15588

    doi: 10.1038/ncomms15588

    CrossRef   Google Scholar

    [25]

    Streubel J, Pesce C, Hutin M, Koebnik R, Boch J, et al. 2013. Five phylogenetically close rice SWEET genes confer TAL effector-mediated susceptibility to Xanthomonas oryzae pv. oryzae. New Phytologist 200:808−19

    doi: 10.1111/nph.12411

    CrossRef   Google Scholar

    [26]

    Tamura K, Stecher G, Peterson D, Filipski A, Kumar S. 2013. MEGA6: molecular evolutionary genetics analysis version 6.0. Molecular Biology and Evolution 30:2725−29

    doi: 10.1093/molbev/mst197

    CrossRef   Google Scholar

    [27]

    Chen HY, Huh JH, Yu YC, Ho LH, Chen LQ, et al. 2015. The Arabidopsis vacuolar sugar transporter SWEET2 limits carbon sequestration from roots and restricts Pythium infection. The Plant Journal 83:1046−58

    doi: 10.1111/tpj.12948

    CrossRef   Google Scholar

    [28]

    Cao X, Shen Q, Ma S, Liu L, Cheng J. 2020. Physiological and PIP transcriptional responses to progressive soil water deficit in three mulberry cultivars. Frontiers in Plant Science 11:1310

    doi: 10.3389/fpls.2020.01310

    CrossRef   Google Scholar

    [29]

    Liu L, Cao X, Zhai Z, Ma S, Tian Y, et al. 2022. Direct evidence of drought stress memory in mulberry from a physiological perspective: Antioxidative, osmotic and phytohormonal regulations. Plant Physiology and Biochemistry: PPB 186:76−87

    doi: 10.1016/j.plaphy.2022.07.001

    CrossRef   Google Scholar

    [30]

    Zhu P, Kou M, Liu C, Zhang S, Lü R, et al. 2021. Genome sequencing of Ciboria shiraiana provides insights into the pathogenic mechanisms of hypertrophy sorosis scleroteniosis. Molecular Plant-Microbe Interactions® 34:62−74

    doi: 10.1094/mpmi-07-20-0201-r

    CrossRef   Google Scholar

    [31]

    He N, Zhang C, Qi X, Zhao S, Tao Y, et al. 2013. Draft genome sequence of the mulberry tree Morus notabilis. Nature Communications 4:2445

    doi: 10.1038/ncomms3445

    CrossRef   Google Scholar

    [32]

    Jiao F, Luo R, Dai X, Liu H, Yu G, et al. 2020. Chromosome-level reference genome and population genomic analysis provide insights into the evolution and improvement of domesticated mulberry (Morus alba). Molecular Plant 13:1001−12

    doi: 10.1016/j.molp.2020.05.005

    CrossRef   Google Scholar

    [33]

    Xia Z, Dai X, Fan W, Liu C, Zhang M, et al. 2022. Chromosome-level genomes reveal the genetic basis of descending dysploidy and sex determination in Morus plants. Genomics, Proteomics & Bioinformatics In Press

    doi: 10.1016/j.gpb.2022.08.005

    CrossRef   Google Scholar

    [34]

    Chao N, Huang S, Kang X, Yidilisi K, Dai M, et al. 2022. Systematic functional characterization of cinnamyl alcohol dehydrogenase family members revealed their functional divergence in lignin biosynthesis and stress responses in mulberry. Plant Physiology and Biochemistry 186:145−56

    doi: 10.1016/j.plaphy.2022.07.008

    CrossRef   Google Scholar

    [35]

    Chen C, Chen H, Zhang Y, Thomas HR, Frank MH, et al. 2020. TBtools: an integrative toolkit developed for interactive analyses of big biological data. Molecular Plant 13:1194−202

    doi: 10.1016/j.molp.2020.06.009

    CrossRef   Google Scholar

    [36]

    Yu CS, Chen YC, Lu CH, Hwang JK. 2006. Prediction of protein subcellular localization. Proteins: Structure, Function, and Bioinformatics 64:643−51

    doi: 10.1002/prot.21018

    CrossRef   Google Scholar

    [37]

    Camacho C, Coulouris G, Avagyan V, Ma N, Papadopoulos J, et al. 2009. BLAST+: architecture and applications. BMC Bioinformatics 10:421

    doi: 10.1186/1471-2105-10-421

    CrossRef   Google Scholar

    [38]

    Wang Y, Tang H, DeBarry JD, Tan X, Li J, et al. 2012. MCScanX: a toolkit for detection and evolutionary analysis of gene synteny and collinearity. Nucleic Acids Research 40:e49

    doi: 10.1093/nar/gkr1293

    CrossRef   Google Scholar

    [39]

    Chao N, Wang R, Hou C, Yu T, Miao K, et al. 2021. Functional characterization of two Chalcone isomerase (CHI) revealing their responsibility for anthocyanins accumulation in mulberry. Plant Physiology and Biochemistry 161:65−73

    doi: 10.1016/j.plaphy.2021.01.044

    CrossRef   Google Scholar

    [40]

    Shukla P, Reddy RA, Ponnuvel KM, Rohela GK, Shabnam AA, et al. 2019. Selection of suitable reference genes for quantitative real-time PCR gene expression analysis in Mulberry (Morus alba L.) under different abiotic stresses. Molecular Biology Reports 46:1809−17

    doi: 10.1007/s11033-019-04631-y

    CrossRef   Google Scholar

    [41]

    Yan P, Zeng Y, Shen W, Tuo D, Li X, et al. 2020. Nimble cloning: a simple, versatile, and efficient system for standardized molecular cloning. Frontiers in Bioengineering and Biotechnology 7:460

    doi: 10.3389/fbioe.2019.00460[PubMed

    CrossRef   Google Scholar

    [42]

    Li R, Liu L, Dominic K, Wang T, Fan T, et al. 2018. Mulberry (Morus alba) MmSK gene enhances tolerance to drought stress in transgenic mulberry. Plant Physiology and Biochemistry 132:603−11

    doi: 10.1016/j.plaphy.2018.10.007

    CrossRef   Google Scholar

    [43]

    Lv Z, Hao L, Ma B, He Z, Luo Y, et al. 2021. Ciboria carunculoides suppresses mulberry immune responses through regulation of salicylic acid signaling. Frontiers in Plant Science 12:658590

    doi: 10.3389/fpls.2021.658590

    CrossRef   Google Scholar

    [44]

    Zhu P, Zhang S, Li R, Liu C, Fan W, et al. 2021. Host-induced gene silencing of a G protein α subunit gene CsGpa1 involved in pathogen appressoria formation and virulence improves tobacco resistance to Ciboria shiraiana. Journal of Fungi 7:1053

    doi: 10.3390/jof7121053

    CrossRef   Google Scholar

    [45]

    Xuan C, Lan G, Si F, Zeng Z, Wang C, et al. 2021. Systematic genome-wide study and expression analysis of SWEET gene family: sugar transporter family contributes to biotic and abiotic stimuli in watermelon. International Journal of Molecular Sciences 22:8407

    doi: 10.3390/ijms22168407

    CrossRef   Google Scholar

    [46]

    Geng Y, Wu M, Zhang C. 2020. Sugar transporter ZjSWEET2.2 mediates sugar loading in leaves of Ziziphus jujuba mill. Frontiers in Plant Science 11:1081

    doi: 10.3389/fpls.2020.01081

    CrossRef   Google Scholar

    [47]

    Pommerrenig B, Müdsam C, Kischka D, Neuhaus HE. 2020. Treat and trick: common regulation and manipulation of sugar transporters during sink establishment by the plant and the pathogen. Journal of Experimental Botany 71:3930−40

    doi: 10.1093/jxb/eraa168

    CrossRef   Google Scholar

    [48]

    Chen L, Lin IW, Qu X, Sosso D, McFarlane HE, et al. 2015. A cascade of sequentially expressed sucrose transporters in the seed coat and endosperm provides nutrition for the Arabidopsis embryo. The Plant Cell 27:607−19

    doi: 10.1105/tpc.114.134585

    CrossRef   Google Scholar

    [49]

    Gao Y, Zhang C, Han X, Wang Z, Ma L, et al. 2018. Inhibition of OsSWEET11 function in mesophyll cells improves resistance of rice to sheath blight disease. Molecular Plant Pathology 19:2149−61

    doi: 10.1111/mpp.12689

    CrossRef   Google Scholar

    [50]

    Gebauer P, Korn M, Engelsdorf T, Sonnewald U, Koch C, et al. 2017. Sugar accumulation in leaves of Arabidopsis sweet11/sweet12 double mutants enhances priming of the salicylic acid-mediated defense response. Frontiers in Plant Science 8:1378

    doi: 10.3389/fpls.2017.01378

    CrossRef   Google Scholar

    [51]

    Li Y, Wang Y, Zhang H, Zhang Q, Zhai H, et al. 2017. The plasma membrane-localized sucrose transporter IbSWEET10 contributes to the resistance of sweet potato to Fusarium oxysporum. Frontiers in Plant Science 8:197

    doi: 10.3389/fpls.2017.00197

    CrossRef   Google Scholar

  • Cite this article

    Kang X, Huang S, Feng Y, Fu R, Tang F, et al. 2023. SWEET transporters and their potential roles in response to abiotic and biotic stresses in mulberry. Beverage Plant Research 3:6 doi: 10.48130/BPR-2023-0006
    Kang X, Huang S, Feng Y, Fu R, Tang F, et al. 2023. SWEET transporters and their potential roles in response to abiotic and biotic stresses in mulberry. Beverage Plant Research 3:6 doi: 10.48130/BPR-2023-0006

Figures(8)  /  Tables(1)

Article Metrics

Article views(6362) PDF downloads(925)

ARTICLE   Open Access    

SWEET transporters and their potential roles in response to abiotic and biotic stresses in mulberry

Beverage Plant Research  3 Article number: 6  (2023)  |  Cite this article

Abstract: Mulberry (Morus spp., Moraceae) is a traditional economic crop plant and is also being gradually utilized as a beverage plant. SWEETs (Sugars Will Eventually be Exported Transporter) are important sugar transporters involved in various biological processes and responses to various stresses. However, SWEETs in mulberry are still poorly studied without a comprehensive functional analysis of SWEETs. In the present study, a total of 24 SWEETs were identified using the Morus alba (Ma) genome. Phylogenetic analysis showed that these 24 MaSWEETs were clustered with SWEETs from Arabidopsis, Populus and Oryza and fell into four clades. MaSWEETs in the same clade are likely to pose similar intron/exon patterns. These MaSWEETs distributed on 12 chromosomes and tandem duplication and block duplication were responsible for the expansion of SWEETs in mulberry. Transmembrane domains and conserved active sites of Tyr and Asp were observed in MaSWEETs. Cis-elements in promoter regions of MaSWEETs indicated the possible function of MaSWEETs in response to hormones and environment stimulus. MaSWEETs showed quite different expression preference in tissues and organs indicating the possible function divergence. In addition, most MaSWEETs showed a disturbed expression levels in response to various abiotic stresses and Ciboria shiraiana infection. MaSWEET1a was functionally characterized as a negative regulator of resistance to C. shiraiana infection based on in vivo transient overexpression of MaSWEET1a in tobacco and down-regulation of MaSWEET1a/b in mulberry. Our results provided foundation for further functional dissection of SWEETs in mulberry and a potential regulator for genetic modification.

    • Sugars are predominant carbon and energy sources and support plant vegetative and reproductive growth[1]. Transport of sugars across the plant bio-membrane needs the assistance of specific transporters[1,2]. These transporters act as bridges that mediate the distribution of sugars between source–sink organs, which is critical for sugar homeostasis and the cellular exchange of sugar efflux in multicellular organisms[1,3,4].

      SWEETs (Sugars Will Eventually be Exported Transporter) and SUTs (sucrose transporters), MSTs (monosaccharide transporters) are the main known sugar transporters in eukaryotes[5] . Unlike SUTs and MSTs, the relatively newly reported sugar transporter SWEETs are pH-independent transporters. SWEETs play important roles in phloem transport and act as bidirectional transmembrane transporters of sugars along the concentration gradient[3,4]. AtSWEET1 was first identified as a glucose transporter with clear functional characterization[2]. In addition, the SWEET multi-gene family was identified and classified into four clades and the functional divergence of these paralogs were also revealed in Arabidopsis at the same time[2]. For example, Clade II AtSWEET8 contributes to pollen viability and Clade III AtSWEET15 is involved in leaf senescence. Thereafter, another SWEET in Arabidopsis Clade III AtSWEET9 was characterized as an important transporter involved in nectar secretion[6]. Besides the SWEET family in Arabidopsis, the SWEET gene family has been identified in many plants including Tea (Camellia sinensis, Cs), tomato (Solanum lycopersicum, Sl), wheat (Triticum aestivum, Ta) barrel medic (Medicago truncatula, Mt), cabbage (Brassica rapa, Br), daylily (Hemerocallis fulva, Hf), grapevine (Vitis vinifera, Vv), rice (Oryza sativa, Os) and poplar (Populus trichocarpa, Pt and P. alba × P. glandulosa, Pag)[713]. These SWEET homologs belong to the MtN3/saliva family and consist of seven α-helical transmembrane domains (TMs): a tandem repeat of three transmembrane domains (TMs) connected with a linker-inversion TM[2,14].

      SWEETs participate in various biological processes including development, flowering, stress responses and plant-pathogen interaction in plants[4,15]. In addition, different SWEET gene family members show functional divergence or redundancy. In Arabidopsis, AtSWEET8 and 13 support pollen development; AtSWEET11 and 12 provide sucrose to the SUTs for phloem loading and play distinct roles in seed filling; and AtSWEET9 is essential for nectar secretion[24,6,16]. BrSWEET9 in Brassica rapa was also reported to be involved in nectar secretion[17]. Overexpression of PagSWEET7 promotes secondary growth and xylem sugar content[12]. OsSWEET11 and 15 have functions affecting pollen development and are key players in seed filling in rice[18]. SWEET homologs also play important roles in abiotic stress responses. Overexpression of AtSWEET16 promotes freezing tolerance in Arabidopsis[19] while AtSWEET11 and 12 mutants exhibit greater freezing tolerance[20]. AtSWEET15 could be induced by various abiotic stresses including osmotic, drought, salinity, and cold stresses and overexpression of AtSWEET15 results in transgenic plants with hypersensitiveness to cold and salinity stresses[21]. CsSWEET1a and CsSWEET16 were also reported to mediate freezing tolerance[22,23]. Recently, more studies have shown that SWEETs are involved in plant-pathogen interaction and are known as susceptibility (S) genes, acting as targets of effector proteins during host–microbe interactions in many plant species[15]. GhSWEET10, is induced by Avrb6, a transcription activator-like (TAL) effectors from Xanthomonas citri subsp. Malvacearum (Xcm) and is responsible for maintaining virulence of Xcm avrb6 and the cotton susceptibility to infections[24]. OsSWEET11–15 belonging to clade III in rice have been shown to be induced by TAL effector from Xanthomonas oryzae and support pathogen growth[25]. In contrast, some SWEETs could also function as resistance genes. Overexpression of IbSWEET10 can promote resistance to F. oxysporum in sweet potato[26]. Mutation of AtSWEET2 resulted in increased susceptibility to the root necrotrophic pathogen Pythium irregulare[27].

      Mulberry (Morus spp., Moraceae) is a traditional economic crop plant and a new beverage plant. In addition, its fruits are rich in nutrient and bioactive components and the ripening process of mulberry fruits along with sugar accumulation and distribution. Mulberry suffers various abiotic stresses and the disasterous fungal disease sclerotiniose which bursts at the early stage of mulberry fruit development[2830]. Mulberry fruits with sclerotiniose lose their color and flavor and turn pale instead of ripening. C. shiraiana is the dominant causal agent of mulberry sclerotiniose in China, and it results in hypertrophy sorosis sclerotiniose. SWEETs as the important transporters involved in sugar homeostasis are expected to be involved in mulberry fruit development and interaction with sclerotiniose pathogens. However, to date, few studies on SWEETs have been reported in mulberry, although the SWEET gene family may play important roles in mulberry fruit development and responses to abiotic and biotic stresses. Mulberry genome information has been released successively since the Morus notabilis genome was reported in 2013[31]. The chromosome-level genome of M.alba (Ma) was released by Jiao et al. and the genome of M. yunnanensis was recently released by Xia et al.[32,33]. Released genome information makes it possible to perform genome-wide characterization of the SWEET gene family in mulberry. In the present study, a total of 24 SWEET genes were identified in the Morus alba genome and their phylogenetic classification, conserved motifs, gene structures, distribution on chromosomes, cis-elements in promoter regions and tissue expression profile were revealed. In addition, the responses of MaSWEETs to various abiotic stresses and sclerotiniose pathogen infection were also detected. MaSWEET1a was functionally characterized as a negative regulator which increased the mulberry susceptibility to C. carunculoides infection.

    • The xylem, phloem, fruits at four different developmental stages (S0, inflorescence; S1, green fruits; S2, reddish fruits; S3, purple fruits)) and diseased fruit infected with C. shiraiana of Morus atropurpurea variety Zhongshen 1 (Mazs) were collected from the National Mulberry Genebank (NMGB) in Zhenjiang, China for expression profiling. Seedlings of the M. alba var. Fengchi and tobacco (Nicotiana Benthamiana) were grown in a chamber at 22 °C with a 16/8 day/night cycle and 40%–60% humidity. C. shiraiana was provided by Professor Zhao and was cultured in potato dextrose agar (PDA) medium.

      Tobacco at the four-leaf stage was used for transient overexpression. M. alba var. Fengchi seedlings at the four-euphylla stage were used for virus-induced gene silencing (VIGS). Four-week-old seedlings with similar growth conditions (~12–15 cm high) were used for treatments under different abiotic stresses. Detailed information for abiotic stress treatments were reported in our previous study[34]. All the above samples were immediately frozen in liquid nitrogen after collection and then stored at −80 °C until use. Three biological replicates were performed for each experiment.

    • The M. alba genome sequences (.fasta) and annotation file (.gff) were generously provided by Professor Jiao, who released this genome information. The Hidden Markov Model (HMM) profiles of the SWEET domain (PF03083) were downloaded from the Pfam database (http://pfam.xfam.org/) and used to search the candidate SWEET proteins in the M. alba proteome with HMMER software. In addition, the protein sequences of AtSWEETs, OsSWEETs and PtSWEETs were downloaded from TAIR (www.arabidopsis.org/), TIGR (http://rice.plantbiology.msu.edu/) and phyto-zome (https://phytozome-next.jgi.doe.gov/) respectively, and used as queries to search against the M. alba proteome. The Toolbox for Biologists v1.098774[35] was used to analyze the sequence length, molecular weight and theoretical isoelectric point (pI) values of each MaSWEET protein. The distributions of TM helices were predicted by the TMHMM Server v. 2.0 (www.cbs.dtu.dk/services/TMHMM). Prediction of subcellular localization of MaSWEET proteins using the online Tool WoLF PSORT (www.genscript.com/wolf-psort.html)[36].

    • Chromosome location information of MaSWEETs was extracted based on the Morus alba genome annotation file. Tbtools v1.098774 were used to identify syntenic blocks and tandem duplication events using default parameters[37,38]. The results were visualized using Tbtools v1.098774 and both the tandem duplication and block duplication gene pairs were marked.

    • MaSWEETs were aligned using clustal W assembled in MEGA11.0. The alignment result was exported and manually speculated for scanning the MtN3 repeats. The online MEME Suite version 5.5.0 was used to identify 7 conserved motifs from 24 amino acid sequences of SWEET genes in Morus alba. The Hidden Markov Model (HMM) profiles of the SWEET domain (PF03083) were downloaded from the Pfam database.

    • The gene structure of each MaSWEET was displayed based on the genome sequence and its annotation file using Gene Structure View assembled in Tbtools v1.098774. The upstream 2000 bp sequences were extracted for in silico promoter region analysis. Cis-acting elements were predicted using PlantCARE (http://bioinformatics.psb.ugent.be/webtools/plantcare/html/).

    • A neighbor-joining (NJ) phylogenetic tree was constructed using full-length SWEETs protein sequences from A. thaliana, P. trichocarpa, O. sativa and M. alba using MEGA11.0[26] with JTT + G model and bootstrap test with 1000 replicates.

    • RNA extraction and cDNA synthesis were performed as in our previous report using Plant RN52 Kit (Aidlab, Beijing, China) and PC54-TRUEscript RT kit (Aidlab, Beijing, China) according to the manual[39]. RT-qPCR (quantitative real-time PCR) was performed to validate the expression patterns of MaSWEETs in different tissues, fruit development stages and stresses using ABI StepOnePlus™ Real-Time PCR System (USA). The primers are available in Supplemental Table S1. Actin was used as a reference gene[40]. Graphpad Prism8.0 was used to visualize the RT-qPCR results and to perform T-test and ANOVA. P value < 0.05 was marked as significant. At least three individuals were used and three technical replicates respectively were performed for RT-qPCR.

    • The recombinant plasmids pNC-1304-35S:MaSWEET1a were constructed using nimble cloning[41]. Both recombinant plasmids pNC-1304-35S: SWEET1a and empty vector pNC-Cam1304-35SMCS, as the negative control, were transformed into Agrobacterium tumefaciens GV3101 and then transferred into N. benthamiana leaves via Agrobacterium-mediated transient transformation, as previously reported[41]. Overexpression of MaSWEET1a was determined using RT-qPCR by comparing the expression levels of target genes in transgenic plants with those in the negative controls.

    • Virus-induced gene silencing (VIGS) was used to obtain MaSWEET1a/b down-regulated mulberry, in accordance with our previous report[42]. Agrobacterium tumefaciens containing recombinant plasmids pTRV2-MaSWEET1a/b, pTRV1 and pTRV2 (negative control) were cultured in transient transformation buffer and then transferred into mulberry leaves by means of pressure injection. Ten independent mulberry plants were injected. The knock-down efficiency for the target genes was determined by RT-qPCR 15 d after injection by comparing the transgenic plants with the negative controls and wild types.

    • Cell death symptoms and the growth condition of C. shiraiana were recorded to estimate the resistance of transgenic plants to C. shiraiana infection[43,44]. C. shiraiana was inoculated at 2 d after infiltration in tobacco and at 10 d after infiltration in mulberry. The cell death symptoms were photographed after inoculation until the sclerotia appeared. The results are representative of at least three biological replicates.

    • A total of 24 SWEET homologs were identified based on the genome information of M. alba and named according to their orthologs in A. thaliana, P. trichocarpa or V. vinifera[45]. These MaSWEETs encode proteins ranging from 197 aa to 304 aa with molecular weight from 21.45 to 34.12 kDa and theoretical isoelectric points from 7.16 to 9.61 (Table 1). Subcellular localization prediction of these MaSWEETs showed that most of them (18/24) distributed on membrane structures such as plasma membrane (PM), tonoplast membrane (TM) and chloroplast thylakoid membrane (CTM). Phylogenetic analysis of MaSWEETs and SWEETs from model plants such as A. thaliana, P. trichocarpa and O. sativa showed that four clades were formed by these SWEETs (Fig. 1). According to previous studies, SWEETs in plants were generally classed as four phylogenetic clades which is in agreement with our results[8]. Major MaSWEETs (10/24) together with AtSWEET1, 2 and 3 belong to clade I. Clade II and IV contain five MaSWEETs each and Clade III contains four MaSWEETs (Fig. 1 and Table 1).

      Table 1.  SWEET gene family in Morus alba.

      CladeGene nameAccession no.Gene IDCDS Size
      (bp)
      Protein physicochemical characteristicsTMHsSubcellularMtN3/Saliva (PQ-Loop Repeat)
      Length (aa)MW (kDa)pIAliphatic indexLocalization*Domain Position
      IMaSWEET1aXM_024170697.1-0M.alba_G001204972924226.559.61113.517CTM6-94, 131-209
      IMaSWEET1bXM_024170698.1-0M.alba_G001204967822521.459.51115.746CTM1-49, 86-164
      IMaSWEET2aXM_024163709.1-0M.alba_G001924470823525.948.39129.797CTM18-104, 137-221
      IMaSWEET2bXM_024164777.1-0M.alba_G001086370523425.928.98122.867TM17-101, 137-218
      IMaSWEET2cXM_024163707.1-0M.alba_G001924477425728.588.58132.337PM54-126, 159-243
      IMaSWEET2dXM_024163712.1-0M.alba_G001924468122625.278.22128.896EX23-95, 128-212
      IMaSWEET2eXM_024163708.1-0M.alba_G001924472924226.728.49128.067CTM39-111, 144-228
      IMaSWEET2fXM_024163703.1-0M.alba_G001924477725828.88.49132.568PM41-127, 160-244
      IMaSWEET2gXM_024163711.1-0M.alba_G001924468422725.497.62129.167PM10-96, 129-213
      IMaSWEET3XM_010099554.2-0M.alba_G000306378326029.018.89115.737PM9-98, 132-216
      IIMaSWEET4aXM_010091939.1-0M.alba_G000927673524427.459.28109.397CTM10-98, 134-218
      IIMaSWEET4bXM_010113461.2-0M.alba_G000153673824527.48.98122.087EX11-97, 134-216
      IIMaSWEET5XM_024168739.1-0M.alba_G001829571123626.647.63120.937CY10-93, 131-127
      IIMaSWEET7aXM_010108966.2-0M.alba_G000511077425728.329.57128.567CTM11-95, 134-218
      IIMaSWEET7bXM_010108964.1-0M.alba_G000510978926228.999124.967CTM10-97, 134-218
      IIIMaSWEET10XM_010095631.2-0M.alba_G001801688829533.128.86120.317CTM11-96, 132-216
      IIIMaSWEET11aXM_010114440.2-0M.alba_G001690180426729.639.47124.087CTM11-99, 135-220
      IIIMaSWEET11bXM_010095633.2-0M.alba_G001801591530434.127.57112.897CTM12-99, 134-219
      IIIMaSWEET15XM_010092381.2-0M.alba_G000676788529433.427.16109.017CTM12-99, 133-219
      IVMaSWEET16XM_024167733.1-0M.alba_G001461790930233.29.08114.277CTM20-92, 129-212
      IVMaSWEET17aXM_024171451.1-0M.alba_G001461470823526.468.71119.875CY5-78, 116-198
      IVMaSWEET17bXM_024167902.1-0M.alba_G001461375325028.058.71120.886PM8-93, 131-213
      IVMaSWEET17cXM_024171286.1-0M.alba_G000819572023926.728.94111.727CY6-92, 129-212
      IVMaSWEET17dXM_024171287.1-1M.alba_G0008196723240279.43122.677CTM6-92, 127-213
      * The subcellular localizations were predicted by WoLFPSORT. PM, plasma membrane; EX, extracellular; CY, cytoplasmic; TM, tonoplast membrane; CTM, chloroplast thylakoid membrane.

      Figure 1. 

      Phylogenetic relationships of the SWEET family genes in Arabidopsis, Oryza sativa, Populus, Vitis vinifera, and Morus alba. The sequences of the 104 SWEET proteins from the above four plant species were aligned by Clustal Omega, and the phylogenetic tree was constructed by the MEGA 11.0 using the NJ method with 1000 bootstrap replicates. The proteins from Arabidopsis, Oryza sativa, Populus, Vitis vinifera, and Morus alba are indicated with the prefixes of At, Os, Pt, Vv, and Ma, respectively.

    • MaSWEETs distributed on 12 chromosomes except chromosome 1 and 3. Chromosome 5 occupied five MaSEETs which formed a gene cluster. Chromosome 5 is the chromosome that had the most SWEETs and the following is chromosome 6 and 8 which had four MaWSEETs each. There was only one MaSWEET locating on chromosome 4 (Fig. 2). In addition, there were three MaSWEETs on chromosome 2 and 7 respectively and two MaSWEETs on chromosome 1 and 5 respectively. Gene duplication including block duplication and tandem duplication is the main cause for gene family expansion. Tandem duplications were found on chromosome 2, 6 and 12 (linked by red lines in Fig. 2). It is also interesting to find several possible gene clusters such as MaSWEET1a-b on chromosome 5, MaSWEET2a-g on chromosome 9, MaSWEET7 a-b on chromosome 13, MaSWEET17a-b on chromosome 12 and MaSWEET17c-d on chromosome 6. Two gene pairs (MaSWEET4b/MaSWEET5 and MaSWEET11b/ MaSWEET15) resulting from block duplications were also marked (linked by black lines in Fig. 2).

      Figure 2. 

      Distribution of MaSWEET genes in Morus alba chromosomes. The tandem gene pairs are linked by red lines. The block duplications gene pairs are marked by black lines. The scale is provided in megabase (Mb).

    • MaSWEETs always located on membrane structures with transmembrane domains. The prediction results of MaSWEETs using DeepTMHMM showed that most MaSWEETs posed seven types of transmembrane helices (TMH) named TMH1-7 (Table 1, Supplemental Fig. S1). Alignment and conserved motif analysis showed that almost all MaSWEETs kept the conserved TMH and active sites Tyr and Asp indicating by red full triangles (Fig. 3). The active residues Tyr and Asp were reported to be involved in forming hydrogen bonds to ensure sugar transport activity[14]. In addition, all MaSWEETs except MaSWEET4a had a conserved Ser in each triple helix bundle (THB) which can be phosphorylated and is important for SWEET activity (Fig. 3). MaSWEET4a replaced Ser with Thr at the first Ser phosphorylation site between TMH1 and TMH2 which may also retain similar activity as Thr is also a common phosphorylation site. All MaSWEETs retained the conserved second Ser phosphorylation site between TMH5 and TMH6 (Fig. 3).

      Figure 3. 

      Multiple sequence alignment of MaSWEET proteins. The positions of the TMHs are underlined. The positions of the active sites of tyrosine (Y) and aspartic acid (D) are indicated by red triangles. The conserved serine (S) phosphorylation sites are indicated by blue triangles.

    • MaSWEETs gene structures were identified based on the annotation information of the M. alba genome. In summary, there were six MaSWEETs with six introns, 12 MaSWEETs with five introns and five MaSWEETs with four introns (Fig. 4). In addition, genes clustered together based on phylogenetic analysis are likely to show similar gene structures and length. For example, MaSWEET2a, c, d, e, f and g with six introns, MaSWEET7a and MaSWEET7b with four introns, and MaSWEET17a and MaSWEET17b with four introns.

      Figure 4. 

      Gene organization of MaSWEETs and cis-elements in promoter regions of MaSWEETs. (a) Phylogenetic tree using 24 MaSWEETs. (b) Exon/intron structures of Morus alba L. SWEETs. (c) Cis-element distribution in the promoter regions of MaSWEETs.

      Promoter region analysis of MaSWEETs indicated the possible function of MaSWEETs in response to hormones and environment stimulus. Among all the 22 types of cis-elements identified, most of them are light response elements accounting for 44% of the total elements (Supplemental Table S2). In addition, hormone response elements were also widely identified in the promoters of MaSWEETs (Fig. 4c). Most MaSWEETs had abscisic acid (ABA), salicylic acid (SA) or methyl jasmonate (MeJA) related response elements in their promoter regions. Especially, MaSWEET1a-b, MaSWEET16, MaSWEET17a-d had cis-elements involved in response to five types of hormones (ABA, SA, MeJA, auxin and gibberellins). Several Myb binding cis-elements were also identified in promoter regions of MaSWEET2a and MaSWEET10 (Fig. 4c, Supplemental Table S2).

    • The tissue or organ expression profiles of MaSWEETs were revealed. The MaSWEETs with high sequence identity (> 91%) are hard to distinguished by RT-qPCR and were determined by common primers to reveal their total transcription levels. MaSWEETs in phylogenetic clade I showed quite similar expression patterns with highest expression levels in leaf and relatively higher expression levels in early stages (S0 and S1) of fruit development except MaSWEET1a/b (Fig. 5ad). MaSWEET1a/b showed higher expression levels in fruits during whole fruit development with highest expression level in fruit at S0 stage (Fig. 5a). However, MaSWEETs (MaSWEET10, 11a-b, and 15) in phylogenetic clade III showed preferential expression in fruits especially at the late stages (S2 and S3) (Fig. 5jm). MaSWEETs in phylogenetic clade II had different expression patterns in different tissues or organs. MaSWEET4a and b showed highest expression level in fruit at the S1 stage while MaSWEET5 showed obviously preferential expression in xylem (Fig. 5eg). MaSWEET7a and b showed similar expression pattern with MaSWEET2 cluster and MaSWEET3 from phylogenetic clade I. MaSWEETs in phylogenetic clade IV showed similar expression pattern with highest expression levels in leaf (Fig. 5b, d, h, i). MaSWEET16 also had higher expression in fruit with similar expression level at four different development stages (Fig. 5n).

      Figure 5. 

      Transcript levels of MaSWEETs in leaves, xylem, phloem, and different development stages of fruit. Three technical replicates were analyzed. Error bars represent SE. Different letters indicate statistically significant differences (Duncan's test, p < 0.05).

    • Mulberry sclerotiniose is a fungal disease resulting from fungal pathogen infection. Most (20/24) MaSWEETs showed positive or negative responses to the fungal infection. MaSWEET1a/b, MaSWEET2 cluster, MaSWEET4b, and MaSWEET17 a-d showed a significant decrease of expression levels in diseased fruits with sclerotiniose compared with the expression levels in healthy fruits (Fig. 6). In contrast, MaSWEET2b, MaSWEET3, MaSWEET7b, MaSWEET10 and MaSWEET11a-b showed significant increases of expression levels in diseased fruits. MaSWEETs also played roles in response to various abiotic stresses including drought, water logging, cold and high temperature. MaSWEET1a/b showed a positive response to drought with significant increasing expression levels while other clade I MaSWEETs, MaSWEET2b and 3 significantly decreased their expression level under detected abiotic stresses (Fig. 7ad). In contrast, MaSWEET16 significantly increased its expression level under detected abiotic stresses. MaSWEET4a-b, MaSWEET5 in phylogenetic clade II and MaSWEET11a-b in clade III showed similar response patterns with a significant increase of expression levels in response to low temperature (4 °C), high temperature (40 °C) or drought (Fig. 7eg). MaSWEET15 showed high sensitivity for drought and significant increase of expression levels under drought stress. MaSWEET17a-d showed a negative response to temperature change with a significant decrease of expression levels under low temperature and high temperature treatments (Fig. 7oq).

      Figure 6. 

      Expression levels of 24 selected MaSWEET genes in response to fungi stress conditions. Three technical replicates were analyzed. Error bars represent SE. Asterisks indicate significant difference as determined by Student’s t-test (* p < 0.05; ** p < 0.01; *** p < 0.001; **** p < 0.0001).

      Figure 7. 

      Expression levels of 24 selected MaSWEET genes in response to 4 °C, 40 °C , drought, and flood stress conditions. Three technical replicates were analyzed. Error bars represent SE. Asterisks indicate significant difference as determined by Student’s t-test (* p < 0.05 ; ** p < 0.01 ; *** p < 0.001 ; **** p < 0.0001 ).

    • MaSWEET1a/b is quite different from other MaSWEETs in clade I based on expression profile analysis, which showed preferential expression in fruits and a negative response to sclerotiniose pathogen (Ciboria shiraiana) infection. Our unpublished data indicated MaSWEET1a as key genes involved in the pathogen infection process based on comparative transcriptome analysis. Transient overexpression of MaSWEET1a in tobacco and VIGS knock-down of MaSWEET1a/b in mulberry were performed. RT-qPCR results validated the successful overexpression of MaSWEET1a in tobacco and knock-down of MaSWEET1a/b in mulberry (Fig. 8b, d). The expression level of MaSWEET1a affected the resistance to C. shiraiana infection in both tobacco and mulberry (Fig. 8a, d). Overexpression of MaSWEET1a decreased the resistance to C. shiraiana infection with more severe cell death symptoms observed in OE-line tobacco (Fig. 8a). Knock-down of MaSWEET1a/b in mulberry could increase the resistance to C. shiraiana infection in mulberry (Fig. 8d). These results proved that MaSWEET1a is an important negative regulator of resistance to C. shiraiana infection in mulberry.

      Figure 8. 

      Functional characterization of MaSWEET1a/b in tobacco and mulberry. (a) Damage of tobacco overexpressing MaSWEET1a after infection by C. shiraiana. (b) Expression levels of MaSWEET1a in tobacco. MaSWEET1a-T1, T2, T3 are three independent treatments with transient overexpressed MaSWEET1a. (c) Damage of mulberry with knock-down MaSWEET1a/b after infection by C. shiraiana. (d) The expression level of MaSWEET1a/b in mulberry. MaSWEET1a/b-T1, T3, T4 are three independent treatments with down-regulation of MaSWEET1a/b using VIGS.

    • Functional studies on SWEETs have revealed that SWEET homologs not only act as loading and unloading transporters of sugars but also play critical roles in various biological processes. Typically, angiosperm genomes contain about 15–25 SWEET genes[4]. In the present study, a total of 24 SWEET genes were identified and clustered into four clades corresponding with the knowledge of SWEETs in angiosperm. Several SWEETs including MaSWEET1a/b, MaSWEET4a-b, MaSWEET10, MaSWEET11a-b and MaSWEET15 showed preferential expression in fruits indicating their possible roles in fruit development. Similar expression preference of AtSWEETs was also reported in Arabidopsis[2]. It is noted that fruit-preferential expressed MaSWEETs still showed temporal expression difference during fruit development indicating time-course regulation of MaSWEETs for fruit ripening in mulberry. Early-stage expressed MaSWEET1a/b and MaSWEET4a/b further differed from late-stage expressed MaSWEET10, MaSWEET11a/b and MaSWEET15 in terms of detailed expression patterns during fruit ripening. MaSWEET2b and 2 cluster genes (MaSWEET2a, c-g) showed preferential expression in leaves which is similar with the ortholog ZjSWEET2.2 in Ziziphus jujuba. ZjSWEET2.2 was reported to be involved in mediating sugar loading in leaves[46]. MaSWEET3, 7a-b,16 and 17a-d also showed highest expression levels in leaves indicating their possible roles in sugar source loading or unloading. MaSWEET11b showed higher expression levels in phloem and its ortholog AtSWEET11 in Arabidopsis was reported to be involved in sugar phloem loading[3].

      Sugar signal is critical for plants in response to various stresses. Previous studies have shown that SWEETs participated in abiotic and biotic responses in many plant species including arabidopsis and rice[21, 25]. AtSWEET11, 12, 15 and 16 were reported to be involved in affecting cold tolerance in Arabidopsis[1921]. HfSWEET17 was also reported as a positive regulator of resistance to cold stress in daylily[11]. Cold environment (4 °C) induced expression of MaSWEET4a, 4b, 5, 11a, 11b and 16 in mulberry. Interestingly, these cold-induced MaSWEETs can also be induced by high temperature. MaSWEET15 which is the ortholog of AtSWEET15 can be induced by drought as well as low or high temperature. MaSWEET15, MaSWEET1a/b. 4a, 4b, 5, 7a, 11a, 11b, and 16 also showed positive responses to drought. It is obvious that some SWEET genes can be induced by different stresses. AtSWEET15 was also reported to be induced by osmotic, drought and salinity[21]. MaSWEET4a, 4b, 11a, 11b and 16 can be induced by low or high temperature and drought indicating their important roles in response to various abiotic stresses in mulberry.

      SWEETs were generally thought to 'support the enemy' during infection. SWEETs especially those that function as exporters generally facilitate the export of sugars out of host cells, which support pathogen growth in the apoplasm[15,47,48]. Clade III SWEETs including AtSWEET11, 12, OsSWEET11 were characterized as negative regulators of resistance to fungal infection and Clade III SWEETs including OsSWEET11, 13, 14 and GhSWEET10 were characterized as negative regulators of resistance to bacterial pathogen infection[2,24,49,50]. In contrast, clade I AtSWEET2, a glucose importer and clade III IbSWEET10 were reported as positive regulators of resistance to fungal infection[48,51]. Therefore, roles of SWEETs in response to pathogen infection may be quite different. Most MaSWEETs were disturbed in diseased fruits that resulted from sclerotiniose pathogen infection. MaSWEET2b, MaSWEET3 in clade I, MaSWEET7b in clade II, MaSWEET10 and MaSWEET11a-b in clade III showed significant increase in expression levels in diseased fruits while MaSWEET1a/b, MaSWEET2 cluster (MaSWEET2a, c-g) in clade I showed significant decrease of expression levels in diseased fruits. MaSWEET1a was further validated as a negative regulator of resistance to C. shiraiana infection. Given the fact that MaSWEET1a/b was repressed in diseased fruits, it is likely that a possible pathway through repression of MaSWEET1a exists in mulberry to defense pathogen infection.

      In conclusion, we have performed a genome-wide investigation of SWEET genes in Morus and a comprehensive analysis including phylogenetic analysis, promoter analysis and expression profile analysis was also carried out. Their possible roles in development and response to abiotic and biotic stresses were addressed. In particular, the functional role of MaSWEET1a in regulation of tolerance to C. shiraiana infection was validated using both VIGS knock-down and transient overexpression in tobacco combined with inoculation of C. shiraiana. The results in this study provides a foundation for studying the function of the SWEET family in mulberry plants and provides a negative regulator of resistance to C. shiraiana infection for further genetic modification.

      • This work was jointly supported by the National Natural Science Foundation of China (32201526), Crop Germplasm Resources Protection Project of the Ministry of Agriculture and Rural Affairs of the People's Republic of China (19200382), National Infrastructure for Crop Germplasm Resources (NCGRC-2020-041), and China Agriculture Research System of MOF and MARA (CARS-18). We thank Professor Aichun Zhao at Southwest University for providing us C. shiraiana strains and Professor Feng Jiao at Northwest University of Agriculture and Forestry who provided us the genome annotation file of M. alba.

      • The authors declare that they have no conflict of interest.

      • # These authors contributed equally: Xiaoru Kang, Shuai Huang, Yuming Feng

      • Copyright: © 2023 by the author(s). Published by Maximum Academic Press, Fayetteville, GA. This article is an open access article distributed under Creative Commons Attribution License (CC BY 4.0), visit https://creativecommons.org/licenses/by/4.0/.
    Figure (8)  Table (1) References (51)
  • About this article
    Cite this article
    Kang X, Huang S, Feng Y, Fu R, Tang F, et al. 2023. SWEET transporters and their potential roles in response to abiotic and biotic stresses in mulberry. Beverage Plant Research 3:6 doi: 10.48130/BPR-2023-0006
    Kang X, Huang S, Feng Y, Fu R, Tang F, et al. 2023. SWEET transporters and their potential roles in response to abiotic and biotic stresses in mulberry. Beverage Plant Research 3:6 doi: 10.48130/BPR-2023-0006

Catalog

  • About this article

/

DownLoad:  Full-Size Img  PowerPoint
Return
Return