Transcriptome Profiling and Functional Validation of RING-Type E3 Ligases in Halophyte Sesuvium verrucosum under Salinity Stress

Owing to their sessile nature, plants have developed a tapestry of molecular and physiological mechanisms to overcome diverse environmental challenges, including abiotic stresses. Adaptive radiation in certain lineages, such as Aizoaceae, enable their success in colonizing arid regions and is driven by evolutionary selection. Sesuvium verrucosum (commonly known as Western sea-purslane) is a highly salt-tolerant succulent halophyte belonging to the Aizoaceae family; thus, it provides us with the model-platform for studying plant adaptation to salt stress. Various transcriptional and translational mechanisms are employed by plants to cope with salt stress. One of the systems, namely, ubiquitin-mediated post-translational modification, plays a vital role in plant tolerance to abiotic stress and other biological process. E3 ligase plays a central role in target recognition and protein specificity in ubiquitin-mediated protein degradation. Here, we characterize E3 ligases in Sesuvium verrucosum from transcriptome analysis of roots in response to salinity stress. Our de novo transcriptome assembly results in 131,454 transcripts, and the completeness of transcriptome was confirmed by BUSCO analysis (99.3% of predicted plant-specific ortholog genes). Positive selection analysis shows 101 gene families under selection; these families are enriched for abiotic stress (e.g., osmotic and salt) responses and proteasomal ubiquitin-dependent protein catabolic processes. In total, 433 E3 ligase transcripts were identified in S. verrucosum; among these transcripts, single RING-type classes were more abundant compared to multi-subunit RING-type E3 ligases. Additionally, we compared the number of single RING-finger E3 ligases with ten different plant species, which confirmed the abundance of single RING-type E3 ligases in different plant species. In addition, differential expression analysis showed significant changes in 13 single RING-type E3 ligases (p-value < 0.05) under salinity stress. Furthermore, the functions of the selected E3 ligases genes (12 genes) were confirmed by yeast assay. Among them, nine genes conferred salt tolerance in transgenic yeast. This functional assay supports the possible involvement of these E3 ligase in salinity stress. Our results lay a foundation for translational research in glycophytes to develop stress tolerant crops.


RNA-Seq Quality Check, De Novo Assembly, and Annotation
We generated~307 million 101 bp PE RNA-Seq reads with an Illumina NovaSeq 6000 from Sesuvium verrucosum root tissue that constitutes three control samples, three salt treatment early samples, and three salt treatment late samples (Table S1). More than 97% of the sequencing reads had a Phred score > Q20. After rRNA contamination removal, adapter trimming, and low-quality read trimming, we retained 233 million paired-end (PE) reads. De novo assembly of the trimmed reads yielded 301,627 transcripts (GC%: 43.1; N50: 2373) ( Table S2). The assembly quality of the transcriptome was further confirmed using BUSCO (~99.5% of single copy plant-specific orthologous genes) and read alignment back to assembled transcripts (98%). Further, assembled transcripts were clustered and 195,255 transcripts were generated (Table S2). By similarity search (BLAST) of these clustered transcripts against the UniPrto-trEMBL protein database, the bacterial and fungal related transcripts were removed manually from the cluster. Blob plots were generated to highlight the contaminants found in the assembly ( Figure S1). After removing the contamination, we retained 131,454 final transcripts (total size~207 Mb, N50 = 2531, G + C% = 39.5, length range= 295-29,964 bp, average transcript length = 1579 bp). The BUSCO analysis predicted~99.3% (single copy plant-specific orthologous genes) expressed genes from the final assembly. The read alignment of individual samples showed that 92% of reads aligned back to the final transcriptome (Table 1). A second quality check was performed by extracting the average read depth relative to transcript length ( Figure 1A) to ensure that the majority of transcripts were represented. The transcripts were similarity searched against UniProt-trEMBL, UniProt-Plant, and Arabidopsis thaliana proteome databases; 69,826 (~53%), 69,816 (~53%), and 62,922 (~48%) transcripts were annotated in S. verrucosum root transcriptome from each of the respective databases (Table S2). S. verrucosum transcripts showed high similarity (as observed from top BLAST hits) with Spinacia oleracea (22,830), Opuntia streptacantha (21,804), Beta vulgaris subsp. Vulgaris (4109), and Vitis vinifera (2296) ( Figure 1B). The relationships between the identified top hit plant species are shown in Figure 1C. We were able to identify Gene Ontology (GO) IDs in 53,575 (~41%) assembled transcripts. We found 42,432 Biological Process (BP) GO ID's, 77,656 Molecular Function (MF) GO IDs, and 45,556 Cellular Compartment (CC) GO IDs from the transcriptome (Table S3). The most abundant GO IDs found in the transcriptome are shown in Figure 1D. Average

Positive Selection and Ortholog-Based Analysis
We performed an orthology analysis of the S. verrucosum transcriptome with other species (S. portulacastrum, S. humifusum, and M. crystallinum) of Aizoaceae and Arabidopsis as an outgroup (Figure 2A). This analysis revealed 59,963 clusters, 59,543 ortholog clusters containing at least two species, and 420 single-copy gene clusters ( Figure S2). The five species shared 8367 orthogroups ( Figures S2 and S3).     A phylogenetic analysis was performed with OrthoFinder [18], and the results confirmed the relationship of the Aizoaceae species with its relative species of the family Cactaceae and Amaranthaceae ( Figures 1B and 2B). The FUSTr positive selection analysis identified 101 gene families under positive selection (Table S4). We examined the protein family (PFAM) enrichment of these genes using dcGO [19] and found that GO terms related to abiotic osmotic and salt stress response, such as response to abscisic acid, response to salt stress, and response to water, were uniquely enriched in the salt treat-ment groups ( Figure 2C). Interestingly, the GO term [32434] 'regulation of proteasomal ubiquitin-dependent protein catabolic process' was also enriched. The single-RING E3 ligase (SIAH1) is a representative gene linked to this GO term. We used the GO terms in Revigo analysis [20] to draw a closer relationship using hygrometric semantic plots of the different GO terms ( Figure 2D). Abiotic stress terms clustered closer to each other than other GO terms, including the regulation of proteasomal ubiquitin-dependent protein catabolic process.

Classification and Comparison of E3 Ligase Abundance
Expressed metabolic pathway enzymes were identified using KEGG-KAAS. A total of 11,484 (~9%) transcripts associated with various metabolic pathways were identified (Table S5) from the root transcriptome of S. verrucosum. From the KEGG annotation, 632 transcripts that were involved in Ubiquitin system (ubiquitin-activating enzyme (E1), ubiquitin-conjugating enzyme (E2), E3 ubiquitin-ligases, deubiquitinating enzymes, and Ubiquitins and ubiquitin-like proteins) were identified (Table S6). Components of the ubiquitin system described herein are based on the KEGG-KO IDs retrieved in the KEGG-KAAS analysis. These results showed that E3 ligase enzymes (433 transcripts) are predominant in the S. verrucosum root transcriptome compared to other classes of ubiquitin enzymes. We further classified E3 ligase enzymes into different sub-categories and found that single RING-finger E3 ligase enzymes (212 transcripts) and multi-subunit RING-finger E3 ligase enzymes (186 transcripts) were predominantly expressed. Other types of E3 ligases were also found, including atypical E3 ligases, HECT-type E3 ligases, U-box type E3 ligases, and UBL E3 ligases ( Figure 3A). Furthermore, we compared the number of single RING-finger E3 ligases with ten different plant species ( Figure 2B). The E3 ligase-based phylogenetic relationship between the ten selected plant species is shown in Figure S4. Overall, the single RING-type E3 ligase showed~37-50% abundance compared to other E3 ligase classes in all the species compared ( Figure 3B).

Differential Gene Expression Analysis
Differential gene expression analysis showed 3513 significantly upregulated genes (log2 fold change ≥ 1 and p < 0.05) and 1967 downregulated genes (log2 fold change ≤ −1 and p < 0.05) in samples from the 'early' salt treatment group. Comparison of the control with the 'late' salt treatment group showed 2367 upregulated genes and 2902 downregulated genes. The E3 ligase differential expression was limited to 13 genes at that magnitude change ( Figure 4A). This is highlighted in the UpsetR intersection plot of different transcript expressions at early and late stages, including E3 ligases, which is represented at the tail of the distribution ( Figure 4B). The majority are significantly (p-value < 0.05) represented by the single RING-type E3 ligase ( Figure 4A,B, Table S6).  To confirm the RNA-Seq-based differential gene expression results, we carried out qPCR experiments for 14 genes. We selected seven single-RING E3 ligase genes and seven highly/moderately expressed other genes for the RT-qPCR analysis (Table S7). The correlation analysis between the RNA-Seq results and qPCR results showed a good correlation (p-value < 0.001; r = 0.8 (Early) and r = 0.7 (Late)) between the studies ( Figure 4C). Some of the E3 ligase transcripts were upregulated in the early stage of salinity stress, which decreased significantly under prolonged salinity stress. Oxidative stress-responsive enzymes such as 2-alkenal reductase and peroxidase 1 were also upregulated (Table S7).

Yeast Functional Analysis of E3 Single-RING E3 Ligases
To determine if the differentially expressed single RING E3 ligase genes affect salinity tolerance, we selected 12 genes whose transcripts were upregulated upon salinity stress for a transgenic yeast system to explore their functionality. The yeast functional analysis was performed using a spot assay on medium containing 1.2 M NaCl. Yeast transformed with orthologs of BAH, SAP5, RHA, RNF141, RNF185, DRIP1, RHF2A1, RZFP, and SDIR1 show higher growth under salinity stress as compared to an empty vector control ( Figure 5A). This result is consistent with the image analysis of colony growth density ( Figure 5B). While many of the selected E3 ligases have not previously been linked to abiotic stress tolerance, several studies have shown that SAP5, SDIR, and DRIP play important roles in regulating the abiotic stress response. The yeast growth assay supports the potential involvement of these E3 ligases in mediating salinity stress response. Further study is needed to elucidate their function in plant stress response.

Discussion
Soil salinity is considered a major threat to plant growth and productivity. To overcome the inhibitory effects of salt, plants have evolved complex physiological and biochemical adaptations. Lately, Halophytes have gained much attention due to their remarkable potential to survive under salinity stress. Their ability not only to survive but also to complete their life cycle make them ideal models to unravel the molecular mechanisms of salinity stress adaptation. In this study, we generated the first high-depth RNA-Seq data for the halophyte S. verrucosum under salinity stress. The de novo assembly generated approximately 131,454 transcripts, which represent~99.3% single copy plant orthologous genes, indicating a high-quality transcriptome assembly compared to other halophytes already reported [21,22]. Approximately 53% of the genes were annotated using the UniProt-trEMBL database.
Based on the evidence of large-scale adaptive radiation in the Aizoaceae family, we investigated the rapidly evolving gene families in S. verrucosum. Our positive selection results identify gene families implicated in abiotic stress tolerance. These families include genes coding for E3 ligase. Our characterization of S. verrusosum E3 ligases revealed RING, UBOX, HECT, and RBR classes. Among the different classes of E3 ligases, the RING-type, especially the single RING-type E3 ligases, are reported to play an important role in tolerance to abiotic stresses, including salinity and drought [23]. The RING-type E3 ligases are the most abundant in S. verrucosum, and they are comparable to the number of RING-type reported in peach, grape, soybean, poplar, and maize [24]. Approximately 48% of the identified E3 ligases transcripts in S. verrucosum are the RING-type. In the Arabidopsis genome, nearly 36% of the E3 ligases are RING-type, which accounts for almost 2% of the predicted proteins [4,25]. Their abundance suggests important and diverse roles of the RING-type E3 ligases in plant growth and development.
In the present study, yeast-based functional characterization of differentially expressed E3 ligases provided evidence of single RING-type E3 ligases for salinity tolerance enhancement in S. verrucosum. Together with orthologs of known E3 ligases such as XERICO and SDIR, we tested several E3 ligases in our transgenic yeast growth assays to elucidate their role in salinity tolerance. Interestingly, yeast colonies transformed with novel BAH-type, RNF 141-like, and RHA-type E3 ligases showed an increased tolerance to salinity stress, indicating their possible role in stress tolerance ( Figure 5A,B).
Stress response mechanisms in plants are classified into two major groups: abscisic acid (ABA)-dependent and ABA-independent mechanisms. In ABA-dependent stress responses, the expression of stress-responsive genes is regulated by cis-acting ABA-responsive elements (ABRE) and multiple transcription factors [26]. Whereas in the ABA-independent mechanism, dehydration responsive element (DRE) and DREB2A transcription factors are deployed for stress tolerance. The expression level of DREB2A is regulated by two RING-type E3 ligases named DRIP1 and DRIP2 [27]. Studies suggest that DRIP E3 ligases act as negative regulators of DREB2A-mediated drought stress responses in Arabidopsis.
Overexpression of DRIP1 resulted in significant delays in the expression of DREB2Adependent drought-responsive genes. On the other hand, Stress-Associated Protein 5 (SAP5), a A20/AN1-type E3 ligase, functions as a positive regulator of stress tolerance by promoting the degradation of DRIP proteins [28,29]. Our functional analysis showed both S. verrucosum DRIP and SAP5 orthologs positively regulated salt tolerance in yeast ( Figure 5); whereas expression data showed the upregulation of a DRIP ortholog in the early stage of the salinity stress, but this expression declined after prolonged exposure to salt stress. Furthermore, the expression of a SAP5 ortholog (TRINITY_DN3310_c0_g1_i1) is significantly (p-value: 0.0060) induced throughout salinity stress in S. verrucosum, suggesting a possible interaction between DRIP and SAP5 orthologs in S. verrucosum.
XERICO, a RING-type E3 ligase, is reported to be involved in stress tolerance through an ABA-dependent pathway [30]. In Arabidopsis, overexpression of XERICO triggered the upregulation of AtNCED3, a key enzyme in ABA synthesis, eventually leading to elevated endogenous ABA levels and promoting increased tolerance to drought stress [30]. Recent studies have also shown the overexpression of AtXERICO in rice as well as the overexpression of ZmXERICO1 in maize conferred improved tolerance to drought and salinity stress [31,32]. The expression of a XERICO ortholog in S. verrucosum was highly induced in early stages of salinity stress, which implies its role as a 'first-responder' in S. verrucosum to salinity stress. Interestingly, the transcript levels of several genes in the ABA synthesis pathway, including genes encoding for NCED enzymes, were upregulated in our data. These findings suggest functional conservation of XERICO genes in monocots and dicots.
Another RING-type E3 ligase, SALT-AND DROUGHT-INDUCED RING FINGER 1 (SDIR1), was reported to play a key role in the ABA-mediated stress signaling pathway by controlling transcriptional regulation of ABI3 and ABI5 TFs [12]. Overexpression of the SDIR1 gene in Arabidopsis yielded increased tolerance to drought stress [12]. The SDIR1 gene also regulates the ABA-depended salt stress responses in Arabidopsis [13]. Similar findings were reported in rice by the overexpression of OsSDIR1 [14]. Recently, SDIR1-INTERACTING PROTEIN1 (SDIRIP1), a negative regulator of ABA signaling, was identified as the substrate for SDIR1 gene. We observed a significant upregulation of SDIR1 transcript (TRINITY_DN6910_c0_g1_i54, p-value: 0.0018) throughout the early and late phases of salinity stress ( Figure 4A, Table S8). The yeast spot assay test also confirms the possible role of the SDIR1 gene in the salinity tolerance of S. verrucosum ( Figure 5).
In plants, Seven In Absentia (SINA) is a group of E3 ligases involved in many biological processes, including autophagy and drought responses. In Arabidopsis, five SINA members were reported, SINAT1-5 (SINA of Arabidopsis thaliana 1-5). SINAT5 is the first of the SINA family to be characterized, and it is shown to negatively regulate auxin dependent lateral root development in Arabidopsis [33]. Whereas in rice, a SINAT5 ortholog, Oryza sativa Drought-induced SINA protein 1(OsDIS1) acts as a negative regulator of drought stress by altering transcriptional regulation of stress responsive genes [34]. Overexpression of the AtSINAT2 gene led to drought tolerance by increasing stomatal closure, resulting in water loss mitigation in Arabidopsis [35]. In S. verrucosum, a SINA E3 ligase transcript closely related to Arabidopsis SINAT2 was significantly upregulated under salinity stress (TRINITY_DN1735_c0_g1_i19) (p-value: 0.0098) ( Figure 4A). The expression level was elevated throughout the salinity stress experiment (Table S7). This indicates that SINA-type ligases may also regulate salinity tolerance of S. verrucosum.
In conclusion, this study provides valuable insight into the molecular programming of E3 ligases in halophytes. In addition to the above-mentioned known E3 ligase orthologs, we have also identified the differential expression of several E3 ligases previously unknown to respond to abiotic stress in our salt treatments. Many of these new salt-responsive E3 ligases increased salinity tolerance in our transgenic yeast growth assay and provide the basis for further study of salinity tolerance and translational research for crop improvement.

Sample Preparation, RNA Isolation, and Sequencing
Salinity stress experiments of Sesuvium verrucosum were performed in a hydroponic system. S. verrucosum cuttings containing three nodes from the apical shoots were grown in 1/10 strength Murashige and Skoog solution for five weeks for rooting. The rooted plants were treated with 300 mM NaCl dissolved in 1/10 strength Murashige and Skoog liquid medium. The root samples were collected after the stress treatment at different time points and were labeled as 'Early' and 'Late' samples. The 'Early' samples include the pooled samples collected at 4 h and 24 h of salinity treatment, whereas the 'Late' samples represent pooled samples collected on the 3rd and 7th day of stress. The control root samples were collected from cuttings grown in 1/10 strength Murashige and Skoog liquid medium without salt. Three biological replicates were included for each sample set. The collected samples were flash-frozen immediately in liquid nitrogen and stored at −80 • C.
Total RNA was extracted from the collected root samples using a modified CTAB method [36]. The samples were homogenized in liquid nitrogen using a mortar and pestle. The homogenized root samples (100 mg) were added to 5 mL of CTAB (Cetyl trimethylammonium bromide) buffer (2% w/v), 2% (w/v) PVP (Polyvinylpyrrolidone K30), 100 mM Tris-HCl, pH 8.0, 25 mM EDTA (Ethylenediaminetetraacetic acid), 2.0 M irrespective of their taxonomic representation among other taxa. Isoforms were identified in the input sequences, and TransDecoder (https://github.com/TransDecoder/TransDecoder, accessed on 15 May 2021) was used to extract the best open reading frame for each transcript. Homology analysis was performed with Diamond [50] using peptide sequences, followed by grouping gene families using SiLiX [51]. Finally, alignments were carried out using MAFFT [52]. At least 15 amino acid sequences representing a gene family were tested at a site-specific level using the branch-site model implemented in FUBAR [53]. Genes under positive selection were compared to the DEG list to discern whether the positive selection and regulatory evolution are acting on the same genes.

RT-qPCR Validation
Real-time quantitative PCR analysis was performed to validate RNA-Seq expression patterns. The extracted RNA samples were converted to cDNA using a QuantiTect Reverse Transcription kit (Qiagen, Hilden, German) following the manufacturer's instructions. Fourteen differentially expressed genes (DEGs) were selected for qPCR detection, and the PP2A gene was selected as an internal reference for the analysis. The qPCR analysis was performed using the SYBR Green PCR master mix (Applied Biosystems, Waltham, MA, USA). Three technical replicates were analyzed for each biological replicate. The relative expression pattern was calculated using the 2 −∆∆Ct method [54].

Functional Analysis in Yeast
Twelve differentially expressed E3 ligase genes were selected for functional characterization in yeast. The selected genes were cloned into the pGADT7 vector and transformed in Saccharomyces cerevisiae (strain Y187) using the Yeastmaker TM Yeast Transformation System 2 (Takara Bio Inc., Kusatsu, Japan) and grown at 30 • C for three days. The presence of specific genes in transformed yeast colonies were confirmed using colony PCR. Spot assays of S. cerevisiae harboring the genes were performed on synthetic minimal media without Leucine (SD-Leu) plates as control, and SD-Leu containing 1.2 M NaCl as test. Transformed yeast cultures were grown in liquid SD-Leu until the culture reached an OD 600 of 1, and then 10-fold serial dilutions up to 10 −4 were spotted on each media [55]. Image analyses were performed according to Petropavlovskiy et al. [56]. The colony growth density of triplicates was compared to controls using the software ImageJ2 [57].