Genome-Wide Evolution and Comparative Analysis of Superoxide Dismutase Gene Family in Cucurbitaceae and Expression Analysis of Lagenaria siceraria Under Multiple Abiotic Stresses

Superoxide dismutase (SOD) is an important enzyme that serves as the first line of defense in the plant antioxidant system and removes reactive oxygen species (ROS) under adverse conditions. The SOD protein family is widely distributed in the plant kingdom and plays a significant role in plant growth and development. However, the comprehensive analysis of the SOD gene family has not been conducted in Cucurbitaceae. Subsequently, 43 SOD genes were identified from Cucurbitaceae species [Citrullus lanatus (watermelon), Cucurbita pepo (zucchini), Cucumis sativus (cucumber), Lagenaria siceraria (bottle gourd), Cucumis melo (melon)]. According to evolutionary analysis, SOD genes were divided into eight subfamilies (I, II, III, IV, V, VI, VII, VIII). The gene structure analysis exhibited that the SOD gene family had comparatively preserved exon/intron assembly and motif as well. Phylogenetic and structural analysis revealed the functional divergence of Cucurbitaceae SOD gene family. Furthermore, microRNAs 6 miRNAs were predicted targeting 3 LsiSOD genes. Gene ontology annotation outcomes confirm the role of LsiSODs under different stress stimuli, cellular oxidant detoxification processes, metal ion binding activities, SOD activity, and different cellular components. Promoter regions of the SOD family revealed that most cis-elements were involved in plant development, stress response, and plant hormones. Evaluation of the gene expression showed that most SOD genes were expressed in different tissues (root, flower, fruit, stem, and leaf). Finally, the expression profiles of eight LsiSOD genes analyzed by qRT-PCR suggested that these genetic reserves responded to drought, saline, heat, and cold stress. These findings laid the foundation for further study of the role of the SOD gene family in Cucurbitaceae. Also, they provided the potential for its use in the genetic improvement of Cucurbitaceae.


INTRODUCTION
In natural conditions, plants are often susceptible to drought, salt, extreme temperatures, heavy metals, and other stresses that have a significant impact on the growth, development, and production of plants (Mittler and Blumwald, 2010;Cramer et al., 2011). When a plant is stressed, it adjusts its homeostatic apparatus by producing more reactive oxygen species (ROS) in its cells and ROS, toxic-free radicals produced under stress by plant cells that can oxidize the protein, destroy the cell membrane, and cause DNA damage (Lee et al., 2007;Karuppanapandian et al., 2011;Feng et al., 2015). These stresses inevitably accompany the development of ROS, including hydroxyl radicals (OH), superoxide anion radicals (O 2 − ), peroxide radicals (HOO − ), hydrogen peroxide (H 2 O 2 ), and singlet oxygen ( 1 O 2 ), which cause damage to the cell membrane, peroxidization and deterioration of macromolecules, and ultimately lead to the death of cells. ROS are also considered signaling molecules in different organisms and can affect various physiological processes in plants, For example, some prominent active oxygen scavengers can resist environmental stresses by regulating the expression of enzyme reaction family genes such as superoxide dismutase (SOD), catalase (CAT), peroxidase (POD), glutathione peroxidase (GPX), and peroxidase (PrxR) (Mittler et al., 2004;Ahmad et al., 2010;Bafana et al., 2011;Filiz and Tombuloğlu, 2015). Mainly, ROS are formed in the apoplast, mitochondria, plasma membrane, chloroplast, peroxisomes, endoplasmic reticulum, and cell walls (Mittler, 2017;Hasanuzzaman et al., 2020). Therefore, to manage ROS noxiousness, plants have established well-organized and composite antioxidant defense systems, including numerous nonenzymatic and enzymatic antioxidants.
SOD is a type of metal enzyme that is first found in bovine red blood cells and then characterized in bacteria, vertebrates, and higher plants (Mann and Keilin, 1938;Rabinowitch and Sklan, 1980;Tepperman and Dunsmuir, 1990;Kim et al., 1996;Zelko et al., 2002). Researchers of different studies have found that SODs can catalyze superoxide O 2 to disproportionate into O 2 and H 2 O 2 (McCord and Fridovich, 1969;Tepperman and Dunsmuir, 1990). SODs are detected in plants in roots, fruits, leaves, and seeds, which provide necessary protection for cells against oxidative stress (Tepperman and Dunsmuir, 1990). As per the binding pattern of metallic cofactors that cooperate with vigorous sites, SOD genes are categorized further into four groups: 1) iron FeSODs, 2) copper/zinc Cu/ZnSODs, 3) manganese MnSODs, and 4) nickel NiSODs (Mittler et al., 2004;Feng et al., 2016;Yan et al., 2016;Su et al., 2021). The different subtypes of SODs have almost comparable functions. However, they have different metallic cofactors and amino acid sequencing, in vitro subcellular location and crystal structure, and different hydrogen peroxide sensitivity (Abreu and Cabelli, 2010;Yan et al., 2016). These SODs are distributed in cell compartments individually and play an essential role in oxidative stress (Alscher et al., 2002). Among the SODs, Cu/ZnSODs are predominantly distributed in the chloroplast, extracellular space, and cytoplasm, and are present in certain bacteria and all eukaryotes, while MnSODs mostly present in plant mitochondria (Pilon et al., 2011;Feng et al., 2015;Wang et al., 2018;Hu et al., 2019;Verma et al., 2019;Lu et al., 2020;Silva et al., 2020). MnSODs in the plant genome play a role in eliminating ROS from mitochondria (Møller, 2001). FeSODs are primarily distributed in protozoa, prokaryotes, cytoplasms, and plant chloroplasts, while NiSODs are found in Streptomyces and also in some cyanobacteria, but not in plants (Youn et al., 1996;Wuerges et al., 2004;Miller, 2012).
Recent studies have shown that SODs can protect plants from abiotic stress such as heat, drought, cold, abscisic acid, salt, and ethylene (Wang et al., 2004;Pilon et al., 2011;Asensio et al., 2012;Feng et al., 2015). Several studies have shown that the SOD genes can be induced and transcribed in different plants under various stress conditions, such as heat, drought, cold, salt, osmotic stress, oxidative stress, and hormonal signal transduction (Wang et al., 2004;Pilon et al., 2011;Feng et al., 2015;Feng et al., 2016;Yan et al., 2016). SOD gene family under different hormones and abiotic stress conditions in rapeseed (Su et al., 2021), Salvia miltiorrhiza (Han et al., 2020), Zostera marina (Zang et al., 2020), and Hordeum vulgare  were recently published articles. Furthermore, diverse forms of SOD genes display different expression patterns under different stress conditions. In tomatoes, for example, SlSOD1 is the only significantly upregulated gene in the nine SlSOD genes, while SlSOD2, SlSOD5, SlSOD6, and SlSOD8 are regulated by salt stress. However, the expression levels of four "SlSOD2, SlSOD5, SlSOD6, and SlSOD8" genes are found high during imposed drought environment (Feng et al., 2016). Furthermore, the expression patterns of the same type of SOD gene were different under stress. For example, the studies found that there was no change in the expression of MnSODs in Arabidopsis under oxidative stress, and the researchers found that there was a significant change in the expression of MnSODs in peas, wheat, and Zostera marina under salt stress (Gómez et al., 1999;Wu et al., 1999;Baek and Skinner, 2003;Zang et al., 2020). These results suggest that different SOD genes have different expression patterns under different environmental stresses. Furthermore, researchers have also discovered that alternative splicing and miRNAs may be involved in the regulation of SOD expression (Srivastava et al., 2009;Lu et al., 2011). Until now, the SOD gene family has been described in many plant species, including Arabidopsis, Sorghum, Musa acuminata, Phaseolus vulgaris, and Populus (Kliebenstein et al., 1998;Zelko et al., 2002;Ballal and Manna, 2009;Fernández-Ocaña et al., 2011;Lin and Lai, 2013;Molina-Rueda et al., 2013;Filiz et al., 2014;Cui et al., 2015;Feng et al., 2015;Feng et al., 2016;Ho et al., 2017;Verma et al., 2019). Cucurbitaceae is a significant economic and nutritional crop in the world and SOD gene exists as a superfamily. The phylogeny line, genome circulation, gene assembly, preserved patterns, and expression profiles of these genes in different tissues have been thoroughly analyzed, laying the foundation for functional identification of SOD genes in Cucurbitaceae.

Retrieval of SOD Gene Family in Five Cucurbitaceae Species
To study the SOD gene family of cucurbit plants, the whole genome of individual species of Cucurbitaceae was searched by Frontiers in Genetics | www.frontiersin.org February 2022 | Volume 12 | Article 784878 the Blastp search method and Arabidopsis sequence of SOD used as the query (Kliebenstein et al., 1998 (Gasteiger et al., 2005) was used to evaluate the physiochemical properties (molecular weight, amino acid length, and isoelectric point) of each SOD protein.

Chromosomal Distribution of SOD Genes
The chromosomal localizations of Cucurbitaceae were obtained from CuGenDB database (http://cucurbitgenetics.org/). For chromosome mapping of SOD gene, map inspect tool (http:// mapinspect.software.com) was used. The aforementioned description is used to map the distribution of SOD gene throughout the Cucurbitaceae family individual species used in this study .

Intron/Exon Structure and Conserved Motif Analysis of Protein
The genome sequence and the coding sequence of the SOD gene were downloaded from the genome of the individual species of the Cucurbitaceae family from CuGenDB database. The Gene Structure and Display Server (GSDS) (http:/gsds.cbi.pku.edu.cn/) (Hu et al., 2015) was used to evaluate intron distribution dynamics and the splicing mechanism of each SOD gene. The Multiple Expectation Maximization for Motif Elicitation (MEME) software (Bailey et al., 2009) identified the conserved SOD protein motif (http:/meme-suite.org/tools/meme). Finally, the subcellular localization was predicted using WoLF PSORT (http://wolfpsort.org) (Horton et al., 2007).

Prediction of the SOD Cis-Regulatory Elements in the Promoter
For the prediction of regulatory elements on the promoter regions of CuSODs, 1,500 kb upstream DNA sequence of each SOD gene was collected from all SOD genes using the online web server PlantCARE (http:/bioinformatics.psb.ugent.be/webtools/ PlantCARE/html/) (Higo et al., 1999;Lescot et al., 2002).

Gene Ontology Annotation and MicroRNA Target Site Analysis
We used the CELLO v.2.5 software functional annotation platform to determine the Cucurbitaceae SOD gene's functional Classification. CELLO (http://cello.life.nctu.edu.tw/ site) platform connects the genes with GO terms through hierarchical vocabularies (Conesa and Götz, 2008). Functional enrichment analysis of SOD genes was performed using DAVID online tools (DAVID 6.8; http://david.ncifcrf.gov/) (Huang et al., 2009). The GO terms were classified into three categories: biological process (BP), cellular component (CC), and molecular function (MF). The upregulated SOD genes and downregulated SOD genes were entered separately and p <0.01 was considered to indicate a statistically significant difference.
To determine the miRNA-mediated posttranscriptional regulation of SODs, we searched the 5′ and 3′ untranslated regions (UTRs), and the coding regions, of the SODs for target sites of Lagenaria siceraria miRNAs obtained from various databases and published articles on the psRNA Target server using default parameters (Dai et al., 2018).

Plant Materials and Abiotic Stresses
The L. siceraria plants were planted in the growth chamber; seeds (variety: Winall 808) were provided by The National  . After treatment, the leaves were collected at 0, 4, 8, 16, and 24 h, and instantly frozen in liquid nitrogen immediately, and the sample was obtained to determine the transcription level of the SOD gene family in treated plant seedlings. Finally, all samples were instantly frozen in liquid nitrogen and stored at −80°C until further use.

RNA Isolation and Quantitative Real-Time PCR Reaction
RNA extraction reagent kit (Trizol) was purchased from Hua and Maike Biotechnology Co. Ltd., Beijing, China. Total RNA was extracted following the manufacturer's instructions. The RNA reverse transcription kit (TaKaRa Company, Japan) and Fluorescent Quantitative Reagent SYBR Green Master (Roche, United States) were used, and qRT-PCR was performed using a previously described method (Zhou et al., 2014). Gene-specific primers of LsiSODs and Lsi-actin for the qRT-PCR system were designed by using GenScript server (www.genscript.com) and synthesized by Sangon Biotech (Shanghai) (Supplementary Table S1). Actin gene of L. siceraria was used as an endogenous control to detect the relative expression of LsiSOD genes based on previous studies . Primers used are given in Supplementary Table S1. qRT-PCR reactions were performed in biological triplicates. The relative expression level was calculated by 2 −ΔCt method and statistical analyses were measured using Microsoft Office 2010.

Genome-Wide Characterization of SOD Genes in Five Cucurbitaceae Species
In total, 43 SODs were identified in five Cucurbitaceae species from CuGenDB, and 8-10 genes were retrieved for individual species such as Citrullus lanatus (8 SODs), Cucurbita pepo (10 SODs), Cucumis sativus (9 SODs), Lagenaria siceraria (8 SODs), and Cucumis melo (8 SODs) (Supplementary Table  S2). The detailed information of Cucurbitaceae species (molecular weight, starting and ending point, and theoretical pI) of the protein are listed in Table 1, and all protein sequences are shown in Supplementary Table S10.
In Cucurbita pepo, ten SOD genes were detected, including seven Cu/ZnSODs and three Fe/Mn-SODs; the molecular weight, length, and pI values of SOD proteins were within the ranges of 15.28-59.17 kDa, 52-530 amino acids, and 5.18-9.73, respectively. The subcellular localization prediction showed that Cu/ZnSODs CpSOD1, CpSOD3, CpSOD4, CpSOD6, and CpSOD7 were localized to the chloroplast, and two Cu/ ZnSODs, CpSOD2 and CpSOD9, were localized to the cytoplasm. Furthermore, CpSOD5 and CpSOD10 were localized to the mitochondrion.
In Cucumis sativus, nine SODs were retrieved-six were Cu/ ZnSODs and three were Fe/MnSODs. The molecular weight, length, and pI values of SOD proteins were within the range of 15.26-86.82, while amino acids range from 73 to 318 and 4.97-9.83 kDa, respectively. CsaSOD2, CsaSOD3, CsaSOD4, CsaSOD6, CsaSOD7, CsaSOD8, and CsaSOD9 were acidic and CsaSOD1 and CsaSOD5 were basic. In the prediction of subcellular localization, CsaSOD1 and CsaSOD2 were localized to the mitochondrion. CsaSOD3, CsaSOD4, CsaSOD6, and CsaSOD8 are localized to the chloroplast; CsaSOD5 and CsaSOD9 are in the nucleus; and CsaSOD7 is localized to the cytoplasm.
In Cucumis melo, a total of eight SODs were identified, including four Cu/ZnSODs and four Fe/MnSODs, and the SOD protein's molecular weight, length, and pI values were observed within the ranges of 15.90-36.95 kDa, 134-347 amino acids, and 5.16-8.49, respectively; three SODs such as MelSOD1, MelSOD4, and MelSOD8 are basic, and MelSOD2, MelSOD3, MelSOD5, MelSOD6, and MelSOD7 are acidic. The prediction of subcellular localizations showed that MelSOD1 and MelSOD8 are in the mitochondrion, and MelSOD3 and MelSOD7 are in the cytoplasm. MelSOD2, MelSOD4, and MelSOD5 are in the chloroplast ( Table 1).

Structural Phylogenetic Analysis of SOD Genes Family in Five Cucurbitaceae Species
To define the evolutionary relationship of SODs in different plant species, we constructed a phylogenetic tree of SODs based on the full length of protein sequences and divided them into eight groups. A total of 106 SOD proteins were obtained from 13 plant Frontiers in Genetics | www.frontiersin.org February 2022 | Volume 12 | Article 784878  Figure S2). Based on the phylogenetic tree, the SOD genes were divided into five subgroups according to Arabidopsis (Silva et al., 2020). However, three unique clades in a phylogenetic tree may have independent evolutionary trajectories from other clades. It is also suggested that these clades may have individual evolution that is different from Arabidopsis (Nijhawan et al., 2008;Wei et al., 2012;Manzoor et al., 2021a). According to the bootstrap values quoted on the nodes, topography, and sequence similarities, all identified SODs from the Cucurbitaceae species were categorized into eight subfamilies (I-VIII). SOD protein sequences from all Cucurbitaceae species contributed in all subfamilies ( Figure 1). These results suggested that there may have been gene loss or gain events that occurred throughout the evolutionary process. The gain and loss of specific SOD gene members caused functional divergence.

Conserved Motif Analysis of SOD Gene Family in Five Cucurbitaceae Species
The conserved motif analysis of the SOD family supported the classification and evolutionary relationships of the five Cucurbitaceae species SOD genes. In total, 20 motifs were detected from 43 SOD genes in Citrullus lanatus, Cucurbita pepo, Cucumis sativus, Lagenaria siceraria, and Cucumis melo.
All SOD genes contained at least two motifs except CpSOD3 and CsaSOD5, which contain only one motif ( Figure 2B): Cu/ZnSOD and Fe/MnSODs. Besides, FeSODs and MnSODs were clustered in the same group and belonged to a similar subcluster, while the Cu/ZnSOD was clustered in a different group. A similar cluster distribution was detected in the SOD proteins of each species, indicating that these SOD genes are highly conserved in different plants ( Figure 2A). In Cucurbitaceae species, 5-9 exons were detected in SOD genes by comparing the genomic DNA and CDS sequences using the GSDS 2.0 utility ( Figure 2C). However, differences were observed in the size and number of exons/ introns in Cucurbitaceae. In the organization of exons/introns, a high degree of conservation has been observed, which is consistent with the high degree of similarity found by multiple alignments between protein sequences, which gives high similarities between them. The evolutionary analysis suggested that structural gene diversity is the primary source for the evolution of multigene families (Xu et al., 2012;Manzoor et al., 2021b). C. lanatus contains 6-8 exons and C. pepo contains 3-9 exons; interestingly, CpSOD3 contains 2 exons with 1 intron. C. sativus contains 5-8 exons while CsaSOD9 has 3 exons and 2 introns. Both CpSOD3 and CsaSOD9 genes are the smallest genes among all 43 SOD genes identified in five Cucurbitaceae species. C. melo consists of 5-8 exons in which MelSOD8 has 5 exons and 4 introns. L. siceraria contains 6-10 exons; interestingly, LsiSOD1 and LsiSOD2 comprise 8/7 exons/introns. On the other hand, LsiSOD3, LsiSOD5, and LsiSOD6 comprise 6/5 exons/introns that have the upstream and downstream sequence. LsiSOD4 contains exons/introns (7/6) with upstream and downstream sequences and LsiSOD7 contains 10 exons with the only downstream sequence. The remaining LsiSOD8 contains 7 exons and 6 introns with both upstream and downstream sequences.

Chromosomal Distribution and Promoter Analysis of SOD Gene Family
The chromosomadal distribution of the SOD gene family of the Cucurbitaceae species was determined, as shown in Figure 3 (Supplementary Table S3), and detailed data are given in Supplementary Table S4. C. lanatus chromosome 2 contains two genes, and chromosome 3 contains three genes, while chromosomes 4, 7, and 10 contain only one gene. In C. pepo, chromosome 0 contains three genes; 1, 2, 5, 11, and 20 contain only one gene while chromosome 8 contains 2 genes. In C. sativus, chromosome 1 contains three genes, and chromosomes 4 and 6 contain two genes each, while only one gene was found on chromosomes 2 and 3. In L. siceraria, chromosomes 1, 2, 10, and 11 contain one gene, while chromosomes 6 and 7 contain two genes, respectively. C. melo contains two genes on chromosome 2, while chromosomes 5, 6, 7, 8, 11, and 12 contain only one SOD gene.
To clearly understand the function and regulation of SOD proteins, cis-elements in the promoter sequence of SOD genes in Cucurbitaceae were retrieved. The 1,500-bp upstream region of the start codon of each SOD gene was analyzed by using the PlantCARE database. According to the obtained results, the ciselements can be divided into three classes: stress related, light related, and hormone response elements ( Figure 4A). Five ciselements of stress response were determined, including LTR, TC rich, MBS, ARE, and box-w1; these elements reflected the response of plants to drought, low temperature, anaerobic induction, stress defense, and fungal inducer. Four hormonesensitive cis-elements (SA, MeJA, GA, and ethylene) were identified, which are associated with ABA, SA, ethylene, and MeJA responses ( Figure 4B). On the promoter region of SOD gene family of Cucurbitaceae, a considerable number of phytoresponse cis-elements were detected (Supplementary Table S4). The results suggested that the cis-elements of SOD promoter had a positive response to abiotic stress and regulation mechanism of plant growth and development.

Gene Duplication Events and Collinearity Relationships in Five Cucurbitaceae Genomes
To further understand the evolution and new function of genes, gene duplication events of the SOD gene family were identified in five Cucurbitaceae species (Citrullus lanatus, Cucurbita pepo, Cucumis sativus, Lagenaria siceraria, Cucumis melo). We analyzed three modes of gene duplications in all SOD gene  Table S5). Our study showed that duplication events play an important role in SOD gene expansion, and TRDs and WGD occurred at high frequency in Cucurbitaceae species ( Figure 6).
We further analyzed the collinearity relationships of SOD genes between five Cucurbitaceae species (Citrullus lanatus, Cucurbita pepo, Cucumis sativus, Lagenaria siceraria, Cucumis melo) as these five plants belong to the Cucurbitaceae family and shared a similar antique (Supplementary Table S6) (Figure 7). A total of 35 collinear gene pairs were found between the five Cucurbitaceae genomes, including 5 orthologous gene pairs between Arabidopsis and Cucurbita pepo, 4 orthologous gene pairs between Arabidopsis and Lagenaria siceraria, 4 orthologous gene pairs between Arabidopsis and Cucumis melo, 8 orthologous gene pairs between Citrullus lanatus and Cucumis sativus, 7 orthologous gene pairs between Citrullus lanatus and Lagenaria siceraria, and 7 orthologous gene pairs between Cucurbita pepo and Lagenaria siceraria, suggesting a close relationship among the five Cucurbitaceae genomes. The results show that the genetic relationship between SOD gene pairs in C. lanatus, C. pepo, C. sativus, C. melo, and L. siceraria is close. No pairs of collinear SODs are shared between Arabidopsis and five Cucurbitaceae genomes, indicating the long distance phylogenetic relationship between two species.

Gene Ontology Annotation Study of SOD Gene Family
Functional enrichment analysis of SOD genes was performed using DAVID. The SOD genes were categorized into three functional groups-biological processes, molecular functions, and cellular components-that are characteristics of genes or gene products, which enable us to understand the diverse molecular functions of proteins (Altschul et al., 1990). These results may be related to protein sequence similarities caused by genomic events (8, Supplementary Table S7). Evaluation of the biological processes mediated by SODs indicated that the same percentage (~27%) of proteins was involved in oxidoreductase activity and ion binding. Among the SOD proteins,~14.58% of members showed potential involvement in protein binding, enzyme binding, and DNA binding, respectively, while nucleic acid binding transcription factor activity involvement is~2.08% during the 8 LsiSOD genes of Lagenaria siceraria life cycle (Figure 8). At the same time, in SOD genes of five Cucurbitaceae species, "Cellular Components" includes 17.40% response to stress and reproduction, respectively. Meanwhile, 15.83% of genes were involved in aging, and 10.73% of genes were involved in homeostatic process and transport, respectively ( Figure 8). Furthermore, the SOD genes of 8 LsiSOD genes of Lagenaria siceraria were involved in the molecular functions, 11.9% of the genes were involved in intracellular, cell, cytoplasm, organelle, and mitochondrion while 10.51% of the genes were involved in plastid and 5.65% were involved in extracellular region, cytosol, nucleus, and extracellular space, respectively (Figure 8). The results corroborated the putative SOD promoter analysis (Jagadeeswaran et al., 2009;Jia et al., 2009;Guan et al., 2013;Manzoor et al., 2021a). Furthermore, analysis of the molecular function annotations revealed that all of the LsiSOD protein functions were enriched in SOD activity and metal ion binding.  Table S8). The analysis indicated that 3 LsiSOD genes (LsiSOD3, LsiSOD6, and LsiSOD7) were targeted on different miRNA. The expression profiles of these miRNAs and their targets are needed to detect and verify in further experiments to determine their biological functions in bottle gourd. The regulation of SOD genes' expression by miRNA needs to be studied further.

Tissue-Specific Expression of LsiSOD Genes
A strong link between gene expression and function has been suggested and that the SOD gene family is primarily involved in plant growth, development, and stress responses. To determine the biological functions in Lagenaria siceraria, the expression profiles of the 8 LsiSOD genes were analyzed in tissues (root, flower, fruit, stem, and leaf) using RNA-seq data downloaded from CuGenDB (http://cucurbitgenetics.org/) (Zheng et al., 2019). These results indicated that LsiSOD genes have tissuespecific expression patterns. However, LsiSOD4, LsiSOD5, and LsiSOD6 exhibited low transcript abundance in root and flower while their highest expression was detected in stems and leaves. This suggested that they play an important role in the development of plants. In addition, LsiSOD1, LsiSOD2, and LsiSOD3 are highly expressed in all parts of plants especially in the early stages of plant development ( Figure 10). All LsiSOD genes were widely expressed and showed tissue-specific expression patterns while LsiSOD8 shows less expression in all different developmental parts.

Gene Expression Analysis of LsiSOD Gene Under Abiotic Stresses
To understand the role of the SODs, qRT-PCR was used to analyze the expression patterns of the SOD gene in response to stress from heat, cold, drought, and NaCl. Significant differences were observed in the expression levels of the LsiSOD genes under different treatments, and complexity was detected in their expression patterns ( Figure 11). The 8 LsiSOD genes were induced by heat treatment, and their expression profiles were significantly enhanced. Among them, at 4 h of treatment one gene LsiSOD7 (2.62-fold), at 8 h of treatment two genes LsiSOD3 (20.35-fold) and LsiSOD4 (14.62-fold), at 16 h of treatment two genes LsiSOD2 (13.22-fold) and LsiSOD3 (16.41-fold), and at 24 h 2 genes LsiSOD2 (18.24-fold) and LsiSOD3 (24.16-fold) reached the highest expression level respectively. Interestingly, the transcription of LsiSOD3 and LsiSOD4 was gradually induced to 8 h, then decreased at 16 h, and finally reached a high level at 24 h, as described in Figure 11. Besides, in the expression of LsiSOD7, a dramatic increase was observed, which reached the maximum level at 4 h and decreased significantly at 8 h, then increased gradually with the treatment ( Figure 11A). During cold treatment at 4 h, many LsiSOD genes were downregulated, LsiSOD3 (1.17-fold) remained unchanged, LsiSOD7 (1.27-fold) slightly increased, and then LsiSOD3 (7.27fold) increased significantly at 8 h ( Figure 11B). After 8 h, except for the increased expression of LsiSOD2 (1.71-fold) at 16 h, the other LsiSOD genes were significantly downregulated at 16 h and continued to decline at 24 h as compared with 0 h. After 8 h, the expression of LsiSOD7 (4.53-fold) and LsiSOD8 (4.39-fold) was progressively induced, and the induction peak appeared at 16 and 24 h, respectively. Under the PEG stress treatment, almost all LsiSOD gene expressions increased significantly at 4 h except LsiSOD1, and then the level of transcription decreased. At the same time, LsiSOD1 decreased somewhat at 4 h, and its transcription remained unchanged ( Figure 11C). At 4 h of treatment with PEG, the expression of LsiSOD8 (4.27-fold) was maximum, whereas the expression of LsiSOD8 (1.36-fold) was highest during 8 h of treatment, while at 12 h of treatment LsiSOD4 (2.17-fold) shows maximum expression. During NaCl stress treatment, mostly the expression pattern of all LsiSOD genes increased tremendously at 8 h, decreased at 16 h, and finally increased at 24 h ( Figure 11D). At 4 h of treatment, expression of LsiSOD3 increased (3.71-fold), and at 8 h of treatment with NaCl, expression of LsiSOD2 increased (8.59-fold), while at 12 h of treatment, the expression of LsiSOD8 increased (3.01-fold) and the expression of LsiSOD4 increased (4.53-fold). It should be worth mentioning that LsiSOD2 transcription level increased significantly in 24 h, which was significantly higher than other LsiSOD genes.
SODs are identified as a crucial enzyme involved in many oxidation processes and protect the plant against ROS (Alscher et al., 2002). Therefore, we systematically analyzed the SOD gene family of Cucurbitaceae and determined the gene expression patterns in L. siceraria under different abiotic stresses (drought, heat, salt, and cold).
A phylogenetic study was conducted on SOD proteins in Lagenaria siceraria and four other cucurbit plants, namely, Citrullus lanatus, Cucurbita pepo, Cucumis sativus, and Cucumis melo. SODs could be classified into eight groups with  good statistical support in line with previous studies (Alscher et al., 2002;Feng et al., 2015;Filiz and Tombuloğlu, 2015;Feng et al., 2016;Wang et al., 2017;Verma et al., 2019). SOD genes are localized to chloroplasts, and mitochondria were classified into the same subgroups. Furthermore, nearly all LsiSODs are systematically grouped with at least one member of all plant species under consideration, highlighting the functional conservation of LsiSODs with other plant species SODs (Figure 1). In this study, SOD gene families were found in Cucurbitaceae genome, including Citrullus lanatus, Cucurbita pepo, Cucumis sativus, Lagenaria siceraria, and Cucumis melo. There are considerable differences in the size of genomes and the number of SODs in these plant organisms, but there is no substantial dependence on the size of genomes. The number of amino acids of SOD proteins ranged from 52 to 530 kDa, showing a significant variation. The difference in the number of SOD genes in different plant species may be due to gene duplication. In addition, the difference between clades might be related to different functions and diversity of exons/introns and conserved motif structure. Intron and exon variations play a Frontiers in Genetics | www.frontiersin.org February 2022 | Volume 12 | Article 784878 13 major role in the evolution of different genes. The variation in introns/exons and motif structure of SOD genes suggests a high level of complexity between Cucurbitaceae species.
We further investigated the diversity between SODs by analyzing their subcellular localization as localization plays a vital role in their functions. The cytoplasm (34%), chloroplast (33%), mitochondria (22%), and nucleus (11%) SOD gene localizations were predicted in five Cucurbitaceae species (Supplementary Table S8), and further study is needed to confirm those localizations. Increasing evidence points to an important role of miRNAs in stress tolerance. Several studies have reported that miRNAs regulate the expression of stressresponsive protein-coding genes at post-transcriptional level, showing a reverse correlation between the miRNA and the target expression (Sunkar et al., 2007;Lomakina, 2018;Manzoor et al., 2020). These miRNAs resulted from computational predictions and deep sequencing, and they are involved in some biological processes reported in plants, including responses to environmental stresses and regulating cell growth, development, and metabolism (Qiu et al., 2007;Li FIGURE 8 | Gene Ontology (GO) annotation results of SOD genes of five species in the Cucurbitaceae family. GO analysis of SOD protein sequences is analyzed for their involvement or function in three important categories: biological process, molecular function, and cellular component.  Table  S8). In Lagenaria siceraria, the miR164 family comprises 3 members that generate four mature products, miR164a/b/c and miR164f, which could target at least three LsiSOD genes. The miR164 family is a highly conserved miRNA that has been found in many plant species. miRNA164a thereby increased the tolerance of the plant to the abiotic stress or increasing the biomass, vigor, or yield of the plant. The miRNA miR164 plays a central role during the development of serrated leaf margins in Arabidopsis. In this study, transcripts of three miRNA families were identified because miRNA expression is highly regulated under different abiotic stress conditions ( Figure 9). Further studies are needed to determine the role of SODs in Lagenaria siceraria.
Gene duplications are an essential mechanism for creating genetic novelty in all plants, which could help organisms to adapt to environmental change (Alvarez-Buylla et al., 2000;Kondrashov, 2012;Manzoor et al., 2020). There are different kinds of gene duplication, including dispersed duplication (DSD), transposed duplication (TRD), and whole-genome duplication (WGD), which differentially contribute to the expansion of plantspecific genes (Qiao et al., 2018, Qiao et al., 2019. Gene replication can cause variation in the number of SOD genes in different plant species, and involves tandem and segmented replication, which plays an important role in SOD gene diversity, expansion, and duplication (Filiz and Tombuloğlu, 2015;Wang et al., 2016;Zhang et al., 2016;Wang et al., 2017). Thus, tandem repetition is likely to play a vital function in the amplification of the LsiSODs, for example, the two neighboring genes LsiSOD3 and LsiSOD6 on Chr6 (Figure 3). Characterization of the gene structure describes that the number of introns from the SOD gene family differed between the observed Cucurbitaceae species (Figure 2); in addition, the number of introns in L. siceraria was 6-8. Previous studies have shown that the plant's SOD gene has a strongly conserved intron pattern and that seven introns are present in most cytoplasmic and chloroplast SODs (Fink and Scandalios, 2002). In our study, only the LsiSOD members were predicted to contain seven introns, as shown in Figure 2.
Analysis of the cis-elements in SOD gene promoters resulted in the detection of three major types of cis-elements associated with light, abiotic stress, and hormone response as well as ciselements related to developmental processes and tissue-specific expression. A relatively large number of light-responsive cis-elements were detected in SOD gene promoters, suggesting that SODs might participate in the abiotic response. Some studies have shown that SOD genes are involved in the response to abiotic stress in different plants like Dendrobium catenatum, Pennisetum glaucum, maize, and Arabidopsis (Divya et al., 2019;Huang et al., 2020;Kliebenstein et al., 2020;Liu et al., 2021). In addition, a series of cis-elements related to abiotic stress responses were identified in SOD gene promoters, such as MBS, ERE, TC-rich repeats, ARE, ABRE, and Box-4, which may regulate gene expression under various stresses. Most of the SOD genes in Arabidopsis, banana, rice, tomato, poplar, cotton, and other different plants can be induced in response to various abiotic stresses such as heat, cold, drought, and salinity (Kliebenstein et al., 1998;Lee et al., 2007;Dehury et al., 2013;Nath et al., 2014;Feng et al., 2016;Wang et al., 2016;Zhang et al., 2016;Wang et al., 2017;Verma et al., 2019). GO annotation analysis further verified these results (Supplementary Table S7). In addition, some researchers have reported similar findings in different crops and found that SOD gene plays an important role under different stress conditions. GO annotation results confirmed the LsiSOD's role in response to different stress stimuli, cellular oxidant detoxification processes, metal ion binding activities, SOD activity, and different cellular components. These results can promote our understanding of LsiSOD genes under different environmental conditions. We also found protein similarities (Supplementary Table  S9) by using protein-protein interactions (PPIs) handle a wide range of biological processes, including metabolic and developmental control and cell-to-cell interactions (Mittler, 2017). Regarding the possible role of SODs from Cucurbitaceae, the expression profile of SODs in different tissues was analyzed based on sequence evidence from RNAseq data. Data analysis suggests that the gene family of SOD was expressed in all tissues; some SODs had tissue-specific expression patterns ( Figure 10).   (Gill and Tuteja, 2010). However, LsiSODs' specific response to heat, cold, drought, and salt is not well understood. Therefore, the analyses of qRT-PCR provide important clues for understanding the possible role of LsiSODs under various stresses. Based on the evaluation of ciselements of LsiSOD gene promoter, three main cis-element types related to light, abiotic stress, and hormone response, as well as cis-elements related to development process, were determined. The cis-elements in LsiSODs promoter include TC-rich motif, LTR motif, MBS motif, ARE motif, and ABRE motif, which is the evidence of abiotic stress response (Supplementary Table S4). These motifs were previously observed in different plant species to cope with abiotic stresses, such as bananas, tomatoes, Brassica, tobacco, and millet (Oppenheim et al., 1996;Kurepa et al., 1997;Kliebenstein et al., 1998;Pilon et al., 2011;Feng et al., 2015;Wang et al., 2018;Hu et al., 2019;Verma et al., 2019). In this study, the expression level of eight LsiSOD genes also changed significantly under different stresses, indicating that these genes play an important regulatory role in stress response and may have a certain functional relationship. Almost all LsiSOD genes were upregulated during heat treatment, and some displays have similar expression ( Figure 11C).
Under cold stress, the expression level of all LsiSOD genes is unregulated significantly, and their patterns of expression are different from each other, as shown in Figure 11D, which means that LsiSOD genes under cold stress may have functional diversity. Similar expression patterns were observed in LsiSOD1 with MaCSD1B and MaCSD2B (Feng et al., 2015). In PEG treatment, almost all Lsi-SOD1 genes have identical expression patterns, reach the highest level at 4 h, and then decline, as shown in Figure 11A, meaning that the function of the LsiSOD gene is related to drought.
The expression of most of the LsiSOD genes changed during NaCl treatment; LsiSOD2 is the only one of the eight LsiSOD genes with a significant increase, as shown in Figure 11B, indicating that LsiSOD2 plays an active role in detoxification of ROS during salt stress. Similar results were observed during salt stress in tomato (Feng et al., 2016). SlSOD1 and LsiSOD2 were grouped in the same subgroup, showing approximately high amino acid sequence homology. Altogether, we conclude that the LsiSOD gene plays a specific role in ROS removal induced by abiotic stress, which enhances plant adaptability to stress. In addition, some LsiSOD genes were correlated to abiotic stress exhibiting different expression patterns. For example, cold treatment reduced the expression of LsiSOD1, heat treatment and NaCl treatment increased the expression of LsiSOD1, while under PEG stress, the expression level remained the same without significant change, as shown in Figure 11, which means that the function of LsiSOD1 in different signaling pathways was different. It needs to be further explored to elucidate the role of the LsiSOD gene in L. siceraria under various abiotic stresses.

CONCLUSION
In this study, we identified SOD genes from the Cucurbitaceae family and analyzed their genomic structure, GO annotation analysis (molecular processes, biological functions, and cellular components), miRNA, gene duplication events (TD, PD, DSD, WGD, and TD), conserved motif patterns, phylogenetic relationships, mode of gene duplications, subcellular localization, and RNA-seq data analysis. qRT-PCR was used to evaluate the LsiSOD gene regulatory response under a variety of abiotic stresses, such as heat, cold, PEG, and NaCl. This research provided insight and further functional identification of the Cucurbitaceae family of SOD genes and laid the framework for understanding the molecular mechanism of the SOD gene in response to stress and plant growth. Genome-wide study of SOD genes provides insights into the evolutionary history and has laid a foundation for gene role, functional characteristics, and molecular mechanism in the plant development process and stress response (Mistry et al., 2021).

DATA AVAILABILITY STATEMENT
The datasets presented in this study can be found in online repositories. The names of the repository/repositories and accession number(s) can be found below: CuGenDB, PRJNA387615.

ACKNOWLEDGMENTS
Thank full for this grant from National Natural Science Foundation of China.