Comparative analysis of SIMILAR to RCD ONE (SRO) family from tetraploid cotton species and their diploid progenitors depict their significance in cotton growth and development

Abstract

Background

SRO (Similar to RCD1) genes family is largely recognized for their importance in the growth, development, and in responding to environmental stresses. However, genome-wide identification and functional characterization of SRO genes from cotton species have not been reported so far.

Results

A total of 36 SRO genes were identified from four cotton species. Phylogenetic analysis divided these genes into three groups with distinct structure. Syntenic and chromosomal distribution analysis indicated uneven distribution of GaSRO, GrSRO, GhSRO, and GbSRO genes on A2, D5 genomes, Gh-At, Gh-Dt, Gb-At, and Gb-Dt subgenomes, respectively. Gene duplication analysis revealed the presence of six duplicated gene pairs among GhSRO genes. In promoter analysis, several elements responsive to the growth, development and hormones were found in GhSRO genes, implying gene induction during cotton growth and development. Several miRNAs responsive to plant growth and abiotic stress were predicted to target 12 GhSRO genes. Organ-specific expression profiling demonstrated the roles of GhSRO genes in one or more tissues. In addition, specific expression pattern of some GhSRO genes during ovule development depicted their involvement in these developmental processes.

Conclusion

The data presented in this report laid a foundation for understanding the classification and functions of SRO genes in cotton.

Introduction

Changing climatic conditions adversely affects plants that directly or in-directly affect their growth, development and thus influence the economic production. Due to the sessile nature, plants cannot move, instead they have to adapt to environmental stresses by making changes at the cellular, molecular, and/or physiological levels. During evolution, plants have evolved advanced mechanisms to cope with multiple stresses. Several phytohormones i.e. salicylic acid (SA), jasmonic acid (JA), ethylene together with reactive oxygen species (ROS) and kinases participate either in synergistic or antagonistic interactions to regulate stress responses (Blazquez et al., 2020). Additionally, plant-specific proteins and transcription factors including WRKY, MYC2, NAC, ARF, AP2, and DOF (Chattha et al., 2020; Du et al., 2017; Gu et al., 2017; Singh et al., 2021; Wang et al., 2018a, 2019; Xiao et al. 2018; Zhang et al., 2020) are also involved in plant development and stress responses.

SRO (similar to radical-induced cell death one) protein family has been recently recognized, which has multiple regulatory roles in plant development and mitigates stress responses (Jaspers et al., 2010a; Jiang et al., 2018; Liu et al., 2014). The SRO genes family is a group of plant-specific proteins which are characterized by the presence of two domains, a PARP [poly (ADP-ribose) polymerase] domain and a C-terminal catalytic RST (RCD1-SRO-TAF4) domain (Jaspers et al., 2010b). Moreover, some SRO proteins also possess an N-terminal WWE domain besides of PARP and RST domains (Webb et al., 2011). It is presumed that WWE domain mediates protein-protein interactions (Aravind, 2001) while RST domain promotes interaction between RCD1 and transcription factors (Jaspers et al., 2009).

SRO proteins from Arabidopsis constitute one RCD1 (AtRCD1) and five SRO (AtSRO1-AtSRO5). AtRCD1 also known as AtCEO, is the first characterized SRO gene in Arabidopsis. It has been reported that AtRCD1 supports and enhances the ability of Saccharomyces cerevisiae to respond to oxidative stress (Belles-Boix et al., 2000). Loss of AtRCD1 function causes pleotropic effects in Arabidopsis, including developmental defects and enhancement of sensitivity to abiotic stress (Ahlfors et al., 2009; Jaspers et al., 2009; Katiyar-Agarwal et al., 2006; Teotia et al., 2009). Further, mutation studies revealed that AtRCD1 regulates several stress related genes and acts as a key regulator in phytohormonal signaling pathways (Ahlfors et al., 2004; Overmyer et al., 2000). AtRCD1 and AtSRO1 are the two closest homologs, play overlapping roles during growth and developmental processes in plants (Teotia et al., 2009). Both genes harbor an additional WWE domain along with PARP and RST domains and are categorized in group I (Jaspers et al., 2009), while other Arabidopsis SRO genes (AtSRO2-AtSRO5) only contain PARP and RST domains and are categorized in group II (Jaspers et al., 2010b). Although, AtSRO1 and AtRCD1 have functional redundancy in plant development but functional variations exist in response to abiotic stresses (Teotia et al., 2009). For example, mutations of AtSRO1 could improve plant tolerance towards osmotic and oxidative stresses (Jaspers et al., 2010a), but mutations of AtRCD1 compromise plant tolerance to salt stress (Zhao et al., 2018).

Moreover, other members of Arabidopsis SRO gene family have also been documented to increase plant tolerance or adaptability against several abiotic stresses. AtSRO2/3 participate in the tolerance of strong light, salt, and ozone stresses (Jaspers et al., 2010b), and the induction of AtSRO5 improves plant tolerance to salt stress by modulating the level of H₂O₂ in roots (Borsani et al., 2005). Apart from Arabidopsis, functional characterization of SRO genes from non-model crops including wheat (Jiang et al., 2020), rice (You et al., 2014), maize (Jiang et al., 2018), tomato (Babajani et al., 2009), cabbage (Qiao et al., 2020), and banana (Zhang et al., 2019) have also been reported. For example, one of the six SRO genes in banana, MaSRO4 was reported to be the key regulator in response to multiple stresses by interacting with NAC6 and MYB4 transcription factors (Zhang et al., 2019). Similarly, overexpression of MdRCD1 increases root growth and improves plant’s ability to tolerate salt, oxidative, and drought stresses by regulating abscisic acid (ABA) signaling pathway (Li et al., 2017). Moreover, BrSRO8 strongly responds to heat stress in cabbage (Qiao et al., 2020) and OsSRO1c gene in rice increases plant tolerance to oxidative and drought stresses by interacting with transcription factors SNAC1, DST and modulating H₂O₂ level (You et al., 2014). Interestingly, SlSRO1 is closely related to AtSRO5 based on the high sequence similarity, and is strongly induced under salt stress and regulates stress responses in tomato (Babajani et al., 2009). Previously, it has also been reported that miRNAs, such as miR444 which is known to involve in the regulation of plant reproductive development, also regulates the expression of two SRO genes (TaSRO2b.3-4 A and TaSRO2b.5-6B) in wheat (Jiang et al., 2020).

Cotton, a natural fiber producing crop, provides major raw material to the textile industry. Sustainable production of cotton is seriously affected by a number of biotic and abiotic stresses during its growth period, such as flooding, drought, heat wave, and pathogen infection (Mahmood and Hussain, 2020). With the increment and advancement of transcriptomic and genomic data of cotton species, it has greatly facilitated the characterization of the functions of specific genes in cotton. Considering the important roles of SRO gene family during the growth and development in plants, comparative analysis of SRO genes among cotton species is worth to be studied in depth.

Here, genomes from four cotton species have been screened for the identification of putative SRO genes. Next, chromosomal mapping, gene duplication, syntenic relationship, phylogeny relationships, conserved domains, exon-intron structure, motifs distribution, cis-elements, and miRNAs have been predicted to dissect the potential roles of GhSRO genes. Moreover, transcript abundances of GhSRO genes in various organs and within different ovule developmental stages have been analyzed to explore the function of potential candidate genes. Our results provide useful information and lay foundation in deciphering the specific roles of SRO genes in cotton, which will be helpful for the growth and stress-related breeding research in cotton.

Materials and methods

Screening and sequence retrieval of SRO genes from cotton species

Firstly, all the reported SRO genes from Arabidopsis were downloaded (Jaspers et al., 2009). Secondly, BLASTP program available on CottonFGD platform (https://cottonfgd.org/) (Zhu et al., 2017) was employed with e-value of 1e⁻¹⁰ to extract sequences of candidate SRO from Gossypium hirsutum (JGI_v1.1), G. barbadense (P-90_HEAU_v1), G. arboreum (CRI-v1.0-a1.0), and G. raimondii (JGIJGI_v2_a2.1) using all the reported AtSRO genes as queries. Further, protein sequences of cotton SRO genes were checked using Pfam (http://pfam.xfam.org/) (El-Gebali et al., 2019) and NCBI conserved domain databases (https://www.ncbi.nlm.nih.gov/Structure/cdd/wrpsb.cgi) (Yang et al., 2020). All redundant sequences without potential domains were manually discarded. The complete genome, proteome, and coding sequences of putative SRO genes of cotton species were obtained through browsing CottonFGD website (https://cottonfgd.net/analyze/) (Zhu et al., 2017) and nomenclature of candidate genes was based on their physical positions on chromosomes. Physicochemical properties of putative SRO genes were obtained from online webtool ExPasy (http://web.expasy.org/protparam/) (Bjellqvist et al., 1994) and CottonFGD website (https://cottonfgd.org/) (Zhu et al., 2017). The CELLO v2.5 (http://cello.life.nctu.edu.tw/) (Yu et al., 2006) and WOLF PSORT (https://www.genscript.com/wolf-psort.html?src=leftbar) (Horton et al., 2007) software were employed to infer the subcellular localization of cotton SRO genes.

Chromosomal mapping, duplication and syntenic analysis of cotton SRO genes

Information obtained from CottonFGD database (https://cottonfgd.org/) (Zhu et al., 2017) was utilized to find out the chromosomal positions of putative cotton SRO genes. Their positions were visualized on chromosomes using PhenoGram webtool (http://visualization.ritchielab.psu.edu/phenograms/plot). Duplication and syntenic analysis were performed using online BLAST program available at CottonFGD website (https://cottonfgd.org/sequenceserver/) (Zhu et al., 2017). Criteria opted for duplication study among homologs of four cotton species (G. hirsutum, G. barbadense, G. arboreum, and G. raimondii) was similarity > 80% and alignment percentage > 80% of full length proteins (Shaban et al., 2021). For visualization of syntenic relations, a circos plot was generated among the four studies species of cotton using Circos program (TBtools software) (Chen et al., 2020). Clustal omega (http://www.ebi.ac.uk/Tools/msa/clusalo/) was used for sequence alignment and aligned sequence file was submitted to PAL2NAL webtool (http://www.bork.embl.de/pal2nal/) (Suyama et al., 2006) for determining of Ka, Ks and Ka/Ks values.

Phylogenetic analysis of SRO genes

Protein sequences from four species of cotton, Arabidopsis (Jaspers et al., 2009), maize (Jiang et al., 2018), tomato (Li et al., 2021), and apple (Li et al., 2017), were aligned using clustalW software (Thompson et al., 2003) with default parameters. Alignment file was submitted in MEGA7.0 software (Kumar et al., 2016) for creation of unrooted tree using both neighbor joining and maximum likelihood methods with 1000 replications.

Structural, motifs and promoter analysis of GhSRO genes

Gene structural diagram was generated in Gene Structural Display Server (GSDS 2.0)( http://gsds.cbi.pku.edu.cn/) (Hu et al., 2015) using the full length protein sequence and the corresponding DNA sequences. Potential conserved motifs in GhSRO, GbSRO, GrSRO, and GaSRO proteins were explored using online MEME software with default setting (Bailey et al., 2015). A 1 500 bp (base pair) upstream promoter region from start codon of each GhSRO gene was obtained from CottonFGD website. The potential cis-elements in these promoter regions were identified via PlantCare webtool (http://bioinformatics.psb.ugent.be/webtools/plantcare/html/) (Lescot et al., 2002).

Prediction of putative miRNAs and their target GhSRO genes

Putative miRNA sequences were obtained from various miRNA databases, including plant miRNA database (http://bioinformatics.cau.edu.cn/PMRD/), the cotton EST database (http://www.leonxie.com/), the miRBase (http://www.mirbase.org/) and published articles. Targets of miRNAs in GhSRO genes were predicted using the online psRNATarget server (http://plantgrn.noble.org/psRNATarget/home) with default parameters (Dai et al., 2018).

Transcriptomic data and expression analysis of GhSRO genes

In the current study, publically available RNA-sequencing data related to various tissues (root, stem, leaf, sepal, petal, anther, filament, pistils, bract, and torus) and in various ovule developmental stages (0, 1, 3, 5, 10, 15, 20, and 25 DPAs) of G. hirsutum accession TM-1(Bioproject PRJNA490626) were downloaded from CottonFGD (https://cottonfgd.net/analyze/) (Zhu et al., 2017). The normalized relative expression of candidate genes were analyzed, and the differentially expressed genes and hierarchical clustering was ploted based on pearson coeffecient method in Genesis software (Version 1.7.7) (Sturn et al., 2002) as described by Zhang et al. (2021). Further, digital data were visualized in the form of heatmap using Genesis software (Version 1.7.7) (Sturn et al., 2002) as proposed by Ren et al. (2017).

Results

Identification, classification and properties of SRO genes of cotton species

The Arabidopsis SRO genes were used as query in BLASTP search against whole genomic databases of four species of cotton. As a result, a total of 36 SRO genes were identified, including 12, 11, 6 and 7 SRO genes from G. hirsutum, G. barbadense, G. arboreum, and G. raimondii, respectively (Table S1). All these putative SRO genes were checked for characteristic SRO domains using Pfam and NCBI conserved domain databases. According to the composition of characteristic SRO domains, all the putative SRO genes were classified into three groups. Group I included genes consisted of WWE, RST super family, and PARP super family domains. Group II genes lacked WWE domain but contained other two domains, while Group III harbored only the RST super family domain. Among 36 SRO genes, 13 genes belonged to group I, 11 genes belonged to group II and 12 genes belonged to group III (Table S2). The detailed information related to their physio-chemical characteristics was determined using ExPASy and CottonFGD webtools. The number of amino acids of SROs ranged from 305 to 616, molecular weight varied between 34.263 kDa and 68.973 kDa, the isoelectric point spaned from 6.379 to 9.077, and the grand average of hydropathy fluctuated between − 0.479 and − 0.215 (Table S2).

Phylogeny and evolutionary analysis of cotton SRO genes

To explore the phylogeny among cotton SRO genes with other well characterized SRO genes from Arabidopsis, maize, tomato, and apple, an unrooted tree was generated using their protein sequences. In this phylogenetic tree, all 45 SRO genes were scattered among three groups (Group I, Group II, and Group III) which consisted of 18, 12, and 15 genes, respectively (Fig. 1). Cotton SRO genes were randomly distributed in all three groups. Among the GhSRO genes, GhSRO12, GhSRO6, GhSRO10, and GhSRO4 were clustered in group III with most of the Arabidopsis SRO genes (AtSRO2, AtSRO3, AtSRO4, and AtSRO5). Four GhSRO genes (GhSRO8, GhSRO2, GhSRO7, and GhSRO1) were clustered in group II. Maximum number of SRO genes were clustered in Group I, including GhSRO11, GhSRO5, GhSRO3, and GhSRO9 from upland cotton, two Arabidopsis SRO genes (AtRCD1 and AtSRO1), maize SRO gene (ZmSRO5), apple SRO gene (MdRCD1), and tomato SRO gene (SlSRO2). As seen in Fig. 1, many orthologous genes from four cotton species were clustered in one branch. Seven orthologous pairs were found between G. hirsutum and G. barbadense (GhSRO6/GbSRO6, GhSRO8/GbSRO8, GhSRO7/GbSRO7, GhSRO1/GbSRO1, GhSRO11/GbSRO10, GhSRO3/ GbSRO3, and GhSRO9/GbSRO9). Two orthologues were found between G. arboreum and G. barbadense (GaSRO6/GbSRO4 and GaSRO2/GbSRO2), One orthologue between G. raimondii and G. hirsutum (GrSRO3/GhSRO10) and one orthologue between G. arboreum and G. hirsutum (GaSRO5/GhSRO5). Moreover, GaSRO1/GrSRO6, GaSRO2/GrSRO7, GaSRO3/GrSRO1, GaSRO4/GrSRO5, GaSRO5/GrSRO4, and GaSRO6/GrSRO3 were also considered orthologous pairs with more than 94% amino acid sequence identity (Table S3).

Fig. 1

Phylogenetic analysis of cotton SRO genes. The maximum likelihood phylogenetic tree among SRO genes of G. arboreum, G. raimondii,G. barbadense, G. hirsutum, Arabidopsis, maize, potato, and apple. The tree was constructed using MEGA 7.0 (1 000 bootstrap value) using the full length amino acids sequences of SRO genes. Three groups were represented with different color

Chromosomal distribution and syntenic study of cotton SRO genes

To explore the evolutionary dynamics and syntenic relations among cotton SRO genes, a circos plot was generated among two diploid cotton species and two allotetraploid cotton species. Results showed that putative GaSRO, GrSRO, GhSRO, and GbSRO genes were unevenly distributed among three chromosomes of A2 (Ga-A05, Ga-A08, Ga-A12), D5 (Gr-D04, Gr-D08, Gr-D09), At (Gh-A05, Gh-A08, Gh-A12), Dt (Gh-D05, Gh-D08, Gh-D12), At (Gb-A05, Gb-A08, Gb-A12), and Dt (Gb-D05, Gb-D08, Gb-D12) genomes, respectively (Fig. 2). The gene number varied from one to four, with one gene on Ga-A08, Gr-D04, Gh-A08, Gh-D08, Gb-A08, and Gb-D08, two genes on Ga-A05, Gr-D09, Gh-A05, Gh-D05, Gh-D12, Gb-A05, Gb-D05, Gb-D12, and three genes on Ga-A12, Gh-A12, Gh-D12, and Gb-A12. The maximum number of genes (four) were found on Gr-D08 chromosome (Table S4, Fig. 2). Gene duplication analysis of GhSRO genes revealed six duplicated pairs which shared more than 95% similarity in nucleotide sequences. These pairs are GhSRO1-7, GhSRO2-8, GhSRO3-9, GhSRO4-10, GhSRO5-11, and GhSRO6-12 (Fig. S1). In order to investigate the duplication mechanism, synonymous (Ks), non-synonymous (Ka) substitution rates and respective ratios were calculated. Ka and Ks values of GhSRO genes were in the range of 0.0129 to 0.0588 and 0.0154 to 0.0461, respectively. A Ka/Ks ratio equal to 1 signifies neutral selection, whereas a ratio greater than 1 suggests positive selection, and a ratio less than 1 indicates purifying selection (Shaban et al., 2021). Interestingly, all the duplicated gene pairs underwent segmental duplication with the purifying selection pressure (Ka/Ks < 1) during evolution (Table S5).

Fig. 2

Synteny analysis of cotton SRO genes. Syntenic relationship among SRO genes of two diploid (G. arboreum, G. raimondii) and two allotetraploid species (G. hirstutum, G. barbadense) of cotton was shown in the form of circos plot. The chromosomes from A and D sub-genomes were highlighted with blue and green, respectively

Structure and motif analysis among SRO genes of cotton species

Gene structural analysis comprising of introns/exons distribution and numbers showed that Group I members had the highest number of exons and introns (Exon = 6 and Intron = 5), except GrSRO2 which comprised of 5 exons and 4 introns. All members of Group II and Group III possessed 4 exons and 3 introns, while two members of Group III contained 5 exons and 4 introns (Fig. 3A). Collectively, structural diagram depicts the conservation in number and distribution of exons/introns among closely related SRO genes within the group.

Conserved domain/motif analysis provides clues about gene duplications and functional conservation during evolution. Using online MEME tool, seven motifs were found among SRO proteins of four cotton species and their distributions were displayed in Fig. 3B. Different SRO proteins contained different motifs, except for motif 1, which was conserved among all members of cotton SRO proteins. Moreover, cotton SRO genes within subgroups harbored similar motifs. More intriguingly, genome specific homologs of GhSRO genes, such as GhSR01/GhSR07 and GhSRO2/GhSRO8 shared similar clustering pattern, exon/intron distribution and motif conservation.

Fig. 3

Structure and motif distribution of Cotton SRO genes/proteins. A Phylogenetic tree and exon/intron structure among cotton SRO proteins, un-rooted tree using neighbor-joining method was constructed with protein sequences using MEGA 7.0 software. Exon/intron structure, length and distribution were analyzed by comparing genomic and cDNA sequences of cotton SRO genes. B Position of seven specific motifs was schematically displayed using online MEME tool

Cis-element analysis

Promoter regions of every gene possess certain elements that regulate gene expression. These regulatory elements are known as cis-elements. To analyze specific functions of GhSRO genes during cotton development and in response to different stress factors, we comprehensively investigated the cis-regulatory elements of GhSRO genes. As shown in Fig. 4A, mainly four types of elements were found in GhSRO genes, including light responsive, growth or development related, hormones-responsive, and stress related. Light responsive elements were predominant (59%) in GhSRO genes, followed by hormones (17%), stress related elements (17%), and plant growth responsive (7%). Multiple cis-elements were found in each GhSRO genes and the patterns of these cis-elements varied among GhSRO genes. Comparing with other GhSRO genes, GhSRO4 had the maximum number of low temperature-responsive, development-related, and GA-related elements. Interestingly, only GhSRO5 and GhSRO7 had cis-elements related to defense and stress responses (Fig. 4B).

Fig. 4

Identification of cis-elements at the promoter region of GhSRO genes. A The percentage of light, hormones, stress and growth/development related cis-elements of all GhSRO genes. B Types of cis-elements in each GhSRO gene was represented in the form of bar graphs

Identification of potential miRNA targeting sites in GhSRO genes

MiRNAs are 20–25 nt small non-coding RNAs known to participate in certain regulatory function on target genes in eukaryotes (Hua et al., 2018; Wang et al., 2018b). Recently, the regulatory roles of these small RNAs in plant developmental processes and stress responses were widely studied (Ma et al., 2015; Wei et al., 2013; Xie et al., 2015; Xie and Zhang, 2015). To further elucidate the function of GhSRO genes in cotton and to predict potential miRNAs targeting sites, coding DNA sequences of GhSRO genes were submitted in psRNATarget server (http://plantgrn.noble.org/psRNATarget/home). The results revealed 23 miRNAs targeting 12 GhSRO genes in upland cotton (Fig. 5, Table S6). The predicted miRNAs were related to fiber development, growth stages and environmental stresses in plants (Wang et al., 2020). These miRNAs were identified in numerous RNA-sequencing and bioinformatics based studies. Some identified miRNAs were validated through experiments, while some were only predicted without wet-lab validation. The results showed that ghr-mir2868 targeted GhSRO2, GhSRO8 genes and target sites was present at RST domain. Also, ghr-mir482d, ghr-miR414f, ghr-n72, ghr-miR2595, ghr-mir482a, and ghr-miR482h targeted GhSRO1, GhSRO3, GhSRO5, GhSRO7, GhSRO9, and GhSRO11, respectively. Generally, homologs of GhSRO were targeted by the same type of miRNAs. For example, GhSRO4, GhSRO10, GhSRO6, and GhSRO12 were targeted by ghr-miR2529 and novel-miR-27, respectively. The complete details comprising predicted miRNAs and their potential targets in GhSRO genes were provided in Table S6.

Fig. 5

Prediction of miRNA-mediated targeting site of GhSRO genes. Light grey box represented the total length of CDS of each GhSRO gene. WWE, PARP and RST domains were shown with different colored boxes. Light dark color filling (boxes) showed the target position of miRNA in each GhSRO gene. Sequences and length of each miRNA and their targeted RNA sequence in GhSRO genes were also displayed along the CDS of each GhSRO gene

Expression profiling of GhSRO genes in various tissues and ovule developmental stages of cotton

Expression patterns of genes partly reflect their functions. SRO genes had been reported to regulate various developmental processes in plants (Jaspers et al., 2010b; Zhang et al. 2019). To preliminary investigate the biological functions of GhSRO genes in various tissues and ovule development, we analyzed transcriptomic data of GhSRO genes (Zhu et al., 2017). Differential expression pattern from all the candidate GhSRO genes were observed in most of the studied tissues (roots, stem, leaf, sepal, petal, anther, filaments, pistils, bract, and torus). Of the total 12 candidate genes, six GhSRO genes (GhSRO6, GhSRO4, GhSRO12, GhSRO10, GhSRO7, and GhSRO1) were predominantly induced in anthers. Among other GhSRO genes, GhSRO3 was highly induced in roots, while GhSRO8 showed higher expression in sepals and torus, lower expression in anthers and filaments, and barely detectable expression in roots and bracts (Fig. 6A).

As shown in Fig. 6B, all candidate SRO genes induced differentially during various ovule developmental stages. Two genes (GhSRO4 and GhSRO10) more specifically induced in ovules at the later developmental stages, while two genes (GhSRO8 and GhSRO2) with higher expression during mid-developmental stages of ovules, while others were differentially induced during initial and later stages of ovules development, such as GhSRO11, GhSRO5, GhSRO3, and GhSRO9. Interestingly, homologs of GhSRO genes depicted similar expression pattern in different tissues and during ovule development. Such as, GhSRO3/GhSRO9 and GhSRO5/GhSRO11 were preferentially induced in root and stem tissues and showed similar expression during initial and mid ovule developmental stages. Remarkably, GhSRO10 most preferentially induced during ovule development at 25 DPA, indicating its specific function at this stage (Fig. 6B). Collectively, our expression profiling data suggested the important roles of GhSRO genes at various growth and developmental stages of cotton.

Fig. 6

Expression pattern of GhSRO genes. A Tissue-specific expression pattern. B Expression pattern during different ovule developmental stages. Transcriptomic data related to tissue specificity and ovule developmental stages were accessed from CottonFGD database and used for generation of heatmap

Discussion

Cotton being a natural fiber and oil producing crop contributes a major share in the economy of several developing countries. Due to global climate change, various environmental stresses have become a major obstacle in sustainable production of cotton (Jans et al., 2021; Rahman et al., 2018). Therefore, it is very important to discover stress responsive candidate genes under changing environmental conditions for their best utilization in cotton breeding programs. Availability of genomics, transcriptomic data of important cotton species and advancement of bioinformatic tools have provided the opportunity to characterize key genes that participate in growth/developmental processes and regulation of various stress responses in plants.

SRO is a unique family of small proteins that regulate multiple growth related processes in plants (Qiao et al., 2020). However, cotton SRO genes have not been comprehensively characterized before. Here, comprehensive analysis of SRO genes from four species of cotton, including two diploid (G. arboreum, G. raimondii) and two allotetraploid (G. barbadense, G. hirsutum) were performed, with the aim to explore their diverse functions in cotton growth and developmental processes.

The identified SRO genes in tetraploid species is roughly twice to that of diploid species, this suggests that duplication of SRO genes in diploid species occur before the emergence of tetraploid cotton. In addition, variations in the number of SRO genes among diploid and tetraploid species corresponds to variations in genome size of these species, consistent with the study by Chattha et al. (2020). Domain analysis of 36 putative cotton SRO genes has revealed that 13 cotton SRO genes harbor all three characteristic SRO domains (PARP, RST, WWE), 11 SRO genes possesses two domains (PARP, RST), while 12 SRO genes constitute only RST domain. In general, presence or absence of these domains is somehow representative of their functional similarity or diversification. As previous studies reported that PARPs are specific enzymes that regulate various processes in plants including chromatin modification, transcription, DNA damage repair and cell death pathways (De Block et al., 2005; Kim et al., 2005; Lamb et al., 2012). RST domain permits interactions between RCD1 and other transcription factors to regulate important developmental processes in plants (Jaspers et al., 2009), while WWE domain is thought to be required in facilitation of interactions between proteins by forming specific globular structure (Aravind, 2001). The specific domain architecture of cotton SRO genes is representative of their special functions in cotton growth and developmental processes.

In phylogenetic tree, cotton SRO genes of four species are randomly distributed in three groups (Group I, Group II, Group III) along with SRO genes from other plant species including Arabidopsis, maize, tomato, and apple. One way to predict the functional similarity of genes is to find the close association through phylogenetic tree. According to phylogenetic analysis, Group I cotton SRO genes are found to be more related to AtSRO1, AtRCD1, MdRCD1, and ZmSRO5. In this group, AtSRO1 and AtRCD1 which have been studied comprehensively, have been reported to regulate several growth/developmental processes and stress responses in Arabidopsis (Jaspers et al., 2009; Song et al., 2019). MdRCD1 in apple and ZmSRO5 in maize have multiple functions related to root growth and regulate the tolerance against several abiotic stress factors (Jiang et al., 2018; Li et al., 2017). Group III cotton SRO genes showed more resemblance with four Arabidopsis genes (AtSRO2, AtSRO3, AtSRO4, and AtSRO5), depicting their functional similarity.

Gene distribution pattern of cotton SRO genes on chromosomes is uneven. Out of 13 chromosomes, all the GaSRO, GrSRO, GhSRO, and GbSRO genes are distributed only across three chromosomes of A2, D5 genomes and Gh-At, Gh-Dt, Gb-At, and Gb-Dt subgenomes, respectively. Gene duplication or replication is an important evolutionary mechanism for organisms to expand gene family and adopts novel gene functions. Segmental and tandem duplication are the ways through which genes are more frequently duplicated in plants. Replicated gene pairs located on the same chromosome represent the tandem duplication, whilst duplicated gene pairs from different chromosomes are consider under segmental duplication event (Song et al., 2019). In the current study, segmental duplication has been found to be the main cause for expansion of GhSRO gene family in cotton, which corroborates the previous findings (Chattha et al., 2020; Shaban et al., 2021). Generally, synonymous (Ks) and non-synonymous (Ka) substitution values are computed to explain the duplication mechanism (Hurst et al., 2002). All the duplicated pairs of SRO genes have Ka/Ks ratio < 1, meaning have undergone the purifying selection pressure but not much functional diversification during evolution.

Comparative structural analysis including exon/intron architecture and predicted motifs among cotton SRO genes has revealed the variations between subgroups and sequence conservation within subgroups. This further supports the phylogenetic clustering and evolutionary relationship among closely related cotton SRO genes. High structural conservation of SRO genes within groups has also perceived in other plant species (Jiang et al., 2018; 2020). Promoter regions located cis-elements fine tune gene expression and corresponding analysis provide clues about functional characteristics of genes (Hernandez-Garcia and Finer, 2014). PlantCARE database contains a collection of previously reported cis-elements in plants (Lescot et al., 2002), helpful to investigate promoter regions of candidate GhSRO genes in details. The discovery of several growth/development related, light responsive, hormones responsive, and stress responsive cis-elements in the promoter regions of GhSRO genes predicts their multiple functions in different biological processes of plants.

Recently, a number of miRNAs are characterized for their potential functions in post-transcriptional regulation during development and stress related responses in plants (Hua et al., 2018; Wang et al., 2018b). SRO genes have been extensively characterized of their role at various growth and developmental stages and in stress responses in multiple plant species (Jiang et al., 2018; Qiao et al., 2020; Zhang et al., 2019), and their regulation mediated by miRNAs requires futher detailed analysis. In this study, 23 miRNAs of upland cotton have been found to target 12 GhSRO genes. These miRNAs, for example ghr-miR482d, ghr-miR414f, ghr-n72, and novel_miR_27 regulate various developmental processes and are responsive to abiotic and biotic stresses in cotton (Wang et al., 2020; Xie et al., 2015; Zhu et al., 2013). Our results have revealed that most closely related genes are targeted by the same miRNAs. The above mentioned miRNAs is identified based on computational approach, further experimental validation of these miRNAs' regulatory roles on their target genes using molecular techniques are needed.

Generally, the pattern of gene expression predicts their functions. Expression pattern of GhSRO genes in different organs and at various ovule developmental stages provide clues for their functional specificity in cotton growth and development. Previously, the SRO genes have been found to be tissue-specific in multiple plant species, such as maize (Jiang et al., 2018), wheat (Jiang et al., 2020), and banana (Zhang et al., 2019), corresponding to their functional specificity in these species. In this study, half of the candidate GhSRO genes are dominantly expressed in anthers, signifying their functional redundancy in cotton reproductive development. While the other half of GhSRO genes have showed differential expression in the studied tissues, implying their functional diversity across different tissues. Interestingly, GhSRO4 and GhSRO10 are dominantly induced during later stages of ovule development, GhSRO8 and GhSRO2 are induced during mid-developmental stages and some GhSRO genes are induced during initial and later stages of ovule development indicating the spatio-temporal roles in ovule development. More intriguingly, GhSRO genes expression pattern analysis justify the phylogenetic clustering. GhSRO5, GhSRO11, GhSRO3, and GhSRO9 clustered in the same group in the phylogenetic tree have shown similar expression pattern in various tissues and at ovule developmental stages. Collectively, GhSRO genes expression pattern results suggest that GhSRO genes have important regulatory roles at various growth and developmental stages in cotton.

Conclusion

In conclusion, a total of 36 putative SRO genes from four species of cotton are distributed among three chromosomes of A2, D5, and At and Dt genomes, respectively. SRO genes in upland cotton have been found to expand through segmental duplication during evolution. Genes belonging to same groups or subgroups have revealed similar phylogeny, gene structure, and motif distribution. Wide variety of cis-elements in the promoter region of GhSRO genes suggest their roles at various growth and developmental stages of cotton. The transcriptomic data based expression pattern reveals the differential expression of GhSRO genes in various tissues and during ovule development. This provides evidence of the possible involvement of specific GhSRO genes in cotton growth and development. The data presented here provides fundation for the further exploration of the specific function of GhSRO genes in cotton growth and development.