Dose effects of restorer gene modulate pollen fertility in cotton CMS-D2 restorer lines via auxin signaling and flavonoid biosynthesis

Dose effects of Rf1 gene regulated retrieval mechanism of pollen fertility for CMS-D2 cotton. Cytoplasmic male sterility conditioned by Gossypium harknessii cytoplasm (CMS-D2) is an economical pollination control system for producing hybrid cotton seeds compared to artificial and chemical emasculation methods. However, the unstable restoring ability of restorer lines is a main barrier in the large-scale application of “three-line” hybrid cotton in China. Our phenotypic investigation determined that the homozygous Rf1Rf1 allelic genotype had a stronger ability to generate fertile pollen than the heterozygous Rf1rf1 allelic genotype. To decipher the genetic mechanisms that control the differential levels of pollen fertility, an integrated metabolomic and transcriptomic analysis was performed at two environments using pollen grains of four cotton genotypes differing in Rf1 alleles or cytoplasm. Totally 5,391 differential metabolite features were detected, and 369 specific differential metabolites (DMs) were identified between homozygous and heterozygous Rf1 allelic genotypes with CMS-D2 cytoplasm. In addition, transcriptome analysis identified 2,490 differentially expressed genes (DEGs) and 96 unique hub DEGs with dynamic regulation in this comparative combination. Further integrated analyses revealed that several key DEGs and DMs involved in indole biosynthesis, flavonoid biosynthesis, and sugar metabolism had strong network linkage with fertility restoration. In vitro application of auxin analogue NAA and inhibitor Auxinole confirmed that over-activated auxin signaling might inhibit pollen development, whereas suppressing auxin signaling partially promoted pollen development in CMS-D2 cotton. Our results provide new insight into how the dosage effects of the Rf1 gene regulate pollen fertility of CMS-D2 cotton.


Introduction
Cotton (Gossypium hirsutum L.) plays a key role to promote social and economic development worldwide . As an important cash crop in China, cotton productivity is frequently influenced by several challenges, while utilization of hybrid vigor can increase cotton productivity and improve fiber quality Shahzad et al. 2019). As an economically ideal pollination control system, cytoplasmic male sterility (CMS) can reduce seed production costs and improve seed purity, and it has been applied to generate hybrid cotton seeds (Yu et al. 2016). However, the ability to restore pollen fertility in CMS-D2 is generally influenced by the type of fertility restorer (Rf) genes, the nuclear background of restorer lines, and the external environment (Wu et al. 2017a;Zhang et al. 2020;  1 3 2022), which seriously hinders the large-scale application of "three-line" hybrid cotton in production.
With the development of molecular biology and highthroughput sequencing technology, several major Rf genes have been mapped and cloned for functional validation in rice (Hu et al. 2012;Huang et al. 2015;Jiang et al. 2022;Tang et al. 2014), maize (Qin et al. 2021), wheat (Melonek et al. 2021) and rape (Liu et al. 2016). Moreover, some studies have confirmed that Rf genes may exhibit dose effects on fertility restoration ability (Cai et al. 2013;Jiang et al. 2022;Melonek et al. 2021;Zhang et al. 2021). In rice, Rf 5 and Rf 6 genes can restore the fertility of HL-type indica CMS lines, and their dosage effects contribute to the revival of pollen fertility . Another study revealed that multiple alleles of Rf 3 and Rf 4 appeared to be responsible for variation in the pollen fertility of rice (Cai et al. 2013). In cotton, Rf 1 can restore the pollen fertility of both CMS-D2 and CMS-D8 lines, while Rf 2 can only restore the pollen fertility of CMS-D8 line Wu et al. 2011;Zhang and Stewart 2001). However, the cloning verification and fertility restoration mechanisms of these two restorer genes have not been reported so far. Our previous studies have found that the CMS-D2 sterile cytoplasm had negative effects on pollen fertility as well as seed cotton yield, and the homozygous Rf 1 Rf 1 allelic genotype SR showed stronger fertility restoration ability than the heterozygous Rf 1 rf 1 allelic genotype SH under high-temperature (HT) stress (Zhang et al. 2022a;Zuo et al. 2022). Unfortunately, the molecular basis for how the different alleles of the Rf 1 gene regulate pollen fertility in CMS-D2 cotton remains largely unclear. Therefore, further research on the dosage effects of the Rf 1 gene will help strengthen the selection and breeding of stable restorer lines in the future.
Considering the complexity of the molecular mechanism underlying pollen fertility restoration for CMS in flowering plants, many studies have already investigated the key metabolites and regulatory factors related to the growth of pollen tubes, double fertilization, and seed development (Gomez et al. 2015;Guo and Liu 2012). Major endogenous phytohormones including auxin (Cecchetti et al. 2008;Min et al. 2014), gibberellin acid (GA) (Chhun et al. 2007), and jasmonic acid (JA) (Fu et al. 2015;Khan et al. 2020) were reported to be involved in regulating pollen development and anther dehiscence. Starch accumulation is crucial for pollen maturation and viability during the late stages of pollen development (Wu et al. 2016). The differential levels of starch along with Cys proteases and AMS proteins are often associated with the reduction of pollen viability (Datta et al. 2002;Li et al. 2006;Sorensen et al. 2003). In addition, flavonoids are free radical scavengers, components of the pollen coat, and contribute to anther fertility (Filkowski et al. 2004;Hsieh and Huang 2007). Besides, an imbalance in lipid metabolism may disrupt anther cuticle and pollen development in plants (Ariizumi and Toriyama 2011;Shi et al. 2015;Zhang et al. 2022a). Research on the fertility restoration mechanism of CMS determined that excess ROS accumulation most likely caused pollen sterility in plants (Wan et al. 2007;Yang et al. 2018;Zhang et al. 2019Zhang et al. , 2022b. Recently, metabolomic and transcriptome analyses have become an important research strategy for developmental biology, and their integrated datasets provide useful insights into the regulatory mechanism of pollen development in CMS crops (Zhang et al. 2022a). However, the genetic determinants regulating pollen fertility restoration are not yet investigated more thoroughly in CMS-D2 cotton.
In this study, the pollen fertility of the SH genotype with one dominant Rf 1 allele (heterozygous Rf 1 rf 1 ) was found to be significantly lower than that of the SR genotype with two dominant Rf 1 alleles (homozygous Rf 1 Rf 1 ) in CMS-D2 cotton. An integrated metabolomic and transcriptome analysis was then performed at two environments using four cotton genotypes differing in Rf 1 alleles or cytoplasm. Our results identified how Rf 1 gene dosage influences the profiles of metabolites and transcripts in the pollen grains of various Rf 1 genotypes in cotton. Furthermore, the potential genetic determinants that revive pollen fertility in homozygous and heterozygous Rf 1 allelic genotypes of CMS-D2 cotton were also proposed. This study provides a new perspective for further elucidating the genetic mechanism of fertility restoration in CMS cotton.

Plant materials and growth conditions
In this study, four cotton near-isogenic lines (NILs) with homozygous and heterozygous Rf 1 alleles carrying normal upland cotton (AD1) and CMS-D2 sterile cytoplasm (denoted N and S, respectively) were used to investigate the dose effects of restorer gene on pollen fertility. Specifically, these four genotypes were named NR [N(Rf 1 Rf 1 )], NH [N(Rf 1 rf 1 )], SR [S(Rf 1 Rf 1 )], and SH [S(Rf 1 rf 1 )], and the breeding details were described in our previous studies Zhang et al. 2020Zhang et al. , 2022aZuo et al. 2022). All harvested seeds were conserved at the Cotton Heterosis Utilization Laboratory, the Institute of Cotton Research of Chinese Academy of Agricultural Science (ICR-CAAS). During the last week of April 2020, the seeds of selected materials were planted at the Baibi East breeding base of ICR-CAAS located in the Yellow River Basin cotton region (Anyang, Henan, China) as well as in the experimental field of the Cotton Research Institute of Jiang Xi Province located in the Yangtze River Basin cotton region (Jiujiang, Jiangxi, China), respectively. During the cotton full-bloom stage in summer, mature pollen samples were collected from ten plants and combined for each biological replication from NR, NH, SR, and SH both in Anyang and Jiujiang. All harvested samples were quickly frozen in liquid nitrogen and then stored at − 80℃ before further utilization.

Untargeted metabolomics data acquisition and mass spectrometry analysis
For metabolic profiling, frozen pollen for each sample with six biological replicates was first thawed on ice. About 1 g of mature pollen for each sample was used to make powder form and extracted with 120 μL of precooled 50% methanol, vortexed for 1 min, and then incubated at room temperature for 10 min. Later, the extraction mixture was stored overnight at -20 ℃. After centrifugation at 4000 g for 20 min, the supernatants were transferred into new 96-well plates. Metabolite accumulation was measured six times utilizing fully independent tissue samples. In addition, pooled QC samples were also prepared by combining 10 μL of each extraction mixture. All samples' quantitative metabolic data were acquired by the LC-MS system (www. lc-bio. com). The online Kyoto Encyclopedia of Genes and Genomes (KEGG) and HMDB database were used to annotate the metabolites' physical and chemical properties and biological functions. Meanwhile, an in-house fragment spectrum library of metabolites was further used to validate the metabolite identification. The significant differential metabolite features (DMFs) were screened with the fold change ≥ 2 or ≤ 0.5 between the target pollen samples, and the importance in projection (VIP) value ≥ 1 combined with q-value < = 0.05 based on the Benjamini-Hochberg test.

RNA extraction, library construction, and data analysis
Transcriptome sequencing was performed with three biological replicates on the same materials used for metabolic profiling (LC-Bio, Hangzhou, China). Total RNA from pollens for each sample was extracted using the TIANGEN RNAprep Pure Plant Plus Kit (Polysaccharides & Polyphenolics-rich; DP441) according to the vendor's protocol. After assessing the purity, quantity, and integrity of total RNA, the cDNA was then synthesized by SuperScript ™ II Reverse Transcriptase (Invitrogen, cat.1896649, USA). The final cDNA libraries were constructed following the protocol for the mRNA-Seq sample preparation kit (Illumina, San Diego, CA, USA). Finally, 2 × 150 bp paired-end sequencing (PE150) was performed on an Illumina NovaSeq ™ 6000 platform (LC-Bio, Hangzhou, China) following the manufacturer's recommended protocol.
After removing the low-quality reads that contained adaptor contamination and undetermined bases, the remaining clean reads were then mapped to the upland cotton TM-1 reference genome ) using HISAT2 (https:// ccb. jhu. edu/ softw are/ hisat2). After the comprehensive transcriptome was generated using gffcompare (https:// github. com/ gpert ea/ gffco mpare/), and then StringTie was used to assess the expression level for mRNAs via calculating FPKM. Further identification of differentially expressed genes (DEGs) between samples was performed using R package edge R (https:// bioco nduct or. org/ packa ges/ relea se/ bioc/ html/ edgeR. html) with fold change ≥ 2 or ≤ 0.5 along with an adjusted P-value < 0.05. The GOseq R package (Young, t al. 2010) and KOBAS software (Mao et al. 2005) were used for Gene Ontology (GO) functional categories analysis and to test the statistical enrichment of the DEGs in the KEGG pathways, in which GO terms or KEGG pathways with an adjusted P-value < 0.05 were considered to be significantly enriched.

Quantitative real-time PCR analysis
Total RNA was isolated from pollen samples using the TIANGEN RNAprep Pure Plant Plus Kit (Polysaccharides & Polyphenolics-rich; DP441), and one microgram (μg) of total RNA from individual replications was used for cDNA synthesis using a PrimeScript ™ RT reagent Kit with gDNA Eraser (TaKaRa, Dalian, China). Then, qRT-PCR analysis was performed on a Mastercycler ep realplex instrument (Eppendorf, Germany) using TransStart Top Green qPCR SuperMix (TransGen Biotech, Beijing, China). The relative expression level of genes was calculated using the 2 −ΔΔCt method as previously described (Wu et al. 2017b;Zhang et al. 2019). Expression levels were normalized by the G. hirsutum Histone3 ((GhHis3)) as an internal control to standardize RNA content, and each gene in each sample was analyzed with three biological replicates. Gene-specific primers of qRT-PCR are listed in supplementary Table S1.

Exogenous application of the IAA inhibitor and promoter
To explore the potential role of IAA in regulating cotton anther development, the flower buds of all four cotton NILs were sprayed with 10 and 100 μM NAA, and deionized water was sprayed as a negative control (Mock). Besides, SR and SH were also pretreated with 20 μM IAA inhibitor Auxinole. The above different concentrations of IAA inhibitor and promoter solutions were sprayed onto all buds of four cotton lines in the afternoon once every 2 days, with at least ten plants per treatment. During spraying, the anther morphology of each material for each treatment was observed and recorded simultaneously after the first treatment; that is, the number of flowers with fewer pollen grains or dehiscence anther, and the total flowers were counted, and thus finally the percentage of flowers with fewer pollen grains was calculated.

Phenotypic analysis and determination of pollen viability
Pollen grains from representative flowers were stained with 0.5% 2,3,5-triphenyltetrazolium chloride (TTC) (Solarbio, Beijing, China) solution to observe pollen viability, as described in detail in our previous study . Specificity, TTC reacts with dehydrogenase in normal anther tissue and turns red. Thus, viable pollen appeared red, and partially viable pollen appeared reddish, whereas dead as well as sterile pollen appeared colorless. Images of TTC staining pollen grains were captured under a bright field using an Olympus SZX16 research stereo microscope system (https:// lifes cience. evide ntsci entifi c. com. cn/ en/ micro scopes/ stereo/ szx16/).

Statistical analysis
Results in this study were presented as the means of three biological replicates ± the standard deviation (SD) and bar graphs were displayed with the GraphPad Prism 8 software. Statically significance was calculated among the different treatments and Mock using the two-way ANOVA followed by the least significant difference (LSD) test, and values shown with different superscript letters were considered significantly different (P-value < 0.05). The statistical significance analysis of gene expression levels was conducted using a two-tailed unpaired Student's t-test, and a P-value < 0.05 was considered a statistically significant difference.

Phenotypic comparison of male fertility in four cotton near-isogenic lines (NILs) at two environments
Our previous study has found that the anther fertility of the sterile cytoplasmic restorer lines SH and SR was significantly lower than that of the normal upland cotton cytoplasmic restorer lines NH and NR under HT stress. Specifically, sterile cytoplasmic restorer line SR with two dominant Rf 1 (Rf 1 Rf 1 ) genotypes presented higher pollen fertility than SH with a single dominant Rf 1 (Rf 1 rf 1 ) genotype (Zuo et al. 2022). Here, the representative anther morphology and pollen viability among four cotton NILs were compared at both Anyang (AP) and Jiujiang (JP) environments. Obviously, SH and SR showed a decreased filament length and pollen viability, and an increase in the exposed length of stigma in both environments, especially in the Yangtze River Basin cotton region where the summer temperature is higher (Fig. 1). However, there was no obvious difference in the external morphology and size of the intact flower among these four cotton NILs (Fig. 1a, b). Considering that in the CMS-D2 sterile cytoplasm, SH with heterozygous Rf 1 allele showed an obvious reduction in pollen amount and pollen viability than SR with homozygous Rf 1 allele, whereas there was no obvious difference in pollen fertility between NH and NR with the same normal upland cotton (AD1) cytoplasm but different Rf 1 alleles (Fig. 1c-j). These findings suggest that the dose of the restorer gene may be involved in differences in pollen fertility restoration. Thus, a comprehensive comparative analysis of metabolomic and transcriptome sequencing data was conducted to further explore the molecular basis of how the Rf 1 dosages affect pollen development in cotton CMS-D2 restorer lines.

Overview of metabolite abundance in various Rf 1 genotypes of CMS-D2 cotton
The quantitative and qualitative metabolic data among four cotton NILs were first comparatively analyzed in both environments to identify key metabolic substances. Principal component analysis (PCA) and correlation analysis based on metabolites revealed obvious differences among studied samples (Fig. S1). In PCA, both PCA1 and PCA2 accounted for 27.27% and 16.35% of sample variation, respectively. PCA with all metabolite features showed SH had a distinct cluster from other genotypes. In addition, metabolomic profiles of different samples showed significant differences under both environments (Fig. S1a). The correlation analysis further indicated that metabolomic profiles displayed good repeatability and reliability (Fig. S1b).
A total of 5,391 differential metabolite features (DMFs) were detected in both environments. Of which,530,579,195, and 1,678 metabolite features had significantly higher abundance in AP_NH vs AP_NR, AP_SH vs AP_SR, JP_ NH vs JP_NR, and JP_SH vs JP_SR, respectively. Similarly, totally 664, 864, 152, and 1,449 metabolite features showed lower abundance among the above four various comparisons, respectively (Fig. 2a, Table S2). The distribution of overlapped or specifically accumulated metabolite features indicated dynamic changes in each Rf 1 genotype. A total of 369 metabolites overlapped in SH vs SR under both environments were identified, of which 147 DMFs had higher abundance while 222 presented lower abundance in SH compared with SR (Fig. 2b, c). In addition, heat map analysis also showed the quantitative abundance of these metabolic substances was significantly different in SH in each environment (Fig. 2d). The majority of differential metabolites (DMs) specific in SH vs SR belonged to fatty acyls, carboxylic acid and derivatives, indoles and derivatives, and flavonoids (Fig. 2e). It is noteworthy that specific DMs in SH vs SR were significantly enriched to glycerolipid metabolism, linoleic acid metabolism and plant hormone signal transduction (Fig. 2f), indicating the restorer gene Rf 1 dosages may affect pollen development through interaction with the sterile cytoplasm caused large differences in the composition and concentration of metabolic substances involved in lipid metabolism, flavonoid metabolism, and auxin signaling pathways.

Global dynamics of transcriptional profiling in various Rf 1 genotypes of CMS-D2 cotton
RNA sequencing was also performed with three biological replicates on the same pollen samples of four cotton NILs used for metabolic profiling to analyze the profiles of transcripts in different environments. A total of 1.038 billion raw reads were generated from 24 samples with an average read length of 150 bp. After stringent quality checks followed by data filtering, an average of 6.41 Gb valid data were obtained for each sample, and the ratio of valid data and sequencing Q30 values of all samples were greater than 98%, indicating the accuracy of the data obtained in this study (Table S4). PCA and correlation analysis based on all the identified genes revealed obvious differences among studied samples (Fig. S2). In PCA, SH clustered separately from other Rf 1 genotypes that stated an obvious difference, especially in the JP environment of the Yangtze River Basin cotton region (Fig. S2a). In addition, the correlation analysis among different samples further indicated the reliability and good repeatability of sampling (Fig. S2b).
RNA-seq analyses showed that a large transcriptome reprogramming occurred in all Rf 1 genotypes of CMS-D2 cotton (Table S5). A total of 2,490 differentially expressed genes (DEGs) were identified in the four pairwise comparisons. Specifically, 100 and 235 DEGs were up-and down-regulated in AP_NH vs AP_NR, whereas 243 upand 100 down-regulated DEGs were identified in JP_NH vs JP_NR. Likewise, AP_SH vs AP_SR and JP_SH vs JP_SR comparative combinations showed more DEGs, that is, 561 and 1,484 total DEGs, respectively, of which 473 up-and 88 down-regulated DEGs were identified in the AP environment, while 1,305 and 179 DEGs were upand down-regulated in JP environment (Fig. 3a, Table S5). Furthermore, a total of 96 DEGs unique to SH vs SR comparative combination under both environments were also identified, which may be involved in Rf 1 dosage to regulate pollen fertility of CMS-D2 restorer lines (Fig. 3b, c). Hierarchical cluster analysis showed the expression profiles of Fig. 1 Comparison of anther phenotype and pollen activity of different Rf 1 genotypes at two environments. a, b Representative phenotype of anthers for NR, NH, SR, and SH in Anyang (AP) and Jiujiang (JP), respectively. c-f The 0.5% 2,3,5-triphenyltetrazolium chloride (TTC) stained pollen grains of NR, NH, SR, and SH in AP. g-j The 0.5% TTC stained pollens of NR, NH, SR, and SH in JP. Scale bars: 1 cm (red) in a and b; 250 μm (white) in c-j several specific DEGs were significantly higher or lower in SH than in other Rf 1 genotypes. The majority of genes annotated with LHY, PME, At3g48460, PCMP-E32, SS2, and AGPS1 had shown higher regulation expression in SH. Whereas the HSP genes that regulate stress response exhibited significant down-regulation in SH compared with other Rf 1 genotypes in both environments (Fig. 3d). Subsequently, we further conducted GO functional classification and KEGG pathway enrichment analysis on these 96 specific DEGs (Fig. S3, Tables S6 and S7). Many specific DEGs showed functional annotation to the regulation of transcription, response to abscisic acid, nucleus, cytoplasm, and binding to ATP (Fig. S3a, Table S6), while had pathways enrichment in 'Flavonoid biosynthesis', 'Circadian rhythm-plant', and 'Indole alkaloid biosynthesis' (Fig. S3b, Table S7).

Regulatory network associated with pollen fertility in various Rf 1 genotypes of CMS-D2 cotton
The association between transcripts and metabolites permits the identification of biological networks of target traits, and the final part of the co-expression network is shown in Fig. 4a. Here, several differential genes and metabolites involved in indole biosynthesis, sugar metabolism, and flavonoid biosynthesis had shown strong network linkage with pollen fertility. These can be, therefore, considered as hub genes and metabolites that influenced the mechanism of Fig. 2 Characterization of differential metabolites (DMs) in pollen of different Rf 1 genotypes. a Statistics of differential metabolite features (DMFs) among different Rf 1 genotypes. b Distribution of up-regulated DMFs among various comparisons. c Distribution of downregulated DMFs among various comparisons. d Heat map analysis showing the quantitative abundance of key DMFs in six biological repeats of different Rf 1 genotypes. e Classification of specific differential metabolites (DMs) identified in SH vs SR. f KEGG pathway enrichment analysis of the specific DMs identified in SH vs SR pollen devolvement from pollen germination to maturity. Most of the hub metabolites were higher accumulated in SH than in SR while their abundance in NH was similar to that of NR, especially in the JP environment of the Yangtze River Basin cotton region (Fig. 4b). In particular, auxin pathway compounds such as 1H-Indole-1-carboxamide, 6-chloro-2,3-dihydro-5-methyl-N-[6-[(2-methyl-3-pyridinyl)oxy]-3-pyridinyl] showed significantly higher abundance in SH. Among flavonoid compounds, isorhamnetin, patuletin, and quercetin 3-(6'-malonyl-glucoside) had higher accumulation in SH than that in SR. In addition, erlose and palmitoyl ethanolamide involved in the sugar metabolism pathway presented higher abundance in SH (Fig. 4b, Tables S2 and S3). Further qRT-PCR analysis confirmed that the expression levels of hub genes such as SARF4, G9 and AACT , PME28 and PME58 involved above three key pathways were significantly higher in SH than other genotypes ( Fig. 4c-g). In brief, the disruption of metabolites and their regulatory genes most likely produced discrepancies in biochemical and molecular processes linked with anther development. Furthermore, these changes might cause dysfunction in the interactions among the nucleus Rf 1 gene and sterile cytoplasmic genome. This finally led to the difference in pollen fertility and viability in SH as compared to SR. Further in-depth research on hub genes and metabolites can be helpful to understand the genetic control of pollen fertility in CMS-D2 cotton.

Exogenous auxin treatment inhibits pollen development while inhibitor partially promote pollen development in CMS-D2 cotton
Molecular evidence has shown that auxin regulates pollen germination and pollen tube growth in plants . To further explore the potential role of auxin in regulating pollen fertility of CMS-D2 cotton, exogenous 10 and 100 μM auxin analogue NAA were applied to the flower buds of four cotton NILs, and 20 μM auxin inhibitor Auxinole were also simultaneously applied to SH and SR in vitro. As expected, the percentage of flowers with fewer pollen grains in four cotton NILs showed a statistically significant increase after treatment with different concentrations of NAA, especially in 100 μM NAA treatment, compared with the corresponding Mock (Fig. 5). Conversely, compared with the Mock, the percentage of flowers with fewer pollen grains presented a statistically significant decrease only in SH after the 20 μM Auxinole treatment (Fig. 5b). Moreover, there were significant differences in the percentage of flowers with fewer pollen grains between different Rf 1 genotypes, such as between NR and NH, as well as between SR and SH, under low concentrations of 10 μM NAA treatment, but no significant differences were found between them under higher concentrations of 100 μM NAA treatment (Fig. 5). These demonstrated that exogenous auxin Fig. 4 Identification of key genes and metabolites involved in Rf 1 dosage to influence pollen fertility of CMS-D2 cotton. a Overview of network linkage among the key DMFs and DEGs, and their involved pathways. The circle represents the DMFs; the rectangle represents the pathway; the triangle represents DEGs. The red, green, orange, and blue colors are involved in IAA biosynthesis, flavone biosynthesis, sugar metabolism, and other pathways, respectively. b Heat map showing abundance levels of key metabolic compounds related to IAA biosynthesis, flavonoid biosynthesis, and sugar metabolism. c-g qRT-PCR analysis validating the expression levels of key genes involved in IAA biosynthesis (c), flavonoid biosynthesis (d, e), and sugar metabolism (f, g) treatment inhibited pollen development, whereas the application of auxin inhibitor Auxinole significantly improved the pollen fertility in SH, indicating the balance of auxin may be necessary to decipher pollen sterility in CMS-D2 cotton.
To further reveal the role of auxin signaling on pollen fertility, expression profiles of ten key genes were further confirmed in the pollen grains of treated plants through qRT-PCR analysis. Most of the selected genes were annotated to auxin-responsive genes and encoded as GH3, AUX|IAA, and SAUR family genes (Fig. S4). The analysis showed that most GH3 and AUX|IAA-related genes had significant differential expression in NAA-treated pollen grains compared to the Mock (Fig. 6). In pollens of NR and NH, GH3, AUX|IAA, and TIR1 family genes showed different degrees of higher expression, while SAUR genes presented a significant downward expression (Figs. 6a-e, S4). For SR and SH, the expression of GH3.17 and AUX22D displayed a gradual increase trend, but SAUR32 and IAA14-2 showed a downward trend after different concentrations of NAA treatments. Conversely, the expression of GH3.17 and IAA14-2 decreased after auxin inhibitor Auxinole treatment (Fig. 6f-j). These findings suggest that auxin balance is indeed indispensable for pollen development in cotton. In the case of sterile cytoplasm, the difference in Rf 1 dosage may disrupt the balance of auxin levels essential for pollen fertility restoration. This may be the main reason that over-activated auxin signaling suppresses the pollen development both in SR and SH, whereas suppressing auxin signaling partially promote pollen development in CMS-D2 cotton (Figs. 5,. In addition, we further determined a significant difference in expression profiles of flavonoid-related genes AACT and G9 and sugar metabolism-related genes PME58 and PME28 in response to the above-mentioned auxin treatments (Figs. 4,  S5). However, further research would be necessary to dissect how the interaction among various metabolic pathways involved has synergistic or antagonistic effects on modulating the mechanism of fertility restoration in CMS-D2 cotton.

The Rf 1 gene does have dosage effects on pollen fertility restoration for CMS-D2 cotton
The combination of CMS and restorer lines is indispensable for the development of elite "three-line" hybrid varieties (Kim and Zhang 2018), but the restoring abilities of restorer lines largely depend on genetic background and diversity of fertility restorer alleles (Cai et al. 2013;Jiang et al. 2022;Melonek et al. 2021;Zhang et al. 2021). For example, the Rf 3 and Rf 4 can restore the fertility of the wild-abortive type CMS (CMS-WA) in rice, and the genotype with Rf 3-4 Rf 3-4 /Rf 4-4 Rf 4-4 possessed the strongest restoring ability, and genetic effects of Rf 4 alleles appeared to be strong than that of Rf 3 (Cai et al. 2013). Similarly, the Rf 5 and Rf 6 contribute to the fertility restoration process in Honglian (HL)-type japonica CMS lines. The lowest fertility was observed in lines with the Rf 5 rf 5 r f 6 rf 6 genotype, whereas the Rf 5 Rf 5 Rf 6 Rf 6 genotype showed the highest ability to produce fertile pollen. Therefore, the additive and dosage effects of Rf 5 and Rf 6 controlled the percentage of fertile pollen in CMS rice . Another research determined that rice hybrids harboring two restorer genes have a more stable seed-setting rate than plants containing only one Rf gene (Zhang et al. 2017). Recent research has identified multiple Rf genes in the genome of chili pepper. However, plants homozygous for the recessive Rf 1 (rf 1 rf 1 Rf 2 Rf 2 ) produced lower pollen grains compared to Rf 2 (Rf 1 Rf 1 rf 2 rf 2 ). This may be due to that Rf 1 is the main restorer gene while Rf 2 is the minor restorer gene in chili pepper (Zhang et al. 2022c). In cotton, Rf 1 and Rf 2 are identified as the main fertility restorer genes. The dominant Rf 1 can restore pollen fertility in both CMS-D2 and CMS-D8 systems, whereas the dominant Rf 2 can only restore the fertility of CMS-D8 (Kohel  . 1984;Meyer 1975;Weaver and Weaver 1977;Zhang and Stewart 2001). Our previous field evaluation observed allelic differentiation in the Rf 1 gene caused variation in pollen fertility and pollen germination rate in CMS-D2 cotton (Zuo et al. 2022). Consistently, our investigation further confirmed a homozygous Rf 1 Rf 1 (SR) genotype showed higher pollen fertility than heterozygous Rf 1 rf 1 (SH) genotypes in CMS-D2 cotton at both environments (Fig. 1a, b). These various Rf 1 allelic genotypes can be the ideal material to study the interaction between the nucleus and the cytoplasmic genome. In addition, these restorer lines will be useful to improve the efficiency of "threeline" hybrid breeding in cotton.

The dose effects of Rf 1 altered the dynamics of metabolites and transcripts in pollen of CMS-D2 cotton
The restorer lines differing in Rf 1 alleles with CMS-D2 or upland cotton cytoplasm can offer a better platform to understand the functional mechanism of fertility restoration in CMS-D2 cotton. This study compared the metabolites and transcripts profiles in pollen grains of four cotton genotypes containing homozygous and heterozygous Rf 1 alleles. Specifically, our results revealed a significant difference in metabolite substances along with gene regulation between the heterozygous and homozygous Rf 1 allele genotypes of CMS-D2 (Figs. 2 and 3). Importantly, the Fig. 6 The qRT-PCR analysis validating the transcript levels of auxin-responsive genes in pollen grains of different Rf 1 genotypes after auxin analogue NAA and its inhibitor Auxinole treatments. Here NR and NH were treated with 10 and 100 μM NAA, while SR and SH were treated with 10 and 100 μM NAA, and 20 μM Auxinole integrated metabolomic and transcriptomic analysis further uncovered key DEGs and DMs involved in indole alkaloid and flavonoid biosynthesis pathways between homozygous and heterozygous Rf 1 allele combinations (Fig. 4). Since Rf genes opted for diverse ways to restore fertility in plants, the key roles of auxin and flavonoids are comprehensively discussed in the pollen fertility of CMS-D2 cotton. Auxin, as an important hormone, is an important regulator of growth and development in plants and can influence sexual reproduction including the development of stamens, gynoecia, and ovary. It further promotes the maturation of egg cells along with the polar development of the embryo (Aloni et al. 2006;Mol et al. 2004;Nemhauser et al. 2000). IAA functional activities were also found to influence pollen tube growth in Torenia fournieri (Wu et al. 2008). In addition, it can initiate anther dehiscence (Cecchetti et al. 2008), affect stamen development, and its flow can improve the elongation of stamen filament (Hirano et al. 2008). In this study, IAA-related compounds had a higher abundance in the heterozygous Rf 1 allele genotype than homozygous Rf 1 allele genotype of CMS-D2, but there was no significant change in two cotton lines with the normal upland cotton cytoplasm (Fig. 4b). This means allelic differentiation in the Rf 1 gene might generate various species of auxin compounds in sterile cytoplasm. These qualitative and quantitative metabolic changes ultimately lead to shorter filaments and lower fertility in Rf 1 rf 1 genotypes of CMS-D2 compared to the other genotypes (Fig. 1).

Auxin signaling modulates the retrieval of pollen fertility for CMS-D2 cotton
Auxin-responsive genes were broadly grouped into three major classes including AUX|IAA, SAUR , and GH3 (Guilfoyle 1999), and have been shown to participate in the regulation of anther dehiscence and pollen development (Min et al. 2014;Zhou et al. 2015). Previous studies revealed ARF17 regulated the expression of CalS5 and this gene was found to be essential for pollen wall formation (Yang et al. 2013). The overexpression of GH3.9 appeared to reduce plant height, silique size, and stamen length in Arabidopsis . Importantly, SAUR39 acts as a negative regulator for auxin synthesis and transportation, and overexpression of SAUR39 promotes the senescence of leaves and inhibits growth and yield in rice (Kant et al. 2009). Overexpression studies on AtIAA31 reported that it may cause early development arrest of the shoot apical meristem and ApAux/IAA3 plants exhibited similar auxin-related aberrant phenotype and delay growth (Sato and Yamamoto 2008;Yang et al. 2019). Besides, the genes such as GH3, SAUR , and AUX|IAA have been reported to induce negative effects on pollen development (Bemer et al. 2017;Kant et al. 2009;Sato and Yamamoto 2008;Zhou et al. 2015). In this study, many auxin-responsive genes such as GH3, SAUR , and AUX|IAA were up-regulated in the heterozygous Rf 1 allele genotype (Fig. S4). These results suggest that the allelic differentiation in the Rf 1 gene may contribute to pollen fertility by modulating the appropriate auxin level. In brief, the allelic differentiation in Rf 1 most probably mediated dynamic changes of GH3, SAUR , and AUX|IAA encoding genes in CMS-D2 cotton. On the other way, excessive accumulation of auxin caused the up-regulation of auxin-responsive genes in the heterozygous Rf 1 allele genotype, which ultimately resulted in reduced pollen viability. This finding was consistent with the previous study in cotton (Min et al. 2014). Our results further confirmed that exogenous auxin treatment not only produced a higher number of flowers with fewer pollen grains in the heterozygous Rf 1 allele genotype but also an imbalance of the expression of auxin-responsive genes. In contrast, the application of Auxinole reduced the percentage of flowers with fewer pollen grains in the heterozygous Rf 1 allele genotype as well as altered the expression of auxin target genes (Fig. 5, 6). In addition, it has been stated that circadian rhythms control auxin response genes in the plant (Covington and Harmer 2007). Wu determined that the circadian rhythm pathway differs between CMS-D2 and its fertile lines (Wu et al. 2017b). Both auxin and the circadian clock annotated genes play pervasive roles to mediate various metabolic functions linked with the mechanism of flowering and pollen fertility in CMS cotton (Hocq et al. 2017;Kim et al. 2017;Sanchez et al. 2011). In the Rf 1 rf 1 genotype of CMS-D2, the higher expression of circadian clock LHY genes might cause flowering growth dysfunction and finally result in a higher percentage of sterile pollen phenotype (Fig. 3d). Taken together, we, therefore, infer that the Rf 1 allelic effects may regulate pollen fertility of CMS-D2 cotton via auxin signaling.

Flavonoids along with sugars facilitate the revival of pollen fertility for CMS-D2 cotton
Flavonoids are a larger group of secondary metabolites, abundant in mature pollen, and it has been hypothesized that flavonoids protect nucleic acids in pollen (Pacini and Hesse 2005;Schijlen et al. 2004;Winkel-Shirley 2001). Moreover, flavonoids are important signaling molecules, fertility regulators, and auxin transporters in plants. Previous studies have shown the roles of flavonoids in male fertility and sexual reproduction in many plant species (Kong et al. 2020;Mo et al. 1992;Wang et al. 2020). MYB transcription factors play important roles in the regulation of gene expression during plant growth and mainly participate in primary and secondary metabolism, including anthocyanin and flavanols biosynthesis (Gonzalez et al. 2008;Stracke et al. 2007). Overexpression of cotton GhMYB24 in Arabidopsis caused flower malformation, shorter filaments, 1 3 non-dehiscent anthers, and fewer viable pollen grains (Li et al. 2013). Consistent with previous research, our results detected that several metabolites, as well as genes linked with flavonoids components, had significant differential regulation between homozygous and heterozygous Rf 1 allelic genotypes of CMS-D2 (Fig. 4b). In particular, the dosage effects of Rf 1 alleles may cause changes in the regulation of flavonoids related genes as well as the composition of flavonoids substances (Fig. 4b, d, e). This most likely breaks the ROS balance that results in complex biological disorders during anther development and produces lower fertile pollen in heterozygous Rf 1 allelic genotypes of CMS-D2 (Fig. 1). As high oxidative stress is prevailing in response to disruption in flavonoids, therefore, comprehensive research on flavonoid compounds could be crucial to explore the genetic architecture of fertility restoration in CMS-D2 cotton. In plants, sugar derivatives such as pectin, starch, and cellulose are the basic source of energy and structural constituents for plant cells (Shi et al. 2015;Yu et al. 2015). The pectin polymers, cellulose, and hemicellulose are considered the main component of pollen walls (Hasegawa et al. 2000). The inhibition of pectin formation and degradation leads to a delay in pollen development, partial male infertility, and reduced fruiting rates (Wei et al. 2019;Zhang et al. 2010). In the heterozygous Rf 1 allele genotype, the higher regulation of pectin-encoding genes including PME21, PME28, PME41, and PME58 maybe weaken the impact of Rf 1 rf 1 to restore complete pollen fertility via disruption in the level of pectin during pollen wall formation (Fig. 4 a, f, g). Collectively, we deduce that the Rf 1 dosage may regulate pollen fertility restoration for CMS-D2 cotton through flavonoid biosynthesis along with sugar metabolism.

Potential mechanism of Rf 1 dose effects on pollen fertility restoration for CMS-D2 cotton
Based on our findings of this research along with those of previous studies, we proposed a potential regulatory model showing how allelic differentiation of Rf affects the pollen fertility in CMS-D2 cotton (Fig. 7). In the CMS-D2 cotton system, nucleo-cytoplasmic interaction between Rf 1 and orf610a mediate the balance of ROS production and energy homeostasis to restore normal pollen development (Zhang et al. , 2022b. In addition, the Rf 1 gene may significantly alter the profiles of key genes and metabolites, such as auxin biosynthesis, flavonoid biosynthesis, and sugar Fig. 7 A potential model explaining how allelic differences in the Rf 1 gene influence pollen fertility restoration for CMS-D2 cotton. The instability in nucleo-cytoplasmic interaction between orf610a and Rf 1 alleles may happen in response to an imbalance of auxin, flavonoid, and sugar substances. However, some interactions in this net-work remain still unclear. The lines with arrows and blunt ends in the figure indicate the promotion and inhibition modes, respectively, and the accompanying question marks represent unknown action modes or connections metabolism that may be tightly interlinked with pollen fertility restoration for CMS-D2 cotton. In the heterozygous Rf 1 rf 1 genotype SH with CMS-D2 cytoplasm, excessive accumulation of auxin caused over-activated auxin signals, which may ultimately lead to lower pollen fertility by promoting flavonoid synthesis and inhibiting sugar metabolism. Comparatively, the homozygous Rf 1 Rf 1 allelic genotype SR probably has a strong ability to maintain stable nucleo-cytoplasmic interaction with orf610a via modulating appropriate auxin signaling, and the various compounds related to auxin, flavonoids, and sugars possibly maintain energy homeostasis to generate normal fertile pollen. However, some hypothetical interactions in the regulatory network of how Rf 1 dosage regulates pollen fertility remain largely indistinct, and further in-depth experiments are still needed to explore.

Conclusions
In summary, the present study performed an integrated metabolome and transcriptome analysis in diverse Rf 1 genotypes and uncovered that Rf 1 allelic differentiation affected pollen fertility by altering the landscape of transcripts and metabolite substances. In CMS-D2 sterile cytoplasm, the predominant changes between homozygous Rf 1 Rf 1 and heterozygous Rf 1 rf 1 genotypes were found in pathways including auxin biosynthesis, flavonoid biosynthesis, and sugar metabolism. Further, in vitro application of auxin promoter and inhibitor validated that over-activated auxin signaling could inhibit pollen development while reducing auxin signaling partially promoted pollen development in CMS-D2 cotton. Our results revealed how the dosage effects of the Rf 1 gene regulate pollen fertility of CMS-D2 cotton, and it will help strengthen the selection and breeding of restorer lines with stable fertility in production.