Genome-Wide Analysis of MDHAR Gene Family in Four Cotton Species Provides Insights into Fiber Development via Regulating AsA Redox Homeostasis

Monodehydroasorbate reductase (MDHAR) (EC1.6.5.4), a key enzyme in ascorbate-glutathione recycling, plays important roles in cell growth, plant development and physiological response to environmental stress via control of ascorbic acid (AsA)-mediated reduction/oxidation (redox) regulation. Until now, information regarding MDHAR function and regulatory mechanism in Gossypium have been limited. Herein, a genome-wide identification and comprehensive bioinformatic analysis of 36 MDHAR family genes in four Gossypium species, Gossypium arboreum, G. raimondii, G. hirsutum, and G. barbadense, were performed, indicating their close evolutionary relationship. Expression analysis of GhMDHARs in different cotton tissues and under abiotic stress and phytohormone treatment revealed diverse expression features. Fiber-specific expression analysis showed that GhMDHAR1A/D, 3A/D and 4A/D were preferentially expressed in fiber fast elongating stages to reach peak values in 15-DPA fibers, with corresponding coincident observances of MDHAR enzyme activity, AsA content and ascorbic acid/dehydroascorbic acid (AsA/DHA) ratio. Meanwhile, there was a close positive correlation between the increase of AsA content and AsA/DHA ratio catalyzed by MDHAR and fiber elongation development in different fiber-length cotton cultivars, suggesting the potential important function of MDHAR for fiber growth. Following H2O2 stimulation, GhMDHAR demonstrated immediate responses at the levels of mRNA, enzyme, the product of AsA and corresponding AsA/DHA value, and antioxidative activity. These results for the first time provide a comprehensive systemic analysis of the MDHAR gene family in plants and the four cotton species and demonstrate the contribution of MDHAR to fiber elongation development by controlling AsA-recycling-mediated cellular redox homeostasis.


Introduction
Reactive oxygen species (ROS) are important signaling molecules that play significant roles in plant development, cell growth and in response to environmental stresses and phytohormone signaling [1][2][3][4][5]. To counteract excessive ROS and the detrimental effects of cellular oxidation, plants have evolved corresponding defense mechanism for detoxification of ROS by controlling cellular reduction/oxidation (redox) balance. In this, the ascorbate-glutathione (Asc-GSH) cycle, also namely Asc recycling pathway, is regarded as the critical pathway that mediates the antioxidative molecule of Asc as a coenzyme appeared with an ethylene-induced expression and participated in fiber cell fast elongation development by modulating ROS homeostasis [4,48]. Moreover, the Asc de novo biosynthesis gene GhVTC1 is the key factor to promote cell elongation under control of ethylene [50], while the myo-inositol-1-phosphate synthase (MIPS) located in the MI-dependent D-glucuronate pathway participated in cell elongation by providing the precursor MI to contribute to Asc biosynthesis [51]. Ascorbate oxidase (AO) is a positive regulator to enhance tobacco cell elongation through generation of H 2 O 2 -mediated apoplast oxidation [16]. Together, these results suggest an intimate connection between Asc metabolic genes and cellular ROS balance in fibers. Currently, information regarding MDHAR genes and their functions, as well as the regulatory mechanism involved in fiber development in Gossypium, are largely unknown. In the present study, by a genome-wide identification in four cotton species of G. arboreum, G. raimondii, G. hirsutum and G. barbadense, a total of 36 MDHAR genes were obtained. A comprehensive systematic bioinformatic analysis, including chromosomal location, gene structure, motif distribution and evolutionary relationships, was performed indicating rapid evolution and close evolutionary relationships. G. hirsutum GhMDHARs that showed diverse expression patterns under abiotic stress, and phytohormone treatment demonstrated a close correlation with fiber growth by increasing MDHAR activity, AsA content and AsA/DHA ratio, possibly via involvement of AsA-mediated H 2 O 2 signaling pathway.

Identification and Chromosome Location of MDHAR Genes of Four Gossypium Species
By searching the genome of G. arboreum, G. raimondii, G. hirsutum and G. barbadense using arabidopsis MDHAR proteins as queries with the BlastP procedure, a total of 36 MD-HAR genes were identified, with six in diploid cotton (G. arboreum and G. raimondii), and twelve in allotetraploid cotton (G. barbadense and G. hirsutum) (Supplementary Table S1). The allotetraploid cotton species have twice the MDHAR members in comparison to the diploid cotton species. The MDHARs possess open reading frame (ORF) lengths ranging from 1305 to 1677 base pair (bp), the encoded proteins have 434 to 558 amino acids with theoretical molecular weight (MW) from 46.83 to 62.25 kDa and theoretical pI from 5.9 to 9.06, respectively (Supplementary Table S1).

Phylogenetic and Evolutionary Analyses of MDHAR Members
The MDHAR sequences from four cotton species of G. arboreum, G. raimondii, G. barbadense and G. hirsutum, and from another eight plants including Glycine max, Oryza sativa, Sorghum bicolor, Zea mays, Arabidopsis thaliana, Theobroma cacao, Vitis vinifera and Cicer arietinum, were used for construction of a neighbor-joining (N-J) tree, resulting in three typical groups (Group A-C). The orthologous MDHARs in cotton and cacao indicated a close relationship, which might be the reason that they are in the same family of Malvaceae. Cotton MDHAR homologs were classified into the same branches, with MDHAR1, 3, and 6 in group A, MDHAR 2 and 5 in group B, and MDHAR 4 in group C (Figure 2), demonstrating again the close relationship of gene duplication of MDHARs in the same groups.
The evolutionary relationship could be assessed by genetic synteny across related plant species. The syntenic analysis of MDHAR genes of the four cotton species was executed and visualized using Circos software (Supplementary Figure S1). Two paralogs in G. arboreum, one paralog in G. raimondii, nine paralogs in G. barbadense and nine paralogs in G. hirsutum had syntenic relationships. Meanwhile, the number of 47, 27 and 15 had syntenic relationships in Group A, Group B and Group C, respectively (Figure 2), showing that Group A had many syntenic genes compared with the other two groups, indicating the close evolutionary relationships of the MDHAR genes among the four cotton species. The evolutionary rate of MDHAR duplicated gene pairs was evaluated by the Tajima relative rate test, which indicated that MDHAR1/6 and MDHAR2/5 gene pairs

Phylogenetic and Evolutionary Analyses of MDHAR Members
The MDHAR sequences from four cotton species of G. arboreum, G. raimondii, G. barbadense and G. hirsutum, and from another eight plants including Glycine max, Oryza sativa, Sorghum bicolor, Zea mays, Arabidopsis thaliana, Theobroma cacao, Vitis vinifera and Cicer arietinum, were used for construction of a neighbor-joining (N-J) tree, resulting in three typical groups (Group A-C). The orthologous MDHARs in cotton and cacao indicated a close relationship, which might be the reason that they are in the same family of Malvaceae. Cotton MDHAR homologs were classified into the same branches, with MDHAR1, 3, and 6 in group A, MDHAR 2 and 5 in group B, and MDHAR 4 in group C ( Figure 2), demonstrating again the close relationship of gene duplication of MDHARs in the same groups.
The evolutionary relationship could be assessed by genetic synteny across related plant species. The syntenic analysis of MDHAR genes of the four cotton species was executed and visualized using Circos software (Supplementary Figure S1). Two paralogs in G. arboreum, one paralog in G. raimondii, nine paralogs in G. barbadense and nine paralogs in G. hirsutum had syntenic relationships. Meanwhile, the number of 47, 27 and 15 had syntenic relationships in Group A, Group B and Group C, respectively (Figure 2), showing that Group A had many syntenic genes compared with the other two groups, indicating the close evolutionary relationships of the MDHAR genes among the four cotton species. The evolutionary rate of MDHAR duplicated gene pairs was evaluated by the Tajima relative rate test, which indicated that MDHAR1/6 and MDHAR2/5 gene pairs showed accelerated evolutionary rate in the four Gossypium species, and the MDHAR3/6 gene indicated accelerated evolutionary rate only in G. hirsutum (Supplementary Table S2), implying the potential functional divergence of these MDHAR gene pairs. Evolutionary selection analysis by nonsynonymous substitution rate/synonymous substitution rate (Ka/Ks) indicated that the ratio of all duplicated MDHAR gene pairs showed a value of <1, indicating all these MDHAR gene pairs were purifying selection and suggesting their rapid evolution (Supplementary Table S3).
showed accelerated evolutionary rate in the four Gossypium species, and the MDHAR3/6 gene indicated accelerated evolutionary rate only in G. hirsutum (Supplementary Table  S2), implying the potential functional divergence of these MDHAR gene pairs. Evolutionary selection analysis by nonsynonymous substitution rate/synonymous substitution rate (Ka/Ks) indicated that the ratio of all duplicated MDHAR gene pairs showed a value of <1, indicating all these MDHAR gene pairs were purifying selection and suggesting their rapid evolution (Supplementary Table S3).

Analyses of Gene Structure, Conserved Motif, and Promoter Cis-Acting Element of MDHAR Genes
To obtain more insights into the evolutionary relationship and potential functional diversity of MDHAR genes, exon/intron organization and conserved motif analyses were performed by Gene Structure Display Server (GSDS) and Multiple Em for Motif Elicitation) (MEME) programs respectively. The MDHAR genes could be divided into different branches ( Figure 3A). The same homologous MDHAR genes in different cotton species were classified into the same branch and indicated a highly similar exon/intron arrangement, with the exception of MDHAR6 that contained only one individual exon without intron ( Figure 3B), suggesting conservation of the evolutionary relationship. All MDHARs presented high sequence similarity and contained a typical flavin adenine dinucleotide/nicotinamide adenine dinucleotide (FAD/NAD)-binding domain (Supplementary Figure S2), implying a conserved function in catalysis. As a result of nine discovered conserved motifs, nine, eight and seven conserved motifs were presented in MDHAR1, 3, and 6, MDHAR2 and 5, and MDHAR4, respectively, showing the corresponding predicted subcellular locations in the cytosol, membrane-bound peroxisomes, and chloroplast/mitochondria ( Figure 3C). The detailed amino acid composition of the conserved motifs is presented in Supplementary Figure S3, suggesting the functional diversity of the different MDHARs. Cis-acting element analysis showed that the distributions in different MDHAR promoter regions are largely diverse, with the exception of the highly consistent composition of MDHAR promoters in A-subgenome but not in D-subgenome in the two allotetraploid species G. hirsutum and G. barbadense. Besides the core elements, regulatory elements related to development and responsiveness to stress and phytohormone were also discovered (Supplementary Figure S4) with detailed information of these elements in Supplementary Table S4, suggesting their potential function for plant growth and regulation by these factors.  Figure 4B). These results suggest that these GhMDHAR genes might play important function for fiber development.
As MDHARs have a main function as catalytic enzymes, MDHAR activity was measured during fiber development stages. MDHAR activity indicated significant accumulation in 5-21 DPA fibers and reached a peak value at 15 DPA, compared to that in 0-DPA ovules ( Figure 4C). Further, AsA content was measured to validate the key role of MDHAR in AsA regeneration. The results showed that AsA content maintained a high level in 10-21 DPA fibers (the fiber fast elongation and secondary wall deposition stages) and reached the highest level in 15-DPA fibers. The AsA/DHA ratio showed a similar tendency with significant increase at 10-21 DPA with the peak in 15-DPA fibers ( Figure  4D). These results demonstrated that a higher MDHAR activity is associated with fiber growth by catalyzing AsA recycling and maintaining the AsA-mediated cellular redox status.

Tissue-Specific and Fiber-Developmental Expression Analysis of G. hirsutum GhMDHARs
The upland cotton G. hirsutum provides the majority of the fibers for textile industry. We analyzed the expression pattern of MDHARs during fiber development and did tissuespecific and fiber-developmental expression analyses in G. hirsutum GhMDHARs. The results indicated that GhMDHAR1A/D and GhMDHAR3A/D were preferentially expressed in all detected tissues including roots, stems, leaves, petals, stamens, pistils, ovules and fibers, with the most significant accumulations in fiber development, suggesting their potential important roles in tissue development especially in fiber growth ( Figure 4A). GhMDHAR2A/D indicated moderate expression in all tissues, with the exception of their expressions in fibers and of GhMDHAR2A expression in flower tissues, implying a possible functional Plants 2021, 10, 227 7 of 20 difference in flower development. GhMDHAR5A/D also showed moderate expression with an especially high enrichment in flower tissues, presenting their potential function in flower development. GhMDHAR4A/D and GhMDHAR6A/D displayed consistent expression with a moderate and low expression, respectively. Meanwhile, in the other three cotton species, G. barbadense, G. arboreum and G. raimondii, the expression of MDHARs showed highly similar levels with that of G. hirsutum GhMDHARs (Supplementary Figure S5), implying the diverse and important roles of MDHAR expressions for both fiber growth and tissue development.
To further investigate the possible role of GhMDHARs in fiber development, the expressions of GhMDHAR genes during different fiber development stages were detected by quantitative real-time polymerase chain reaction (qRT-PCR). The results indicated that, during different fiber developmental stages, GhMDHAR3A/D showed the most significant accumulations in 5-21 days post anthesis (DPA) fibers (the fast elongation stages), with the peak value at 15 DPA. GhMDHAR1D also presented moderate upregulated expression levels in 5-21 DPA fibers. GhMDHAR1A and GhMDHAR4A/D were mainly enriched in 15-DPA fibers ( Figure 4B). These results suggest that these GhMDHAR genes might play important function for fiber development.
As MDHARs have a main function as catalytic enzymes, MDHAR activity was measured during fiber development stages. MDHAR activity indicated significant accumulation in 5-21 DPA fibers and reached a peak value at 15 DPA, compared to that in 0-DPA ovules ( Figure 4C). Further, AsA content was measured to validate the key role of MDHAR in AsA regeneration. The results showed that AsA content maintained a high level in 10-21 DPA fibers (the fiber fast elongation and secondary wall deposition stages) and reached the highest level in 15-DPA fibers. The AsA/DHA ratio showed a similar tendency with significant increase at 10-21 DPA with the peak in 15-DPA fibers ( Figure 4D). These results demonstrated that a higher MDHAR activity is associated with fiber growth by catalyzing AsA recycling and maintaining the AsA-mediated cellular redox status.

Expression Analysis of GhMDHARs under Abiotic and Phytohormone Treatments
Expression of MDHAR is regulated by environmental stress, especially drought and salt [25,52], and AsA biosynthesis and metabolism are controlled by phytohormone indole-3-acetic acid (IAA) and ethylene (ETH) [4,16,50]. Thus, to investigate the GhMDHAR expression features under treatments of abiotic stressors, polyethylene glycol (PEG) and NaCl and plant hormones of IAA and ETH treatments were used, and expression abundance of GhMDHAR genes was detected by qRT-PCR, using 0-h expression level of untreated samples as control. The GhMDHAR genes exhibited diverse expression patterns under different treatment conditions. GhMDHAR2A/D, GhMDHAR5A/D, and GhMDHAR6A/D indicated significant accumulated expression under these four treatments, at least in one treated time point, with the exceptions of GhMDHAR5A/D that showed moderate enrichment and of GhMDHAR6A/D that also showed moderate regulation under IAA treatment ( Figure 5A). GhMDHAR3A/D displayed a medium accumulation upon PEG, NaCl, and IAA treatment, but was enriched highly under the 24-h ETH treatment. GhMDHAR4A/D illustrated significant increased expression under the treatments of 24-h IAA and 12-h ETH without obvious difference following PEG and NaCl stresses ( Figure 5A). The MDHAR family demonstrated significant increased enzyme activity under these four conditions ( Figure 5B). The contents of AsA, together with the AsA/DHA ratio, indicated a highly consistency with MDHAR activity, showing significant accumulations under conditions of PEG, NaCl, IAA, and ETH treatments ( Figure 5C). These results imply that, GhMDHARs

Analysis of GhMDHAR Involving in Fiber Elongation Development
Our previous studies reported AsA biosynthesis and metabolism as positive regulators involved in fiber growth and cell elongation [16,50]. Therefore, to investigate the GhMDHAR potential function for fiber development, cotton cultivars that have short-(G2), medium-(LO5R59 and X142), and long-fiber (J12 and XH36) lengths were utilized to analyze the connection between MDHAR expression and fiber growth. As presented in Figure 6A,B, these cotton cultivars showed significant differences in fiber length. Meanwhile, MDHAR enzyme activity significantly increased and the AsA content and the AsA/DHA ratio were also significantly increased in the medium and long-fiber cultivars, compared to that in the short-fiber G2 ( Figure 6C,D). These results provide a close link between high MDHAR expression and fiber development, suggesting a possible important role for MDHAR in fiber cell elongation via regulation of AsA generation and control of AsA-mediated redox status.

Analysis of AsA Antioxidative Enzymes under H2O2 Stimulation
With regard to the key role of ROS in plant cell and fiber cell elongation [4,16,48,53,54], and to the important function of the reductive molecule AsA and the corresponding antioxidative enzymes in redundant ROS elimination during the fiber development process [4], expressions of GhMDHAR and the antioxidative enzymes were measured following H2O2 treatment. The results indicated that, at the transcription level, the expressions of GhMDHAR2A/D and GhMDHAR4A/D were significantly induced after 6-h H2O2 stimulation (Figure 7A), implying a potential crucial function in response to ROS. Contrary to the expression level of untreated material, the H2O2-treated materials showed a significant promotion of MDHAR activity after 24-h supplement, together with high levels of AsA content and AsA/DHA ratio after 6-h addition of H2O2, with a continuous increased level thereafter (Figure 7B,C). Other antioxidative enzyme activities were measured. Superoxide dismutase (SOD) showed a prompt increase after 6-h H2O2 treatment and APX and catalase (CAT) indicated significant promotions after 12-h H2O2 The values of the short-length cultivar G2 were selected as control for statistics difference analysis by one-way ANOVA. Asterisks indicate the significant differences, * p < 0.05, ** p < 0.01, *** p < 0.001.

Analysis of AsA Antioxidative Enzymes under H 2 O 2 Stimulation
With regard to the key role of ROS in plant cell and fiber cell elongation [4,16,48,53,54], and to the important function of the reductive molecule AsA and the corresponding antioxidative enzymes in redundant ROS elimination during the fiber development process [4], expressions of GhMDHAR and the antioxidative enzymes were measured following H 2 O 2 treatment. The results indicated that, at the transcription level, the expressions of GhMD-HAR2A/D and GhMDHAR4A/D were significantly induced after 6-h H 2 O 2 stimulation ( Figure 7A), implying a potential crucial function in response to ROS. Contrary to the expression level of untreated material, the H 2 O 2 -treated materials showed a significant promotion of MDHAR activity after 24-h supplement, together with high levels of AsA content and AsA/DHA ratio after 6-h addition of H 2 O 2 , with a continuous increased level thereafter (Figure 7B,C). Other antioxidative enzyme activities were measured. Superoxide dismutase (SOD) showed a prompt increase after 6-h H 2 O 2 treatment and APX and catalase (CAT) indicated significant promotions after 12-h H 2 O 2 stimulation ( Figure 7D). These results demonstrate potential GhMDHAR candidates to participate in H 2 O 2 response, and indicate that MDHAR, its catalyzed AsA regeneration and the successive cellular redox situation, combined with the significantly increased activity of the antioxidative enzymes, might play important role in detoxification of ROS.
Plants 2021, 10, x FOR PEER REVIEW 12 of 21 stimulation ( Figure 7D). These results demonstrate potential GhMDHAR candidates to participate in H2O2 response, and indicate that MDHAR, its catalyzed AsA regeneration and the successive cellular redox situation, combined with the significantly increased activity of the antioxidative enzymes, might play important role in detoxification of ROS.

Discussion
MDHAR is an essential enzyme for the regeneration of AsA that performs multiple roles in cell metabolism and physiological response as antioxidative small molecules [55,56], and thus plays important roles in plant development and cell growth [13,15,16] via a nicotinamide adenine dinucleotide phosphate (NADPH)-dependent redox reaction. A plethora of MDHAR genes have been identified in plants including pea, tomato, spinach, cucumber, soybean and potato [32][33][34][35][36][37]. Despite three PpMDHAR in Physcomitrella patens and six AtMDHARs in A. thaliana being reported, there are rare studies for the MDAHR family genes at the genome-wide level with a comprehensive systematic analysis. In this work, by a genome-wide analysis in four cotton species of G. arboreum, G. raimondii, G. hirsutum and G. barbadense, a total of 36 MDHAR members (Supplementary Table S1) were identified that located onto different chromosomes, with six members in diploid cotton G. arboreum and G. raimondii, and 12 members in allotetraploid cotton G.  barbadense and G. hirsutum, respectively (Figure 1).
Plant MDHARs are typical FAD monomeric enzymes that catalyze redox reactions using FADH2 as substrate to reduce MDHA to AsA via a GSH-coupled pathway [6,9].

Discussion
MDHAR is an essential enzyme for the regeneration of AsA that performs multiple roles in cell metabolism and physiological response as antioxidative small molecules [55,56], and thus plays important roles in plant development and cell growth [13,15,16] via a nicotinamide adenine dinucleotide phosphate (NADPH)-dependent redox reaction. A plethora of MDHAR genes have been identified in plants including pea, tomato, spinach, cucumber, soybean and potato [32][33][34][35][36][37]. Despite three PpMDHAR in Physcomitrella patens and six AtMDHARs in A. thaliana being reported, there are rare studies for the MDAHR family genes at the genome-wide level with a comprehensive systematic analysis. In this work, by a genome-wide analysis in four cotton species of G. arboreum, G. raimondii, G. hirsutum and G. barbadense, a total of 36 MDHAR members (Supplementary Table S1) were identified that located onto different chromosomes, with six members in diploid cotton G. arboreum and G. raimondii, and 12 members in allotetraploid cotton G. barbadense and G. hirsutum, respectively (Figure 1).
Plant MDHARs are typical FAD monomeric enzymes that catalyze redox reactions using FADH 2 as substrate to reduce MDHA to AsA via a GSH-coupled pathway [6,9]. Cotton MDHARs have a representative FAD/NAD binding domain (Supplementary Figure S2) and conserved motifs ( Figure 3C), which endow them with the function to catalyze the AsA recycling. Phylogenetic tree analysis constructed by a total of 70 MDHAR sequences from the four cotton species and eight other plant species resulted in three groups, implying a consistent function of the plant MDHARs in the same groups. The 36 cotton MDHARs were classified into different branches with the distinct gene structure organization and motif distribution in different branches (Figure 3), suggesting their potential diverse functions. It was established that the MDHARs that target to different cell compartments perform different functions such as cytosolic PpMDHARs for stress protection [39], peroxisomal PsMDHAR for high light intensity adaptation [57], and chloroplastic LeMD-HAR for photoinhibition alleviation [23]. Cotton MDHARs showed different predicted subcellular locations (Figure 3C), maintaining a high consistency with the groups of the phylogenetic tree (Figure 2), implying, again, possible multiple functions of the different targeted MDHARs.
There is abundant evidence that MDHAR genes are expressed correlatively with plant development. Induced expression of pea MDHAR performed important function for seed germination and seedling growth through regulating AsA content [58]. AtMD-HAR4 expression is essential for seedling growth through supply of carbon skeletons and energy, via protection of the seed storage oil hydrolysis process by controlling peroxisome ROS levels [29]. In this study, GhMDHARs displayed diverse expression features in different cotton tissues ( Figure 4A), suggesting their potential diverse functions in tissue development and plant growth. Considering that AsA is the crucial substance for plant cell enlargement and expansion [13,56], as the important AsA recycling enzyme, several GhMDHARs presented significantly promoted expressions during fiber cell development stages, with the highest level in the fast-elongating stage of 15 DPA, at both mRNA and enzyme activity levels ( Figure 4A,B), suggesting there exists a close correlation between MDHAR expression and fiber growth.
MDHAR has been widely studied, especially in the field of plants that resist environmental stress [59]. Eight cotton genes GhMDHAR2A/D, 3A/D, 5A/D and 6A/D indicated significant induced expressions, at both transcriptional and enzymatic levels, under the stresses of PEG and NaCl (Figure 5 A,B). P. patens PpMDHARs were differentially induced with diverse expression patterns under different stress conditions [39]. Avicennia marina AmMDHAR is a salt-inducible gene located in the chloroplast, and its overexpression in tobacco enhanced the tolerance of transgenic plants by increasing the activity of AsA-GSH antioxidative enzymes [52]. Both the mRNA and enzyme expressions of Acanthus ebracteatus AeMDHAR were significantly elevated in the overexpression transgenic rice plant lines that showed increased tolerance to salt stress at both germination and seedling stages [60]. The transgenic tobacco plants overexpressing arabidopsis AtMDHAR1 exhibited strengthened resistance to salt and PEG stresses, showing an increase of MDHAR activity and reductive AsA content coupled with a higher AsA/DHA ratio and lower ROS level as the key factors for the enhanced stress tolerance [25]. In this study, under PEG and NaCl conditions, an elevation of MDHAR enzyme activity and the promotion of corresponding AsA content and AsA/DHA ratio were observed ( Figure 5B,C), which displayed the positive connection between stress tolerance and AsA level increase by promotion of MDHAR expression, indicating consistency with reported studies. Phytohormones have long been established as crucial factors for plants to sense and resist the environmental stress [61]. Our results detected the rapid significant induced expressions of GhMDHAR2A/D and 4A/D under IAA treatment and of GhMDHAR2A/D, 3A/D, 4A/D, 5A/D and 6A/D under ETH stimulation, together with the coincident synchronous increases of AsA and AsA/DHA value, validating again the close correlation between phytohormone and stress perception in plants.
AsA plays multifunction roles for redundant ROS detoxification, cell expansion or enlargement, and plant development [13,15,56]. A higher AsA content that might be catalyzed by increased MDHAR activity in longer fibers was observed (Figure 6), which indicated a similar result of AsA accumulation in elongating fibers in cotton [50,62]. Our former studies discovered that AsA accumulation, biosynthesis and metabolism play a key role for fiber cell growth [4,50]. The cell wall is the crucial factor to decide cell length and mechanical strength and is affected by wall-localized hydroxyproline rich glycoproteins (HRGPs) [7]. AsA could influence cell wall alterations as a cofactor of HRGP to promote cell elongation [63][64][65][66], suggesting the possibility of MDHAR participating in cell elongation through the involvement of cell wall assembly indirectly by catalyzing AsA regeneration. ROS, as oxidative molecules, are decisive factors for both plant cell growth and fiber elongation development by affecting redox homeostasis of the cells [4,48,53,67]. Under H 2 O 2 treatment, the MDHAR expression, AsA content and AsA/DHA value, together with AsA-mediated antioxidative enzymes, were significantly elevated (Figure 7), indicating a close link between MDHAR and ROS signaling. Despite the known role of AsA as a reducing agent to respond to ROS, other oxidative molecules were generated by AsAmediated oxidation reactions, suggesting another possibility for MDHAR involvement in cell elongation via redox signaling pathways [4,16].
Plant hormones such as ethylene and IAA have been reported to control plant root hair, hypocotyl development, and fiber cell elongation [68][69][70][71][72][73][74][75][76][77]. The content of IAA significantly accumulated in the fiber initiation stage and performed important functions for fiber quality [74]. An ethylene signal was the crucial element to decide fiber growth [71,72,75]. Cotton GhMDHARs were significantly induced by IAA and ETH ( Figure 5A), coupled with the existing hormone responsive cis-elements (Supplementary Figure S4), showing their potential regulatory role as downstream signaling factors. As a cofactor of 1-aminocyclopropane-1-carboxylic acid oxidase (ACO), which is the key enzyme for ethylene biosynthesis, AsA may control plant cell development [78][79][80], suggesting again the correlation of MDHAR and ethylene signaling. With regard to the existing crosstalk between ROS and phytohormones, and to different subcellular distributions of MDHARs, further investigation into their synergies and functions in different cellular compartments are expected. In conclusion, our results for the first time provide comprehensive systemic analysis of the MDHAR gene family in plants and the four Gossypium species and indicate that MDHAR may be involved in fiber elongation development via controlling AsA-recycling-mediated cellular redox homeostasis.

Evolutionary Analysis and Promoter Cis-Element Identification
Multiple sequence alignment of the MDHAR proteins was performed using ClustalW and the phylogenetic tree of the MDHAR proteins from four cotton cultivars and another eight plant species was constructed by Molecular Evolutionary Genetics Analysis (MEGA 7.0, Institute for Genomics and Evolutionary Medicine, Temple University, Philadelphia, USA) with the parameters of neighbor-joining (N-J) method of 1000 bootstrap replicates, model of p-distance, and Gaps/Missing data treatment by Pairwise deletion. The exonintron distribution and conserved motif analyses were recognized using the online software of the gene structure display server (GSDS, http://gsds.cbi.pku.edu.cn/index.php) and Multiple EM for Motif Elicitation (MEME, http://meme-suite.org/tools/meme) respectively [88,89]. The Ka/Ks was calculated by Dnasp 5.0 software (Version 5.0, Universitat de Barcelona, Barcelona, Spain) [90]. The Ka/Ks ratios for the GhMDHAR genes were calculated to assess the selection pressure on duplicated genes, with Ka/Ks ratio >1, <1, or = 1 to indicate positive selection, negative selection and purify evolution, respectively [91]. The Tajima relative rate test were used to estimate the evolutionary rate by MEGA 7.0 software. The Circos software (Canada's Michael Smith Genome Sciences Center, Vancouver, BC, Canada) was used for collinear analysis and the collinear file was calculated by Multiple Collinearity Scan toolkit (MC Scan, Plant Genome Mapping Laboratory, University of Georgia, Georgia, USA) [92]. The TBtools software (State Key Laboratory for Conservation and Utilization of Subtropical Agro-Bioresources, Guangdong, China) was utilized to extract the 1500 bp promoter region upstream of start codon which was then submitted to the online software PlantCARE (http://bioinformatics.psb.ugent.be/webtools/plantcare/) for analyzing the cis-acting element composition [93].

Cotton Planting and Material Treatment
Upland cotton (G. hirsutum) varieties (X142, G2, LO5R59, J12) and Sea-island cotton (XH36) plants were cultivated in the experimental field of Shihezi University in Shihezi city, Xinjiang, China. The materials of X142 were collected for gene expression analysis, enzyme activity measurement, and treatments of stress and phytohormone. Cotton bolls of X142 were marked with a label on the day post anthesis (defined as 0 DPA) and the ovules and fibers were then collected in different developmental time points of 5, 10, 15, 18, 21 and 28 DPA for cryopreservation and successive use. Cotton seeds were sown in nutrient soil and planted in the greenhouse at 28 • C with a photoperiod of 16 h light/8 h dark. The leaves of three-week-old seedlings of X142 were treated with 20% polyethylene glycol (PEG), 100 mM sodium chloride (NaCl), 100 µm indole-3-acetic acid (IAA), 100 µm ethephon (ETH), and 100 mM H 2 O 2 for 0, 6, 12, 24 h respectively. The treated leaf materials of X142 were collected and stored at −80 • C after flash-freezing in liquid nitrogen for further utilization. The matured fibers of different cotton varieties (G2, LO5R59, X142, J12, and XH36) were collected for length measurement by three independent experiments with 30 samples for each replicate.

RNA Extraction and qRT-PCR Analysis
Total RNA was extracted from different cotton tissues by our previous method as described [4,16]. The cDNA synthesis was obtained through reverse transcription reaction with 5 microgram total RNA as material according to the instruction of SuperScript ® III (Invitrogen, Carlsbad, CA, USA), and then used as a template for qRT-PCR detections with SYBR green real-time PCR master mixes (Applied biosystems, Foster, CA, USA). The specific primers were provided in Supplementary Table S5 and the relative expression levels of the target genes were calculated through the 2 −∆∆Ct method using the GhUBQ7 gene as internal control for normalization [94,95]. The relative expression levels of GhMDHAR genes measured by qRT-PCR were used to generate a heatmap that visualized by MultiExperiment viewer (MeV, version 4.9, J. Craig Venter Institute, La Jolla, CA, USA) software.

Enzyme Assays
The MDHAR enzyme activity was measured by a spectrophotometric method. The materials were ground to powder in liquid nitrogen, and about 50 mg of raw powder was added into 1 mL of pH 6.0 50 mM potassium phosphate buffer (including 1 mM AsA, 40 mM KCl, and 2 mM CaCl 2 ) that was then centrifuged at 14,000 rpm at 4 • C for 10 min. The reaction mix was obtained by the addition of collected enzyme extracts in pH 7.6 50 mM potassium phosphate buffer (containing 0.25 mM NADH, 2.5 mM AsA, and 15 µM FAD), and enough ascorbate oxidase (0.5 unit) was added into the mix solution to trigger the reaction. The values were detected at 340 nm for 2 min, and the absorbance decrease generated by the oxidation of NADH was used to calculate the MDHAR enzyme activity with the extinction coefficient of 6.22 mM −1 ·cm −1 [96]. The other regular antioxidant enzyme activities of ascorbate peroxidase (APX), superoxide dismutase (SOD), and catalase (CAT) were measured by commercial kits provided by the biotechnological company (Solarbio, Beijing, China).

Measurement of Asc, AsA and DHA
Asc measurement was performed as in our previous method [16]. Briefly, 5 g of material was ground in mortar on ice in 20 mL 50 g/L trichloroacetic acid (TCA) solution with successive dilution to 100 mL 50 g/L TCA solution. The collected supernatant by filtration was used for AsA determination. The reaction components included 1 mL 50 g/L TCA solution, 1 mL sample extract, 1 mL absolute ethanol, 0.5 mL 0.4% phosphate-ethanol, 1 mL 5 g/L BP-ethanol, and 0.5 mL FeCl 3 -ethanol solution. The total AsA content was generated by detecting the values at 534 nm through spectrophotometry according to the Asc content standard curve. The 0.5 mL of 60 mmol/L dithiothreitol (DTT) was added to the reaction mix for total Asc determination.

Institutional Review Board Statement: Not applicable
Informed Consent Statement: Not applicable Data Availability Statement: No new data were created or analyzed in this study. Data sharing is not applicable to this article.

Acknowledgments:
We thank Cotton Research Institute of Chinese Academy of Agricultural Sciences for providing the different fiber-length cotton cultivars (G2, LO5R59, and J12).

Conflicts of Interest:
The authors declare no conflict of interest.