Genome-wide identification of the rubber tree superoxide dismutase (SOD) gene family and analysis of its expression under abiotic stress

Introduction

Plants grow in constantly changing environments and are exposed to various abiotic stresses during growth, development, and production, such as extreme temperature, drought, waterlogging, salt, and heavy metal poisoning (Zhu, 2016). Under normal environmental conditions, the levels of reactive oxygen species (ROS) level are low and in dynamic equilibrium in plant cells, although under abiotic stresses, the dynamic equilibrium of ROS becomes disrupted and a large amount of ROS will be accumulated (Hasanuzzaman et al., 2020; Nadarajah, 2020). ROS include hydrogen peroxide (H₂O₂), superoxide anions (·O₂⁻), hydroperoxyl radicals (·HO₂⁻), hydroxyl radicals (·OH⁻), and alkoxy radical (Del Rio, 2015; Waszczak, Carmody & Kangasjarvi, 2018). Enrichment of ROS can cause protein oxidation, DNA damage, lipid peroxidation, and even cell death (Choudhury et al., 2017; Gill et al., 2015; Gill & Tuteja, 2010).

To adapt to adverse environmental conditions, plants have evolved complex and efficient antioxidant defences, including enzymatic and non-enzymatic systems (Rajput et al., 2021). The enzymatic system is a vital component and consists of many enzymes, such as SOD, CAT and POX (Gill & Tuteja, 2010; Pereira, 2016; Rajput et al., 2021). SOD acts as the first line of protection in plant cells and can catalyse superoxide anion (·O₂⁻) discordance to generate O₂ and H₂O₂, effectively clean up ROS and protect cells from harm, and is an antioxidant enzyme for scavenging ROS in plant cells (Rajput et al., 2021). SOD is a vital metalloenzyme and is categorized into three main subgroups according to different metal cofactors (Cu, Zn, Mn, and Fe), namely Cu/Zn-SOD, Fe-SOD, and Mn-SOD (Fink & Scandalios, 2002). SODs localized in various cellular organs. For instance, Cu/Zn-SOD subfamily members are localized in the cytoplasm, chloroplasts, peroxisomes, and extracellular regions. whereas Fe-SODs are localized in chloroplasts, and Mn-SODs are localized in the mitochondria and peroxisomes (Huseynova, Aliyeva & Aliyev, 2014; Rajput et al., 2021; Stephenie et al., 2020).

Plant SOD genes can be induced by abiotic stresses to cope with adverse environments (Choudhury et al., 2017; Gill et al., 2015; Gill & Tuteja, 2010; Rajput et al., 2021; Wang et al., 2016b; Waszczak, Carmody & Kangasjarvi, 2018). For example, in Arabidopsis thaliana, SOD genes increased significantly under oxidative stress (Kliebenstein, Monde & Last, 1998). Under abiotic stresses (salinity, cold, waterlogging, and drought), eight BnSOD genes are significantly up-regulated in rapeseed (Su et al., 2021). Most SOD genes exhibit transcriptional responses to drought stress in poplar tree (Molina-Rueda, Tsai & Kirby, 2013). In tea plants, most SOD genes are induced after exposure to cold and drought stress (Zhou et al., 2019). In addition, according to Shiraya et al. (2015), in rice, OsMSD1 overexpression enhanced tolerance to heat, and in contrast, MSD1-knockdown rice was notably susceptible to heat stress. TaSOD2 overexpression in both wheat and Arabidopsis increased SOD activity and enhanced resistance to salt and oxidative stress (Wang et al., 2016a). After heterologous expression of the Cu/Zn-SOD gene from Jatropha curcas, salt tolerance in transgenic A. thaliana was improved (Liu et al., 2015). Moreover, the joint overexpression MeCu/ZnSOD and MeCAT1 improved cold and drought tolerance in cassava (Xu et al., 2013). In a recent study, Han et al. (2019) over-expressed LkSOD2, LkSOD4, and LkSOD6 genes from the Japanese larch in A. thaliana and improved the resistance to salt stress. Taken together, these studies show that SOD genes play vital roles in the plant responses to abiotic stresses. Recently, with the rapid development of genome sequencing, the SOD gene family members of different species were comprehensively identified at the genome-wide level in many plant species, e.g., Brassica napus (Su et al., 2021), A. thaliana (Kliebenstein, Monde & Last, 1998), Brassica rapa (Iqbal Qureshi et al., 2021), Salvia miltiorrhiza (Han et al., 2020), Dendrobium catenatum (Huang et al., 2020), Camellia sinensis (Zhou et al., 2019), Zostera marina (Zang et al., 2020), Vitis vinifera (Hu et al., 2019), Gossypium hirsutum (Wang et al., 2017) and Solanum lycopersicum (Feng et al., 2016).

Natural rubber is a vital and renewable raw industrial material (Yamashita & Takahashi, 2020). To date, at least 2,500 species of plants can synthesise natural rubber (van Beilen & Poirier, 2007), but only the rubber tree (Hevea brasiliensis) is capable of producing high-quality commercial natural rubber and constitutes more than 98% of the global rubber production (Supriya & Priyadarshan, 2019). However, rubber trees are frequently subjected to numerous abiotic stresses that affect the plant growth and development, ultimately reducing the latex yield (Mohamed Sathik et al., 2018). To date, there are no reports regarding the identification of the SOD gene family in rubber tree. Therefore, in this study, a genome-wide comprehensive analysis of SOD genes in the rubber tree genome was performed and physicochemical properties and bioinformatics were analysed. Furthermore, tissue-specific expression and expression patterns under various abiotic stresses are comprehensively analysed. This study lays the foundation for further insight into the biological functions of the SOD gene family members and molecular breeding for stress resistance.

Materials and Methods

Results

Gene structures and conserved motifs

The study of gene exon/intron patterns assists understanding the evolution of the development of rubber tree SOD family genes. Therefore, we comprehensively analysed the gene structure and conserved motifs of HbSOD genes. The results displayed that the number of exons and introns of HbSOD genes ranged from four to eight and three to seven, respectively (Fig. 1B). The Cu/Zn-SOD group contained six exons and five introns, except HbCSD3, which contained seven introns and eight exons, and HbCSD4 which contained four exons and three introns. The Fe-SOD subfamily had six to seven exons and five to six introns. The HbMSD1 contained seven exons and six introns, and the number of exons/introns of the HbMSD2 gene contained three introns and four exons. In summary, the number of exons and introns of HbSOD genes was significantly different but had similar intron/exon patterns in the same evolutionary branch. To further investigate the gene structure and function of HbSOD genes, we calculated the motifs of nine HbSOD proteins using MEME software and visualized them using the TBtools software (Fig. 1C). The results showed that eight motifs were identified in HbSODs (Fig. 1D; Table S3). Motif 3 was widely distributed across almost all HbSOD genes, except for HbFSD1. Motif 5 was presented in most proteins of HbSOD gene family, except for HbFSD2 and HbMSD2. This indicates that Motif 3 and 5 were relatively conserved in the evolution of HbSOD proteins in the rubber tree. Motif 1 is a unique conserved domain of Cu/Zn-SOD subfamily members, and motif 6 only exists in HbCSD4 and HbCSD5. Members of the Fe-SOD and Mn-SOD subfamilies contained motifs 2 and 4.

Phylogenetic relationship analysis of HbSOD genes

A phylogenetic relationship of 32 SOD proteins from H. brasiliensis (nine), A. thaliana (eight), R. communis (seven), and P. euphratica (eight) (Table S2) was constructed using the NJ method of MEGA6.0 software. As we can see from the phylogenetic tree (Fig. 2), SODs were clustered into three major groups: Cu/Zn-SOD, Fe-SOD, and Mn-SOD. The Cu/Zn-SOD group contained 16 SOD members (three AtSODs, five HbSODs, five PeSODs and three RcSODs), Fe-SOD nine members (three AtSODs, two HbSODs, two PeSODs and two RcSODs), and Mn-SOD seven members (two AtSODs, two HbSODs, one PeSODs and two RcSODs). Further analysis revealed that the Cu/Zn-SOD subfamily was divided into four branches, where HbCSD4, HbCSD5, RcCSD3, and PeCSD5 were clustered together individually. HbSODs were more closely related to RcSODs and PeSODs in each subclade, particularly in the Fe-SOD and Mn-SOD subfamilies.

Prediction of cis-elements in promoters of HbSOD genes

The analysis of cis elements in the promoter region of HbSODs can assist exploration of the potential mechanisms of genes in response to multiple stresses. We predicted the cis elements of the 2,000 bp upstream sequences of the nine HbSOD genes, and the details of the cis-acting elements are shown in Table S4. One hundred and fifty-eight cis elements were identified from all HbSOD promoters, which were related to hormones and stress responses (Fig. 3; Table S4). Five types of cis-acting elements are related to plant hormone responses to abscisic acid (ABA), auxin, gibberellin (GA), methyl jasmonate (MeJA), and salicylic acid (SA) (Fig. 3). Further analysis showed that most HbSOD genes contained MeJA-responsive elements CGTCA-motif and TGACG-motif (except for HbCSD2 and HbMSD2) and ABA-responsive element ABRE (except for HbCSD2 and HbCSD5), GA-responsive elements TATC-box, P-box, and GARE-motif (except for HbCSD2, HbCSD3, and HbFSD1). In summary, MeJA-, ABA- and GA-responsive elements are largely distributed in most HbSOD genes, suggesting that HbSOD genes may be involved in growth and development as well as stress response regulation mediated by the above plant hormones. Moreover, five different stress-related response elements were identified, including drought, low temperature, stress resistance, defence responce, anaerobic induction, and light response elements. Furthermore, most HbCSD genes contain drought stress response elements (MBS), except for HbCSD2, HbFSD1, and HbMSD2. Low-temperature stress response element (LTR) was identified in HbCSD4, HbFSD1, and HbMSD2. HbCSD3, HbFSD1, and HbMSD2 contain defence and stress response elements (TCA-rich repeats). Nearly all HbSOD genes had an anaerobic induction response element, except for HbCSD1.

Functional annotation analysis of HbSOD genes

GO function annotation of genes allows us to better understand the molecular functions of proteins (Ashburner et al., 2000; Wimalanathan et al., 2018; Wimalanathan & Lawrence-Dill, 2021). Here, we uploaded the HbSOD protein sequences to the eggNOG website for annotation, and TBtools were used for GO enrichment analysis (Fig. 4; Table S5). The results of biological process annotation showed that all nine HbSODs were involved in cellular processes (GO:0009987) and biological processes (GO:0008150). Most HbSODs were involved in response to superoxide (GO:0000303), response to oxygen radical (GO:0000305), response to stress (GO:0006950), cellular oxidant detoxification (GO:0098869), response to an inorganic substance (GO:0010035), response to oxidative stress (GO:0006979), response to a toxic substance (GO:0009636), response to stimulus (GO:0050896), response to chemicals (GO:0042221), and response to abiotic stimuli (GO:0009628) (Fig.4; Table S5). Molecular function (MF) annotation showed that most HbSODs were related to copper ion binding (GO:0005507), superoxide dismutase activity (GO:0004784), zinc ion binding (GO:0008270), oxidoreductase activity (GO:0016491), antioxidant activity (GO:0016209), and metal ion binding (GO:0046872) (Fig. 4; Table S5). Results of cell component annotation showed that HbSODs mainly involved cellular components (GO:0005575), cytoplasm (GO:0005737), chloroplast (GO:0009507), and mitochondrion (GO:0005739) (Fig. 4; Table S5). The results of gene function annotation showed that HbSODs are enriched not only in the process of responding to various stress and scavenging superoxide radicals but also in cell composition, SOD activity, and metal ion complexes.

Expression patterns of HbSOD genes in various tissues

Expression pattern of HbSOD genes in different tissues will assist understanding the function of SOD gene family members in rubber trees. We investigated the expression patterns of nine HbSOD genes in budburst, copper-drown, light green mature leaves, roots, and bark using qRT-PCR. As shown in Fig. 5, nearly all HbSOD genes increased with the development of the leaf and showed higher expression in the mature leaves than in the other developmental stages leaves. We also found that HbSOD2 and HbSOD5 were highly expressed in the bark, whereas expression levels of the other HbSODs were low. Most HbSOD genes had low expression in the roots, whereas HbSOD5 was high.

Expression profiles of HbSOD genes under abiotic stress

Climate change and the associated adverse abiotic stress conditions (extreme temperature, drought, and salinity) influence rubber tree growth and development, ultimately decreasing rubber yield. To investigate the response of HbSOD genes to abiotic stress, we used qRT-PCR to analyse the expression patterns of HbSODs under low temperature (4 °C), salt (300 mmol/L NaCl), heat (42 °C), and drought stress. Different genes exhibit different expression patterns under various abiotic stresses conditions. Under low-temperature stress, the expression of HbCSD2 and HbCSD4 was significantly upregulated at almost all treatment time points, but not HbCSD2, at 1 h (Fig. 6; Table S6). HbCSD2 and HbCSD4 reached their maximum expression levels at 6 h after cold treatment and were 2.81 and 16.03 times that of the control, respectively. HbCSD1, HbCSD3, HbCSD5, HbFSD1, and HbMSD1 were significantly upregulated at 48, 48, 48, 3, and 12 h, respectively, and there were no significant changes at other time points. HbCSD5 was strongly induced by heat stress and maintained a high expression level, reaching a peak at 24 h (36.7 times than of the control) (Fig. 6; Table S6). The HbCSD4 expression at 1 and 3 h was 9.7 and 5.6 times that of the control, respectively, and there was no significant change in other time points. HbCSD2, HbCSD3, HbFSD2, and HbMSD1 were downregulated under heat stress (Fig. 6; Table S6). Under drought stress, nearly all HbSOD genes had relatively high expression levels, except for some genes (Fig. 6; Table S6). For example, HbCSD2, HbCSD3, HbCSD4, and HbCSD5 expression was significantly up regulated, except at 48 h time points. Among them, HbCSD2 and HbCSD4 reached a peak at 3 h and were 14.0 and 7.8 times that of the control, respectively (Fig. 6; Table S6), and HbCSD3 and HbCSD5 reached the top at 1 h and 5.8 and 4.7 times that of the control, respectively. Under salt stress, all HbSOD genes were significantly induced or suppressed. The expression of five HbSOD genes, including HbCSD2, HbCSD4, HbCSD5, HbFSD1, and HbMSD2 showed similar trends, a ‘first increased and then decreased’ trend (Fig. 6; Table S6). HbCSD1, HbCSD3, HbFSD2, and HbMSD1 were significantly downregulated and showed relatively low expression.

Discussion

The rubber tree is a vital tropical capital crop that produce rubber. Recently, with the steady growth in global natural rubber consumption, rubber tree planting has gradually developed into non-traditional rubber planting. Therefore, rubber trees often suffer from various abiotic stresses that affect their growth and development, ultimately decreasing rubber yield (Carr, 2011; Silpi et al., 2006). SOD genes are a vital gene family that encodes SOD proteins, which can scavenge ROS from cells and serve as the first line of defense in the antioxidant defense system (Fink & Scandalios, 2002; Rajput et al., 2021). In the current study, we identified nine SOD gene family members in the whole genome of the rubber tree, including five Cu/Zn-SODs, two Fe-SODs, and two Mn-SODs (Table 1). The number of SOD gene family members of the rubber tree was mostly similar to that in A. thaliana (8), S. miltiorrhiza (8), foxtail millet (8) (Wang et al., 2018), tomato (9) (Feng et al., 2016), cucumber (9) (Zhou et al., 2017), tea plants (10), and less than that in wheat (26) (Jiang et al., 2019), upland cotton (18), and rapeseed (31). Gene duplication events can occur in various ways, including fragment duplication, whole-genome duplication, and tandem duplication, which is one of the main molecular mechanisms of genome evolution to produce genes with new functions (Cannon et al., 2004; Moore & Purugganan, 2003). Therefore, differences in the number of SOD genes among plant species may be triggered by gene duplications, and the duplication events might play a vital role in SOD gene family expansion.

Numerous research results showed that the SOD genes of different species were clustered into three subfamilies (Fink & Scandalios, 2002; Rajput et al., 2021; Sheng et al., 2014). We analysed the evolutionary relationship of SOD proteins in H. brasiliensis, A. thaliana, R. communis, and P. euphratica, which were clustered into three subfamilies (Fig. 2): Cu/Zn-SOD, Fe-SOD, and Mn-SOD. Further analysis revealed that Cu/Zn-SOD proteins could be divided into four branches, with HbCSD4, HbCSD5, RcCSD3, and PeCSD5 clustered together (Fig. 2). Conserved domain analysis revealed that HbCSD4 and HbCSD5 contained a HMA domain adjacent to the SOD-Cu domain (Fig.1; Table 1). Similar results have been reported in other species, such as SlSOD4 and BjuACSD4/6-8 (Verma, Lakhanpal & Singh, 2019). The structure and number of introns and exons play vital roles in the evolution of gene families (Xu et al., 2012; Xu et al., 2019). Here, the number of exons and introns of HbSOD genes was significantly different; however, they had similar intron/exon patterns in the same evolutionary branch. It may be due to differences in SOD genes caused by exon/intron gain/loss, exonisation/pseudo-exonisation, or insertion/deletion in response to different environmental stresses during the evolution of rubber trees (Xu et al., 2012). Moreover, we found eight motifs in HbSODs (Fig.1; Table S3). Motif 1 is a unique conserved motif of Cu/Zn-SOD subfamily members. Motif 6 existed only in HbCSD4 and HbCSD5, and motifs 2 and 4 were distributed in the members of the Fe-SOD and Mn-SOD subfamilies. The HbSODs in each evolutionary branch had similar conserved motifs. In summary, similarities in conserved patterns and gene structures confirmed the inferred evolutionary relationship. Similar results have also been reported for crops such as corn (Liu et al., 2021) and rapeseed (Su et al., 2021).

Gene expression is regulated by cis elements and transcription factors in the promoter region. An analysis of the cis elements of the promoter of HbSOD genes proves helpful for further understanding the function of HbSOD genes in coping with abiotic stress. Here, the promoters of HbSOD genes were analysed using the online analysis software PlantCARE. At least one plant hormone response element was detected in most HbSOD genes (Fig. 3; Table S4). Therefore, HbSODs may be involved in the regulation of rubber tree growth, development, and stress response mediated by plant hormones. Nearly all HbSODs contain more than one cis element in response to abiotic stresses (drought, low temperature, defence, etc.). The prediction results of the cis-acting element suggested that HbSOD genes might play a vital role in the abiotic stress response of the rubber tree. These results were further confirmed by GO functional annotation analysis. Moreover, there have been similar reports on many different crops (Jiang et al., 2019; Su et al., 2021).

In the current study, most HbSOD genes increased in expression with leaf development and expressed were highly expressed in mature leaves (Fig. 5). Huang et al. (2020) reported that most DcaSOD genes were highly expressed in the leaves. Similarly, several genes displayed higher expression in poplars (Molina-Rueda, Tsai & Kirby, 2013). In soybean, GmCSD1, GmCSD2, GmCSD3, and GmCSD6 are highly expressed in the leaves (Lu et al., 2020). These results suggest that SODs may be involved in leaf growth and development. At present, the expression patterns of SOD genes under abiotic stress conditions have been widely reported. For instance, under low-temperature stress, The SmCSD2 of S. miltiorrhiza was upregulated during the entire stress treatment; under salt stress, almost all SmSODs upregulated their expression except SmFSD2; under drought stress, the expression patterns of SmCSD1, SmCSD2, SmCSD3, and SmMSD increased first and then decreased (Han et al., 2020). In tea plants, under low-temperature (4 °C) stress, the expression of nine CsSOD genes was significantly up-regulated; CsCSD2, CsCSD3, CsCSD6, and CsCSD7 were upregulated under drought stress, suggesting that CsSODs may be play a vital role in tea plants in response to cold stress and drought stress (Zhou et al., 2019). In soybeans, only GmFSD2 expression was significantly upregulated in 13 SOD genes under low-temperature stress, and GmFSD3, GmFSD5, and GmCSD5 were upregulated under salt stress (Lu et al., 2020). In watermelon and melon, most SOD genes were upregulated by cold, drought, and salt stress; but different members had different response modes under varying stresses (Zhang, Ding & Wei, 2021). These results suggest that SODs play vital roles in plants against the adverse environments and that homologous genes respond differently to stress in different plants. Here, the expression patterns of HbSOD genes under four abiotic stresses (cold, heat, drought, and salt) were tested using qRT-PCR. The results showed that nine HbSOD genes were differentially expressed under different stresses conditions. For example, under cold stress, most HbSODs were upregulated to varying degrees, especially HbCSD4, which was significantly upregulated at nearly all treatment time points (Fig. 6; Table S6). and the expression level was 16 times higher than that of the control after 6 h of treatment (Fig. 6; Table S6). Furthermore, the involvement of the low-temperature-responsive cis element LTR was detected in HbCSD4, indicating that HbCSD4 might play a vital role under cold stress. Under heat stress, all genes had different degrees of induced expression (Fig. 6; Table S6), in which HbCSD4 and HbCSD5 were strongly induced and up-regulated by heat stress. All genes, except HbMSD2 differentially increased their expression levels under drought stress (Fig. 6; Table S6), and HbCSD2, HbCSD3, HbCSD4, and HbFSD1 significantly up-regulated, except at 48 h time points. At least one involved drought stress-responsive cis elements including MBS, ABA-responsive element (ABRE), and defence and stress response were identified in all HbSOD genes, which could explain the significantly upregulated expression of HbCSD2, HbCSD3, HbCSD4, and HbFSD1 under drought stress. The HbCSD2, HbCSD4, HbCSD5, HbFSD1, and HbMSD2 expression was differentially upregulated and showed a ‘up-then-down’ expression trends under salt stress (Fig. 6; Table S6). The HbCSD2 and HbCSD4 levels were highly significant in response to multiple stresses. This suggests that HbCSD2 and HbCSD4 may be vital genes in the rubber tree in response to various abiotic stresses. Further studies are required to elucidate the full function of HbSODs.

Conclusions

Here, we identified nine SOD genes in the whole genome of the rubber tree that could be clustered into three subgroups: Cu/Zn-SOD, Fe-SOD, and Mn-SOD. To gain further insight, we predicted gene structures and conserved motifs, and found similar exon/intron patterns and conserved motifs within the same evolutionary branch. Cis elements and GO annotation analysis suggested that HbSOD genes may play a vital role in abiotic stress response of the rubber tree. Additionally, the expression of the majority of HbSOD genes increased gradually with leaf development and was highly expressed in mature leaves. HbCSD2 and HbCSD4 were significantly upregulated under low-, high-temperature and salt-stress conditions. Most genes expression levels were significantly upregulated under drought stress. The results of our study provide a basis for further studies on the biological functions of HbSOD genes in rubber tree growth, development, and responses to abiotic stresses.

Supplemental Information