Genome-wide identification, evolution, expression, and alternative splicing profiles of peroxiredoxin genes in cotton

Sequence sources

G. raimondii (Gossypium raimondii L.,JGI) (Accession: PRJNA171262) (Paterson et al., 2012) and G. hirsutum (Gossypium hirsutum L., NAU) (Accession: PRJNA248163) (Zhang et al., 2015) genomic data files were obtained from the JGI database (https://genome.jgi.doe.gov/portal/); G. arboreum (Gossypium arboreum L., CRI) (Accession: PRJNA382310) (Du et al., 2018) and G. barbadense (Gossypium barbadense L., Hau) (Accession: PRJNA433615) (Wang et al., 2018) genomic data files were downloaded from the CottonGen database (https://www.cottongen.org/).

The following protein sequence data were obtained from the corresponding databases. Arabidopsis (Arabidopsis thaliana L.) (Accession: SRA009031) (Filichkin et al., 2010) protein sequence data were obtained from The Arabidopsis Information Resource (http://www.arabidopsis.org). Rice (Oryza sativa L.) (Ouyang et al., 2007) protein sequence data were downloaded from the Rice Genome Annotation Project (http://rice.plantbiology.msu.edu/index.shtml). Cacao (Theobroma cacao L.) (Accession: PRJNA51633) (Motamayor et al., 2013) protein sequence data were obtained from the JGI database (https://genome.jgi.doe.gov/portal/), and grapevine (Vitis vinifera L.) (GenBank: CU459218 –CU462737) (Jaillon et al., 2007) protein sequence data were obtained from the Ensembl Plants database (http://plants.ensembl.org/index.html).

Identification and conserved active site analysis of the PRX family in cotton

The PRX Pfam domain ids (http://pfam.xfam.org/family/PF08534.9, http://pfam.xfam.org/family/PF00578.20, http://pfam.xfam.org/family/PF10417.8) in cotton were identified using the PRX protein sequence of G. hirsutum on the EMBL-EBI (https://www.ebi.ac.uk/Tools/hmmer/) website. The Pfam model files were downloaded from the Pfam database (http://pfam.xfam.org/) (Finn et al., 2009). The protein databases of four cotton species (G. raimondii, G. arboreum, G. hirsutum, and G. barbadense) and Arabidopsis, Rice, Cacao, and Grapevine were searched with hmmsearch (v3.2.1) (https://www.ebi.ac.uk/Tools/hmmer/search/hmmsearch) (E <  = 0.001). All possible PRX genes were identified in these eight crop species. The BLAST program was used to compare the sequences to the Arabidopsis protein database (p <  = 1E⁻¹⁰), and then, the false positive PRX genes were deleted on the basis of the Arabidopsis protein annotation file (https://www.arabidopsis.org/). To further ensure the accuracy of each candidate PRX, ClustalW in the MEGA7 (https://www.megasoftware.net/) software was used to compare the candidate PRX protein sequences with the conserved active sites and positions of conserved cysteinyl residues in the PRX proteins in the PRX protein database (http://csb.wfu.edu/prex/) (Dietz, 2011; Nelson et al., 2011) and delete sequences with no conserved active site (PxxxTxxC…S…W/F). Blastp alignment of the identified PRX proteins from the NCBI reference protein database was performed (percent identity ≥90%). The theoretical isoelectric point (pI) and molecular weight (MW) of the PRX proteins were investigated with ExPASy (http://web.expasy.org/protparam/) (Finn et al., 2014).

Phylogenetic and gene structure analysis

ClustalW in MEGA 7.0 software (Kumar, Stecher & Tamura, 2016) was used to compare the protein sequences, and a phylogenetic analysis was then performed. A phylogenetic tree was constructed in MEGA7 using the maximum likelihood (ML) method, the Poisson correction model, complete deletion and bootstrap analysis performed with 1,000 replicates.

The exon/intron structures of the PRXs were extracted from the corresponding cotton GFF files by TBtools software (http://www.omicshare.com/forum/thread-1062-1-1.html), and the PRX structure maps were drawn by GSDS 2.0 (http://gsds.cbi.pku.edu.cn/) (Hu et al., 2014). MEME (4.11.4) (http://meme-suite.org/) (Bailey et al., 2006) was used to identify the conserved motifs of the cotton PRXs with the following parameters: the maximum number of motifs was 10 and the optimal width was 6 ≤250. Then, TBtools was used to construct a vector graph from the xml file generated by MEME.

Chromosomal mapping and gene duplication

To reveal the chromosomal mapping and duplication relationships of PRXs in cotton, genome databases for G. raimondii, G. arboreum, G. hirsutum and G. barbadense were constructed. The physical positions of PRXs in cotton were fetched from the corresponding GFF files. PRX family duplication data were identified by the MCScanX program (p ≤ 1E⁻²⁰) (Ding et al., 2015), and a comparative genome analysis of the PRX family genes in these four cotton species was carried out (p ≤ 1E⁻²⁰) (Wang, Li & Paterson, 2013). Visualization was carried out with the CIRCOS (http://circos.ca/) tool (Krzywinski et al., 2009). The substitution rates of synonymous (Ks) and nonsynonymous (Ka) sites were calculated by the KaKs calculator program (Suyama, Torrents & Bork, 2006). The divergence time was calculated by the formula T=Ks/2^λ (λ = 1.5 × 10⁻⁸) (Zhang et al., 2006), where Ks was the synonymous substitution of each locus and r was the divergence rate of plant genes. For dicotyledonous plants, r is considered to be 1.5 × 10⁻⁸ synonymous substitutions per site per year (Koch, Haubold & Mitchell-Olds, 2000).

Analysis of cis-acting elements in the promoter region

By using the genome data files (GFF3) of the four cotton species, the promoter sequences of the PRX genes (2,500 bp upstream of the initiation codon “ATG”) were extracted from the cotton genome sequences (Wang et al., 2012). PlantCARE (http://bioinformatics.psb.ugent.be/webtools/plantcare/html/search_CARE.html) (Lescot et al., 2002) was used to predict the cis-acting elements of the promoter sequences, and the abiotic stress response, plant growth and development, and then the hormone-induced cis-acting elements were analyzed (Table S2).

Based on the cis-acting elements identified in the promoter sequence, we used GSDS 2.0 (http://gsds.gao-lab.org/) (Hu et al., 2014) to display the physical sites of the cis-acting elements in the G. hirsutum promoter sequence.

Plant material treatments and expression analysis

To analyze the expression patterns of the PRX genes, PEG, NaCl and salicylic acid (SA) were used to simulate drought stress, osmotic stress and SA hormone induction, respectively. The seeds of G. hirsutum sGK9708 were soaked in flasks overnight, germinated in fine sand at 28 °C in the dark for 2 days, and transferred to a greenhouse for hydroponic growth. The hydroponic greenhouse conditions were as follows: 28 °C day/25 °C night, 14-hour photoperiod and 70% relative humidity. Plants were grown in Hoagland nutrient solution (Xing et al., 2019). At the three-leaf stage, the seedlings were subjected to stress treatments. The cotton seedlings were equally divided into four groups and transferred to nutrient solutions of the hydroponic box supplemented with 200 mM sodium chloride (NaCl), 15% PEG-6000, 0.1 mM SA and a blank control for the salt, drought, hormone induction and control treatments, respectively. A total of 20 cotton seedlings received each treatment, and each treatment was repeated three times. At 0, 1, 3, 6, and 12 h of treatment, four seedlings were randomly selected from four treatment groups at each time point, and the leaves (second true leaf) and rhizomes were removed from the seedlings and immediately frozen in liquid nitrogen. The leaves were stored at −80 °C before RNA extraction.

Total RNA was isolated using an EasyPure Plant RNA Kit (TransGen, Beijing, China), and RNA reverse transcription was performed using a TransScript reverse transcription system (AT341). Gh-PRX-specific primers (Table S3) were used to find candidate-specific primers on the qPrimerDB-qPCR Primer Database (https://www.ncbi.nlm.nih.gov/) website. To ensure the specificity of the Gh-PRX gene primers, the candidate-specific primers were subjected to BLAST homologous comparison in the Primer-BLAST database of NCBI (https://www.ncbi.nlm.nih.gov/) (National Center for Biotechnology Information) to ensure that the specific primers amplified only their respective target gene fragments (Table S3). The general fluorescent dye mixture used was SYBR Green I (TransStart), and the reactions were performed with a 7500 Rapid Real-time PCR system (Roche). Ubiquitin 7 (UBQ7) (GenBank: DQ116441) was used as the internal standard reference to measure the expression levels of the cDNA genes. The volume of each reaction was 10 µl, and the reaction conditions were as follows: 94 °C for 5 min and 40 cycles of 94 °C for 5 s, 55 °C for 30 s and 72 °C for 30 s. During the extension step, the fluorescence signal was measured, and the acquisition time was set to 30 s. Each cDNA sample was repeated three times, and the results were analyzed by the 2^−ΔCT method (Khan-Malek & Wang, 2011).

Analysis of tissue expression and AS profiles of PRXs

We used published G. hirsutum (TM-1) (Zhang et al., 2015) transcriptome data, including data from the roots, stems, leaves, cotyledons, petals, stamens, ovules, seeds, and 5, 10, 15 and 20 DPA fibers. The expression levels of Gh-PRX genes in various tissues were calculated using log2 (FPKM) values (Zhang et al., 2015). The expression values were normalized by Genesis software and illustrated with a heatmap by HemI 1.0 - Heatmap illustrator software (Sturn, Quackenbush & Trajanoski, 2002).

Samples of G. hirsutum sGK9708 RNA without treatment at −80 °C were obtained. Total RNA was isolated using the EasyPure Plant RNA Kit (TransGen, Beijing, China), and RNA reverse transcription was performed using the TransScript reverse transcription system (AT321). The obtained first-strand cDNA was used for subsequent RT-PCR amplification. The primers for synthesizing full-length cDNA were designed using Primer Premier 5.0 software (Premier Biosoft International, Palo Alto, CA). After RT-PCR, the amplification products were detected by 1.2% agarose gel electrophoresis. After recovery and purification, the PCR products were ligated with the pEASY-Blunt cloning vector (TransGen, Beijing, China) (The pEASY®-Blunt cloning kit was purchased from TransGen Biotech, catalog number: CB101-01) (https://www.transgenbiotech.com/) and transformed into E. coli (Escherichia coli)-competent DH5 α cells (TransGen, Beijing, China). Thirty individual colonies were selected for positive PCR identification. The positive colonies were sequenced (Sangon, Shanghai), and the full-length sequences were obtained. These sequences were compared by Splign (https://www.ncbi.nlm.nih.gov/sutils/splign/splign.cgi) and GSDS 2.0 (http://gsds.gao-lab.org/), and the AS events were displayed on the gene structure.

Subcellular localization of PRX proteins

The TargetP 1.1 (http://www.cbs.dtu.dk/services/TargetP/) server website (Emanuelsson et al., 2007) was used to predict and analyze the localization of G. hirsutum PRX proteins in plant cells.

To verify the distribution of PRXs in plant cells, five pairs of homologous genes from among 15 Gh-PRX genes were selected, and the full-length coding sequences (CDSs) of the PRX genes with no terminator codon and restriction site were amplified and subcloned into the pCAMBIA2300-35S-eGFP (CAMBIA) transient expression vector (The pCAMBIA2300-35S-eGFP vector was assembled by inserting the 35S-eGFP fragment into the pCAMBIA2300 vector skeleton in our laboratory). The accuracy of the 10 recombinant vectors was ensured by primer-based PCR and sequencing analysis. Then, we cultured Nicotiana benthamiana seedlings at 28 °C day/25 °C night with a 14-hour photoperiod and 70% relative humidity for 1 month. When the third or fourth leaf was dark green and sufficiently thick, Agrobacterium tumefaciens harboring the recombinant vector was injected into the tobacco leaves (Poulsen, Palmgren & López-Marqués, 2016) to introduce the recombinant vector into tobacco leaf epidermal cells, and the plants were cultured in darkness at 25 °C for 24 h and then under light for 24 h. A square slice of 1 cm² was cut next to the injection area, placed on a slide with PBS buffer, and then gently covered with a cover glass. The position of the fluorescent protein in tobacco epidermal cells was observed by laser confocal microscopy (Olympus FV1200). The pCAMBIA2300-eGFP (CAMBIA) empty vector without the Gh-PRX gene was used as a control. The nuclei were stained with DAPI solution (Solarbio) as a control (Kapuscinski, 1995).

Genome-wide identification and conserved motif analysis of the PRX family in cotton

To identify the PRX genes of cotton, we used the redoxin domain, 1-cys-PRX_C domain and AhpC-TSA domain as references. The protein databases of four cotton species were searched by Hmmsearch (Bhaduri, Ravishankar & Sowdhamini, 2004). Based on the conserved active sites and the positions of conserved cysteinyl residues, 44 PRX genes were retrieved (Table 1), of which 15, 8, 8 and 13 PRXs were identified from G. hirsutum, G. raimondii, G. arboreum and G. barbadense, respectively. The PRXs were then divided into six subfamilies, named 2-CysPRX, 1-CysPRX, PRQ, PRXIIB, PRXIIE and PRXIIF, according to their conserved active sites, the positions of their conserved cysteinyl residues and the relevant literature (Dietz, 2011; Nelson et al., 2011) (Fig. 1).

We performed multiple sequence alignment of the 44 PRX proteins using ClustalW in MEGA7 and identified in detail the specific and conserved active sites of the PRX proteins (Adak & Begley, 2017). The results showed that the active sites of the PRX protein sequences were consistent with those (Fig. 1) by which the PRX family was classified in the PREX database (http://csb.wfu.edu/prex/). All PRX family proteins have conserved P, T, C, S and W/F amino acid sites, and the fifth conserved amino acid in the PRXQ subfamily protein sequence is F (Phe). In addition, we refer to the position and number of conserved cysteine residues in the classification of subfamilies (Dietz, 2011). 2-CysPRX has 2 Cys residues: Cys_P (FFYPLDFTFVCPTEI) is generally located in the N-terminal part of the protein at approximately amino acid (aa) position 118, and Cys_R (EVCP) is located in the C-terminal part of the protein at approximately aa position 240 (Fig. 1). The two catalytic Cys residues in PRXQ are separated by 4 amino acids (Fig. 1) (Kong et al., 2000). The two catalytic residues of PRXII, Cys_P and Cys_R, are separated by 24 amino acids. The subclassification of PRXII proteins has always been based on their subcellular localization (Dietz, 2011), In addition, we compared the known PRXII proteins in Arabidopsis, rice and cotton and analyzed the differences between their conserved sequences, which facilitates a more accurate classification of PRXs (Fig. S1). PRXIIE proteins are localized in the chloroplast (Mhamdi & Van Breusegem 2018; Romero-Puertas et al., 2007) and share the conserved sequences (FGLPGAYTGVCSQ….C.S…W). PRXIIF proteins are localized to the mitochondria (Horling, König & Dietz, 2002), and display a conserved motif (Fig. S1). Compared with the PRXIIF protein sequence in Arabidopsis, the S (serine) conserved was mutated to A (alanine) in cotton. The second cysteinyl (Cys) in the protein sequences of some genes originally predicted to belong to the PRXIIB subgroup (GbPRX6-A, GhPRX3-A, GhPRX11-D, GrPRX1 and GaPRX1) was mutated to V (Valine) (Fig. 1) (Table S5), resulting in these genes having a cysteine residue only at N-terminal amino acid position 50 and being classified as 1-CysPRX. Typical 1-CysPRX proteins do not have a second Cys (Fig. 1) (Dietz, 2011) and only have Cys_P followed by a peroxide domain.

The MW and pI of the PRXs were calculated using the Compute pI/Mw tool on the ExPASy (https://web.expasy.org/compute_pi/) website. The encoded protein lengths of the 44 PRX genes varied from 100 to 300 amino acids, the predicted MWs were 13-30 kilodaltons (kDa), and the pI values were approximately 4.0–10.0. Subcellular localization analysis predicted that 17 of the 44 proteins localized to the chloroplasts (Table S4); thus, these members of the PRX family may be involved in regulating the redox status of chloroplast proteins.

Phylogenetic and gene structure analysis of PRXs in cotton

We constructed a phylogenetic tree using 74 PRXs from eight dicot species (G. raimondii, G. arboreum, G. hirsutum, G. barbadense, Arabidopsis, rice, cacao, grapevine). The results show that the cotton PRX genes are divided into six subfamilies: 2-CysPRX, PRXQ, PRXIIE, PRXIIB, PRXIIF and 1-CysPRX (Fig. 2). Moreover, the cotton PRX genes were clustered within every clade, indicating that these genes had undergone specific expansions in all six subfamilies. A conserved domain analysis helped to clarify the functional evolution and evolutionary relationships of cotton PRXs. A total of 10 conserved motifs were identified in the PRX proteins of the 4 cotton species. The motif logo of these conserved motifs is shown in Fig. S2, and the flags for the 10 conserved motifs are shown in Fig. 3C. The number of conserved motifs in each PRX protein ranged from 4 to 9. There were four motif fragments containing conserved active structural sites (motif 1, motif 4, motif 5 and motif 6) (Fig. 3C), and the details can be seen in Fig. S3. The conserved structural sites in the PRXIIB, PRXIIE, PRXIIF, 1-CysPRX and PRXQ subfamilies are located in motif 1, while the conserved structural sites in 2-CysPRX and 1-CysPRX are composed of motif 4, motif 5 and motif 6. Compared with the motifs containing conserved active sites, the seven other motifs (2, 3, 7, 8, 9, 10) play important roles in the functional differentiation of PRX family genes.

Exon-intron structural differences play a key role in the evolution of polygene families. We found that the number of introns in the 44 PRX genes ranged from 0 to 6. PRXs in the same subfamily showed a similar exon-intron structure (Fig. 3B). The PRXIIB subfamily contained five intron sequences of different lengths. PRXIIE has no introns, and PRXIIF has 4 introns. PRXQ contained three introns, and most of the 2-CysPRX and 1-CysPRX subfamily genes also contained the same number of introns. However, there was an intron sequence before the translation initiation site (TIS) in GbPRX5-A and GbPRX14-D. Introns before the TIS are integral parts of genes and play important roles in regulating gene expression (Li, Li & Zhang, 2013).

Chromosome mapping and duplication analysis

The PRX genes were physically mapped to the chromosomes of G. hirsutum and G. barbadense. We found that the mapping of PRXs to the chromosomes of these two cotton species was similar (Figs. 4A and 4B), which shows that the genes of this family were conserved during evolution. However, the homologs of three of these PRX genes (GbPRX6-A, GhPRX11-D and GbPRX8-D), which are marked in black in the figure, may have been lost during evolution. There are six pairs of homologous genes in the G. hirsutum PRX family, which are distributed on chromosomes Gh_A01, Gh_A04, Gh_A08, Gh_A10, Gh_D01, Gh_D04, Gh_D08 and Gh_D10. Among the 15 Gh-PRX genes, GhPRX3-A was localized to scaffolds, and its exact position has not yet been determined. There are five pairs of PRX homologous genes in G. barbadense, which were mapped to chromosomes Gb_A01, Gb_A08, Gb_A10, Gb_D01, Gb_D08 and Gb_D10. Compared with those of G. hirsutum, the PRXs on Gb_A04 of G. barbadense were lost; therefore, one pair of homologous genes was missing. The comparative genomic analysis of PRXs between G. hirsutum and G. barbadense showed that GbPRX6-A on Gb_A05 originated from a translocation of GhPRX11-D on Gh_D05 (Fig. 4B and Fig. S4). Based on the physical chromosome mapping analysis and comparison of PRXs between the G. arboreum and G. raimondii genomes, we found that only the location of GaPRX5 and GrPRX6 on Chr02 was relatively conserved in these two cotton species (Figs. 4C–4D and Fig. S4).

We identified eight pairs of duplicates among the 15 PRXs of G. hirsutum (GhPRX8-D:GhPRX1-A, GhPRX9-D:GhPRX2-A, GhPRX10-D:GhPRX12-D, GhPRX10-D:GhPRX4-A, GhPRX12-D:GhPRX4-A, GhPRX5-A:GhPRX13-D, GhPRX15-D:GhPRX7-A, and GhPRX14-D:GhPRX6-A) (Fig. 4A), and these duplicates are located on six pairs of chromosomes (Gh_A01:Gh_D01, Gh_A04:Gh_D04, Gh_D05:Gh_D08, Gh_A08:Gh_D05, Gh_A08:Gh_D08, and Gh_A10:Gh_D10). With the exception of GhPRX10-D:GhPRX12-A and GhPRX10-D:GhPRX4-A, which underwent segmental duplication, these duplicates were produced during the process of genomic polyploidization. Compared with those in G. hirsutum, the duplicated genes on the G. barbadense Gb_A04 and Gb_D04 chromosomes were altered (Fig. 4B). GbPRX8-D on the Gb_D04 chromosome was conserved, although the corresponding homologous gene on Gb_A04 was lost during evolution. The other duplicated genes were the same as those in G. hirsutum and remained conserved. The PRXs of the diploid G. arboreum genome are very similar to the PRXs of the tetraploid G. hirsutum and G. barbadense D subgenome, only GaPRX5 was lost during evolution, and segmental duplication genes on Chr05 and Chr08 chromosomes were highly homologous (Fig. 4D and Fig. S4). However, in G. raimondii, these segmentally duplicated genes were physically located on the Gr_chr04 and Gr_chr09 chromosomes (Fig. 4C). According to a KaKs calculation of the duplicated gene pairs in G. raimondii and G. arboreum (Table S1), the divergence time of the duplicated genes (GrPRX3:GrPRX7) on the G. raimondii chromosomes Gr_chr04 and Gr_chr09 was 24.81 MYA, while the divergence time of the duplicated genes (GaPRX3:GaPRX7) on G. arboreum Ga_chr05 and Ga_chr08 was 23.98 MYA. The divergence time of these genes on Chr05 and Chr08 in he divergence time of G. hirsutum and G. barbadense was 23.2MYA. A comparative analysis of four cotton PRXs also showed that these segmental duplication genes were highly homologous, and the physical location of the PRXs in G. arboreum highly overlapped with that on the D subgenomes of G. hirsutum and G. barbadense.

Analysis of cis-acting elements in promoter regions

In plants, cis-acting elements play an important role in gene expression regulation. We found many cis-acting elements in the promoters of the PRX genes in the four cotton species (Table S2) and identified three main types of cis-acting elements.

The first type of cis-acting element was responsive to abiotic stress (Fig. 5A). There were nine types of elements (ARE, DRE1, LTR, MBS, MYC, STRE, TC-rich repeats, WRE3, and WUN-motif) in this category, among which the most abundant cis-acting elements were related to drought stress. The elements related to drought stress were MBS and MYC, and MYC was more frequent. GhPRX5-A, GhPRX6-A, GhPRX7-A, GhPRX9-D, GhPRX13-D, GbPRX1-A, GbPRX2-A, GbPRX5-A, GbPRX8-D, GbPRX9-D, GrPRX2, GrPRX8 and GrPRX5 contained more than five drought-related elements among the four cotton species. Second in terms of presence, the cis-acting elements involved in the response to osmotic stress and defense were STRE- and TC-rich. In addition, there were elements that responded to anaerobic conditions, hypothermia and wounding.

The second type of cis-acting element was responsive to plant hormones (Fig. 5B). There were ten types of elements (ABRE, AuxRR-core, CGTCA-motif, ERE, GaRE-motif, P-box, TATC-box, TCA-element, TGaCG-motif, and TGa-element) that responded to abscisic acid, auxin, ethylene, gibberellin and SA. Among them, ABRE (ABA) and ERE were the most abundant in the cotton PRX gene promoters and were found in 31 and 36 PRX genes, respectively. We also found three types of elements (GaRE-motif, P-box and TATC-box) involved in the gibberellin response in the PRX promoters. Among them, there were 2 GaRE-motifs, 3 TATC-boxs and 14 P-boxs in 44 promoters. In addition, auxin response elements (AuxRR-core and TGA-element), methyl jasmonate response elements (MeJA and TGACG-motif), and SA response elements (TCA-element) were found in many PRX promoters.

The third type of cis-acting element was related to plant growth and development (Fig. 5C). There were nine such types of elements (AACA-motif, AT-rich element, CAT-box, circadian, GCN4-motif, HD-Zip 1, MBSI, MYB, and O2-site). Among them, only one circadian element was involved in circadian regulation, and this element was found in the promoters of only two PRX genes (GhPRX10-D and GhPRX11-D) in G. hirsutum. PRXs may be less sensitive to circadian rhythms. The largest number of cis-acting elements were MYB protein binding sequence sites, and two types of MYB protein sequence binding site elements (MBSI and MYB) were identified. In addition, some cis-acting elements (GCN4-motif and AACA-motif) involved in endosperm expression, cis-acting regulatory elements (CAT-box) involved in meristem expression, and response elements involved in palisade mesophyll cell differentiation (HD-Zip 1) were also found, and some light-related cis-acting elements were also identified (Fig. S5).

Expression profiles of PRXs in G. hirsutum under stress

To investigate the regulation of PRX expression by cis-acting promoter elements, we analyzed the cis-acting elements that respond to drought (MBS and MYC), osmotic stress (STRE and TC-rich) and SA (TCA-element) in the promoter regions of Gh-PRX genes. The physical positions of these five cis-acting elements in the promoter regions of the Gh-PRX genes were plotted (Fig. 6). An element (MYC) involved in drought induction was present in all the Gh-PRX promoter regions, and the effects of drought on the expression levels of these genes were predicted to be strong.

We used PEG-6000, NaCl and SA to treat hydroponic G. hirsutum seedlings. Gene samples were extracted from G. hirsutum for qRT-PCR verification (Figs. 7 and 8). The expression level of genes in each time period was normalized to that at the 0 h time point. The results showed that nine PRXs were upregulated in plant roots under salt stress, including GhPRX1-A, GhPRX2-A, GhPRX5-A, GhPRX10-D, GhPRX11-D, GhPRX12-D, GhPRX13-D, GhPRX14-D and GhPRX15-D (Figs. 7A, 7B, 7E, 7J, 7K, 7L, 7M, 7N and 7O). The promoter regions of GhPRX11-D, GhPRX14-D and GhPRX15-D were rich in TC-rich elements, while GhPRX5-A, GhPRX10-D, GhPRX12-D and GhPRX13-D were rich in STRE elements. Some genes were repressed under stress for 1 h, such as GhPRX4-A, GhPRX6-A, GhPRX8-A and GhPRX9-D (Figs. 7D, 7F, 7H and 7I), in which GhPRX4-A, GhPRX6-A and GhPRX9-D were rich in STRE elements. However, in salt-treated leaves, the expression of most genes was inhibited, including GhPRX2-A, GhPRX4-A, GhPRX5-A, GhPRX6-A, GhPRX8-D, GhPRX9-D, GhPRX10-D, GhPRX12-D and GhPRX13-D (Figs. 8B, 8D, 8E, 8F, 8H, 8I, 8J, 8L and 8M), including the seven genes rich in STRE elements. Five genes were slightly upregulated under salt treatment, including GhPRX1-A, GhPRX3-A, GhPRX11-D, GhPRX14-D and GhPRX15-D (Figs. 8A, 8C, 8K, 8N and 8O). Among them, GhPRX11-D, GhPRX14-D and GhPRX15-D are rich in TC-rich elements. Although GhPRX14-D was slightly lower than the control at 0 h, it was higher than the expression of the control under stress. Under SA stress, five PRXs were downregulated in roots (GhPRX1-A, GhPRX2-A, GhPRX4-A, GhPRX6-A and GhPRX8-D) (Figs. 7A, 7B, 7D, 7F and 7H), seven PRXs were downregulated in leaves (GhPRX2-A, GhPRX4-A, GhPRX5-A, GhPRX6-A, GhPRX8-D, GhPRX10-D and GhPRX12-D) (Figs. 8B, 8D, 8E, 8F, 8H, 8J and 8L), and no significant upward trend was observed in all PRXs. Only GhPRX8-D in the root and GhPRX13-D in the leaf first decreased and then increased. GhPRX2-A, GhPRX4-A, GhPRX6-A and GhPRX8-D were rich in TCA-element cis-acting elements, and four genes were suppressed and downregulated in the roots and leaves.

MYC/MBS cis-acting elements exist in all Gh-PRX genes. Under drought stress simulated by PEG-6000, the expression of most Gh-PRX genes in roots first increased and then decreased (Fig. 7). Especially at 1 h, the expression levels of 10 of the 15 PRX genes increased significantly and then decreased rapidly (GhPRX1-A, GhPRX2-A, GhPRX4-A, GhPRX5-A, GhPRX6-A, GhPRX8-A, GhPRX10-D, GhPRX11-D, GhPRX13-D and GhPRX15-D) (Figs. 7A, 7B, 7D, 7E, 7F, 7H, 7J, 7K, 7M and 7O), which indicated that the PRXs in roots were highly sensitive to drought stress and induced. However, in leaves, most PRXs showed a downregulation trend, such as GhPRX2-A, GhPRX4-A, GhPRX5-A, GhPRX6-A, GhPRX8-D, GhPRX10-D and GhPRX12-D. Only GhPRX9-D showed the same expression as roots, first upregulated and then downregulated (Fig. 8I).

Expression and AS profiles of PRXs in different tissues of G. hirsutum

It is important to analyze the expression of genes in plant tissues and organs to understand gene function. We analyzed the available transcriptome data of various tissues and organs of G. hirsutum TM-1 (Zhang et al., 2015), including the roots, stems, leaves, cotyledons, petals, stamens, ovules, seeds, and fibers at 5, 10, 15 and 20 DPA. A heatmap was used to show the expression of the PRX genes (Fig. 9A). GhPRX3-A, GhPRX5-A, GhPRX10-D, GhPRX11-D and GhPRX13-D were expressed at low levels or not expressed in all tissues. GhPRX3 and GhPRX12 were highly expressed in all tissues, and the FPKM value of GhPRX14 expression level in roots and leaves was approximately 50–60.

Understanding and verifying the AS events of the G. hirsutum PRX family is an important step in the study of gene functional differentiation. Five main types of AS events were used for further analysis: exon skipping (ES), intron retention (IR), 5′ or 3′ alternative splice sites (A5SS or A3SS), alternative first exon (AFE), and alternative last exon (ALE). Because change in the CDS can modify protein sequence and function, only AS events located in the CDSs of the genes were used for further analysis. Using the full-length primers designed for the amplification of PRX gene cDNA, RT-PCR amplification was performed using the leaf and root cDNA of G. hirsutum sGK9708 as a template. The results showed that the products obtained for GhPRX14-D included multiple fragments (Figs. 9B, 9C and 9E); moreover, s everal bands with different levels of brightness and unequal size were present. Thus, we speculate that GhPRX14-D may have multiple transcripts. The PCR products of GhPRX14-D were recovered, purified and ligated into a cloning vector for replication and sequencing. Analysis of the obtained sequences showed two different transcripts of GhPRX14-D in the leaves of G. hirsutum (Table 2) and three in the roots. GhPRX14-D uses three AS methods, IR, ES and A3SS. GhPRX14-Leaf-AS2 and GhPRX14-Root-AS2 use one nonstandard AS site: 5′-AG. AA-3′ (3 clones and 2 clones), which retains 15 bases at the 3′ end of the first intron (Fig. 9D). GhPRX14-Root-AS3 uses two nonstandard AS sites: 5′-AG. AT-3′ (2 clones) and 5′-GG. TT-3′ (2 clones), which retain 14 bases at the 3′ end of the fourth intron. The second, third, and fourth exons were skipped (Fig. 9D).

Subcellular localization of PRX proteins

Determining the distribution of proteins in cells is an important step to verify protein function. Using the TargetP 1.1 Server (http://www.cbs.dtu.dk/services/TargetP/) website to predict and analyze proteins, we found that the PRX proteins were mainly located in chloroplasts and plasma membranes. However, the specificity of a few gene prediction results was too low to determine protein localization (Table S4).

To verify the predicted subcellular localizations of the PRX proteins, we selected 10 PRXs with high expression levels for subcellular localization analysis (5 pairs of homologous genes). We constructed a transient expression 35s-PRX-eGFP vector for each gene that was used to verify the experimental purpose. The results showed that the eGFP fluorescence signal in the control group was dispersed throughout the tobacco leaf cells, while the fluorescence signals of GhPRX4-A and GhPRX12-D, which are members of PRXIIB, were located in the cell membrane (Figs. 10D–10I), thus providing support for the accurate classification of PRXIIB subfamily members. The protein fluorescence signals of GhPRX7-A and GhPRX15-D, members of 1-CysPRX, were located on the cell membrane and nucleus (Figs. 10J–10M and 10N–10Q). GhPRX1-A and GhPRX8-D of PRXIIE, GhPRX6-A and GhPRX14-D of PRXIIF, GhPRX2-A and GhPRX9-D of 2-CysPRX, and their protein fluorescence signals were located in chloroplasts (Figs. 10R–10OO).