The G. barbadense line HaiR, G. hirsutum line Y18R, and red color cotton line pcc, were used in this study. The red-petaled pcc cotton was selected from over 1,000 varieties in the cotton germplasm resource library. The Y18R line is characterized by white flowers, whereas the pcc line exhibits red flowers. The HaiR line is distinguished by its yellow flowers. For the generation of F₂ mapping populations, HaiR and pcc were selected as the parental lines.

Cotton plants were grown under field conditions at the High-Tech Industry Park of the Chinese Academy of Agricultural Sciences in Langfang, Hebei Province (N39°520, E116°700E). Photographic documentation and petal collection were conducted during and post-flowering for subsequent analyses. Quantitative Reverse Transcription Polymerase Chain Reaction (qRT-PCR) and anthocyanin content determination were performed, with each assay conducted in triplicate.

Map-based cloning was conducted using an F₂ population consisting of a total of 870 cotton plants, which exhibited two distinct phenotypes: yellow-flowered and red-flowered. For genetic linkage analysis, 270 individuals displaying the yellow color flower phenotype were selected at the flowering stage. A total of 720 pairs of insertion/deletion (InDel) primers, designed based on genome-wide InDel markers differentiating G.barbadense and G.hirsutum genotypes, were used for preliminary gene mapping (He et al., 2018). The PCR amplifications were electrophoresed on a 4% agarose gel. Primers utilized for fine mapping are detailed in Supplementary Table 3 .

Building on previous research into the metabolic pathways of anthocyanins and proanthocyanidins (Liu et al., 2021; Pucker and Selmar, 2022; Sunil and Shetty, 2022), and we made minor adjustments and illustrated these pathways using Figdraw to accurately depict the metabolic regulation model.

For expression analysis, samples were collected and immediately flash-frozen in liquid nitrogen, then stored at -80°C until RNA extraction. Total RNA was extracted using the RNAprep Pure Plant Kit (TIANGEN, DP441, Beijing, China) following the manufacturer’s protocol. cDNA synthesis was performed using the TransScript One-Step gDNA Removal and cDNA Synthesis SuperMix (TRAN, AT311, Beijing, China) as per the manufacturer’s guidelines. qRT-PCR analysis was conducted using the Chromo 4 Real-Time PCR Detection System (Bio-Rad, CFX96, Hercules, USA), adhering to the manufacturer’s instructions. The cotton GhHistone gene served as the internal reference. Experiments were independently replicated three times. Primers for qRT-PCR are detailed in Supplementary Table 3 . The relative expression levels of genes were quantified using the 2^-ΔΔCT method (Livak and Schmittgen, 2001).

The raw RNA-seq datasets for red and white flowers of G. hirsutum were obtained from PRJNA878950, while those for brown and white cotton fibers were sourced from PRJNA766762. Following the removal of low-quality reads, the remaining clean reads were aligned to the TM-1 reference genome using Hisat2 software (Kim et al., 2019). Differential gene expression analysis was conducted using featureCounts for the calculation of Transcripts Per Million (TPM) values (Liao et al., 2014; Zhao et al., 2020). Gene expression heatmaps were generated utilizing the Pheatmap package (Hu, 2021).

The determination of anthocyanin content followed the methodology described by Xu et al. (2012). Petal samples (0.1 g) were retrieved from a -80°C freezer, ground to a powder in liquid nitrogen, and then mixed with 1 mL of methanol. The mixture was vortexed thoroughly and sonicated for 15 minutes. After centrifugation at 12,000 rpm for 10 minutes, the supernatant was diluted tenfold with methanol and filtered through a 0.22 μm organic membrane. LC-MS analysis was performed using a Thermo U3000 HPLC System. The peak area of the mixture was quantified by comparison with the peak areas of standard substances, including cyanidin, delphinidin, and pelargonidin (all from Sigma, USA). The formula for calculating anthocyanin content is as follows: [A = (C2 - C1) * V * N where (A) represents the content of various anthocyanins in the sample (mg/kg), (C2) is the concentration of each anthocyanin in the test solution (mg/L), (C1) is the concentration of each anthocyanin in the blank solution (mg/L), (V) is the volume of the solution (mL), and (N) is the dilution factor.

Total 96 accessions derived from four projects (PRJNA349094, PRJNA680449, PRJNA412456, PRJNA720818) were analyzed (see Supplementary Table 4 ). The sequencing reads were aligned to the G. hirsutum TM-1 reference genome using BWA software (Li and Durbin, 2009). The resulting BAM files were processed for sorting and filtering using SAMtools (version 1.15) with default parameters (Danecek et al., 2021). SNP calling at the population level was conducted using GATK software (version 4.0) (Van der Auwera et al., 2013), generating a GVCF file for each sample. These GVCF files were then merged, applying the following parameters for filtering: –filter-expression “QD < 2.0 || MQ < 40.0 || FS > 60.0 || SOR > 3.0 || MQRankSum < -12.5 || ReadPosRankSum < -8.0” –filter-name ‘SNP_filter’.

A total of 120,126 SNPs, with a minor allele frequency (MAF) > 0.01 and missing data < 10%, were identified using Plink (version 1.9). These SNPs were used to construct a phylogenetic tree using the neighbor-joining method in PHYLIP software (version 3.698) (Retief, 2000). The phylogenetic tree was visualized and refined using iToL (https://itol.embl.de/). Population genetic structure was assessed using Admixture (Version 1.3.0), with the number of genetic clusters (K) ranging from 2 to 6 and each run consisting of 10,000 iterations.

Using VCFtools (version 0.1.16) (Danecek et al., 2011), we calculated fixation index (Fst) and nucleotide diversity (π) employing a sliding window approach. Each window spanned 50kb, with successive windows overlapping by 10kb. This methodology facilitates a detailed and nuanced analysis of genetic differentiation and diversity across the genome.

Genomic sequences, annotation files, and protein sequences for Gossypium species including G. herbaceum (A1), G. arboretum (A2), G. raimondii (D5), G. hirsutum (AD1), G. barbadense (AD2), G. tomentosum (AD3), G. mustelinum (AD4), and G. darwinii (AD5) were retrieved from Phytozome (Goodstein et al., 2012) and CottonMD (Yang et al., 2023). Protein sequences for the GST and MYB gene families were downloaded from the TAIR database (https://www.arabidopsis.org/). Hidden Markov Model (HMM) profiles for GST_C (PF00043), GST_N (PF02798), and MYB (PF00249) were acquired from the Pfam database (version 35.0).

To identify homologs of the MYB and GST gene families across the Gossypium species, Arabidopsis thaliana GST and MYB protein sequences served as query sequences in BLASTp searches (version 2.15), with stringent criteria set at e-values < 1e-10 and identity ≥ 30%. Concurrently, HMMER (version 3.4) software, utilizing the hmmersearch function and the aforementioned HMM profiles, was employed to search for protein sequences within the cotton varieties, applying an e-value threshold of ≤ 1e-5. Results from both BLASTp and hmmersearch were integrated, allowing for the exclusion of structurally incomplete sequences from the GST and MYB gene families in the analyzed Gossypium species.

Protein sequences from various Gossypium species were consolidated and subjected to multiple sequence alignment using MUSCLE (version 5.0). Subsequently, a multi-species evolutionary tree was constructed employing the Maximum Likelihood (ML) method via FastTree (Price et al., 2009). The resulting evolutionary tree was visualized and refined using the online tool iTOL (Letunic and Bork, 2021). Within this phylogenetic framework, we specifically identified protein sequences associated with anthocyanin biosynthesis. These sequences were then realigned using MUSCLE (version 5.0), and the alignment was further refined for accuracy and clarity using GeneDoc.

Protein sequences from all varieties were compiled into a database, with each serving as its own query sequence. For sequence alignment, Diamond (version 2.1.8) (Buchfink et al., 2021) was utilized, applying an e-value threshold of < 1e-10. The gene annotation files were then converted into GFF format. Duplication classification of genes was conducted using McscanX (Wang et al., 2012), which produced source files detailing gene family duplications. A Perl script, get_mapchart_tandem.pl, was used to generate configuration files for MapChart (Voorrips, 2002), which depicted the genomic locations and tandem duplications of the GST and MYB gene families on chromosomes.

For inter-species collinearity analysis, MCScanX was used to create collinearity files. These files were then modified using another Perl script, highlight.pl, to emphasize blocks containing GST and MYB genes. The visualization of these results was achieved through the JCVI software (Tang et al., 2008), effectively mapping the positions and tandem duplications of the GST and MYB gene families, as well as highlighting significant collinearity blocks.

To elucidate the relationship between MYB genes on Chromosome A07 and structural genes involved in anthocyanin and proanthocyanidin biosynthesis, we employed the String database (https://string-db.org/) for protein-protein interaction (PPI) prediction. The resultant PPI predictions were then visualized through the construction of a network diagram using Cytoscape software (version 3.10.1) (Otasek et al., 2019).

To compare two groups, the Student’s t-test was utilized. For analyzing significance across multiple groups, one-way ANOVA followed by Tukey’s multiple comparisons test was employed. Data are expressed as mean ± standard deviation (SD), based on three independent biological replicates. Significant differences were denoted as follows: *P < 0.05, **P < 0.01, and ***P < 0.001. All graphs and statistical analyses were performed using GraphPad Prism (version 10.1.2).

The synthesis of plant anthocyanins is significantly influenced by transcriptional regulatory factors, primarily MYB, bHLH, and WD40, which together form a regulatory system that controls anthocyanin biosynthesis by binding to structural gene promoters. Once synthesized on the endoplasmic reticulum surface, anthocyanins are transported into the vacuole for storage, facilitated by Glutathione S-transferase (GST) and recognized by multidrug resistance associated protein transporter and ABC protein transporter on the vacuolar membrane ( Figure 1A ). We previously reported a gene, GbBM, which encodes an R2R3 MYB113 transcription factor. This factor directly targets the promoters of four flavonoid biosynthesis genes (GbCHS, GbDFR, GbANS, and GbUFGT), thereby positively regulating the formation of purple spots on cotton flower petals (Abid et al., 2022). In recent years, leveraging map-based cloning techniques, we identified the gene PCC from a red-flowered cotton resource, which governs petal coloration in cotton ( Supplementary Figure 1 and Supplementary Table 1 ). PCC, previously known as GhTT19, encodes a typical glutathione S-transferase. This enzyme is pivotal in petal pigmentation, as it aids in the transport and accumulation of anthocyanins. Significantly, GbBM and PCC/GhTT19 not only share the function of regulating petal coloration but are also positioned closely at the proximal end of the cotton A07 chromosome, separated by a minimal physical distance of merely 0.85MB ( Supplementary Table 2 ). This proximity underscores a potentially significant genetic linkage in the control of floral pigmentation.

To investigate the potential association between the proximal region of the cotton A07 chromosome and pigment development, we examined genes responsible for the coloration of reproductive tissues. Notably, our findings revealed that GhTT2-3A, a gene implicated in the production of brown fibers through the accumulation of anthocyanins/proanthocyanidins, is also situated in this specific chromosomal region. GhTT2-3A, a homolog of the Arabidopsis PA regulator TT2 within the R2R3 MYB TF family, triggers the activation of downstream PA structural genes (Yan et al., 2018). This activation facilitates the synthesis and accumulation of proanthocyanidins (PAs) in cotton fibers, leading to the emergence of brown-colored fibers. The physical separation between GhTT2-3A and GbBM spans 9.2 MB, whereas the gap between GhTT2-3A and PCC/GhTT19 measures 7.6 MB ( Figure 1B and Supplementary Table 2 ). The genes GbBM, GhTT2-3A, and PCC/GhTT19 are closely located on chromosome A07 in cotton, suggesting that the genes responsible for the coloration of cotton reproductive organs likely constitute a metabolic gene cluster at the distal end of chromosome A07.

Further examining the color-determining components, we found that the anthocyanin levels in pcc cotton petals, including pelargonidin, cyanidin, and delphinidin, were significantly higher than in white-petaled cotton Y18R, both during and post-flowering. Notably, the pronounced accumulation of cyanidin in pcc cotton emerged as the key factor for its red flower phenotype ( Figures 1C, D ). Similarly, the presence of cyanidin, delphinidin, and petunidin was instrumental in forming basal spots on the petals of GbBM-induced upland cotton. Additionally, the brown hue of the fibers in brown cotton was predominantly due to the synthesis and buildup of proanthocyanidins.

To analyze the genetic variation of the distal region of chromosome A07 in cotton, we utilized resequencing data from 98 cotton accessions, including 60 G. hirsutum (AD1), 19 G. arboreum (A2), and 19 G. barbadense (AD2). These data were mapped to the reference genome of G. hirsutum (TM-1). We identified a total of 120,126 high-quality SNPs. Population structure analysis and construction of a phylogenetic tree were conducted to understand the genetic diversity within populations. The 98 accessions were classified into four subgroups (G1, G2, G3, and G4) when the value of k was set to 4 ( Figure 2A and Supplementary Table 4 ). G1 subgroup primarily consisted of G. hirsutum with white fiber, G2 subgroup predominantly included G. hirsutum with brown fiber, G3 subgroup mainly represented individuals of G. barbadense, and G4 subgroup was mainly composed of G. arboreum. The results from the phylogenetic tree were consistent with results of population genetic structure analysis ( Figures 2A, B and Supplementary Table 4 ).

The nucleotide diversity (π) of G1 (0.00038) and G2 (0.00049) subgroups was lower compared to G3 (0.0044) and G4 (0.00087) ( Figure 2C ). Both G1 and G2 subgroup belonged to G. hirsutum, while G3 and G4 subgroup were different species. By comparing the Fst values between each group, we observed that the Fst between G1 and G2 subgroup was the lowest, while the differences in Fst between the other two groups were more significant. These results suggest that upland cotton generally has lower genetic diversity and a more stable population. Furthermore, the data indicated that G2 subgroup had higher genetic diversity than G1 subgroup, indicating that brown-fibered upland cotton maintains greater genetic diversity compared to white-fibered upland cotton. Additionally, the results revealed that prolonged artificial selection for long-fibered upland cotton varieties has gradually reduced fiber color diversity, resulting in decreased genetic diversity in modern cultivated white-fibered upland cotton.

Although the Fst difference between G1 and G2 subgroup is relatively low ( Figure 2C ), we still observed some distinct variations in the genome-wide comparison. Notably, there is a significant peak on chromosome A07 ( Figure 2D ). Subsequently, we conducted an analysis specifically focusing on the chromosome A07 and found a distinct difference in the distal region of chromosome A07 ( Figure 2E ). This observation suggests that distal region of chromosome A07 may play a crucial role in determining the brown fibers in cotton. These findings align with previous research that has identified genes responsible for controlling cotton fiber color on the chromosome A07 (Wen et al., 2018).

Gene mapping was conducted on GST genes in 8 different cotton varieties, namely, G. herbaceum (A1), G. arboretum (A2), G. raimondii (D5), G. hirsutum (AD1), G. barbadense (AD2), G.mustelinum (AD3), G. tomentosum (AD4), G. darwinii (AD5), revealing their distribution across all chromosomes in different cotton species. We identified 70, 59, 73, 110, 113, 134, 157, and 149 GST genes in these cotton varieties, respectively. Interestingly, the number of GST gene family members in tetraploid cotton is approximately twice that in diploid cotton ( Supplementary Figure 2A ), it reveals that genome duplication results in a twofold increase in the number of genes. Notably, chromosome Chr03 exhibited the highest number of tandemly duplicated gene clusters for GST genes, while chromosome 06 had the fewest in all diploid cotton species ( Supplementary Figures 3 , 4 ). However, in tetraploid cotton, a significant number of GST genes were found on chromosome A02 and D02, with some forming tandemly duplicated gene clusters ( Supplementary Figures 5 - 9 ). In contrast, chromosome A07 contained only four GST genes in both diploid and tetraploid cotton species, and no gene clusters were observed resulting from tandem duplication of GST genes ( Figure 3A ). The number and location of GST genes on chromosome A07 are consistent throughout the evolution of cotton from diploid to tetraploid, indicating that the gene count remains unchanged during the whole genome duplication (WGD) and tandem duplication events. This result emphasizes the conservative nature of GST genes on the chromosome A07.

Subsequent synteny analysis of the GST gene family revealed a high collinearity of GST genes on chromosome A07 in G. herbaceum, G. hirsutum, and G. barbadense, suggesting their likely origin from G. herbaceum ( Supplementary Figure 19 ). In G. hirsutum and G. barbadense, homologous GST genes on chromosomes A02 and A03 exhibited positional changes compared to G. herbaceum, indicating a replacement of chromosomes A02 and A03 in the A subgenome during the diploid to tetraploid transition ( Figure 3B ). Furthermore, synteny analysis between tetraploid cotton species and G. arboreum with G. herbaceum revealed a higher homology of the A subgenome of all tetraploid cotton species with G. herbaceum ( Supplementary Figure 19 ).

Observations from the constructed phylogenetic tree reveal that Ghe07G08770, Garb_07G021300, Ghir_A07G008080, Gbar_A07G007920, Godar.A07G086400, Gotom.A07G086000, Gomus.A07G083400 cluster together in the Phi subfamily ( Supplementary Figure 21 and Supplementary Table 5 ). Previous studies have indicated that plant GST genes in the Phi subfamily are involved in the transport and accumulation of flavonoids (Kitamura et al., 2012). Sequence alignments demonstrate a high degree of conservation in their sequences, with only two mutation sites identified in G. herbaceum ( Supplementary Figure 22 ), further supporting the notion of high conservation of GST genes related to anthocyanin accumulation and transport across different cotton species. Real-time PCR and RNA-seq analysis to investigate the expression patterns of GST genes in the pcc and Y18R of G. hirsutum during flowering time, reveal that the expression of only the Ghir_A07G008080 (PCC) gene on the chromosome A07 is detected in flowers. Its expression in pcc is significantly higher than in the control cotton Y18R, reaching 90 times the expression level of Y18R ( Figures 3C, D and Supplementary Figure 17 ) The expression of the other three genes is barely detectable ( Supplementary Figure 17C ). It further confirms the key gene of Ghir_A07G008080 influencing the formation of different flower colors in G. hirsutum.

To gain insights into the regulation of reproductive organ color by MYB transcription factors located on the distal end of chromosome A07 in cotton, we conducted a comprehensive analysis. Firstly, we identified the MYB genes present in various cotton species. We respectively identified 216, 211, 218, 416, 426, 430, 438, and 429 MYB genes in the eight cotton varieties ( Supplementary Figure 2B ). Subsequently, we plotted their genomic locations on a chromosomal map, providing a visual representation of their distribution. The distribution of MYB genes across chromosomes in different cotton species was found to be random and uneven. Specifically, in diploid cotton, the MYB genes were predominantly located on chromosomes 05, 07, 08, 11, 12, and 13. Similarly, in tetraploid cotton, the majority of MYB genes in the A subgenome were found on chromosomes A05, A07, A08, A11, A12, and A13 ( Supplementary Figures 10 - 16 ).

We made some interesting observations in our study. Firstly, we found tandem repeat gene clusters on both chromosome 07 in diploid cotton and A07 in tetraploid cotton ( Figure 4A ), except in G. herbaceum, where only one tandem repeat gene cluster was present on chromosome A07. Secondly, during synteny analysis of the MYB gene family, we discovered a high degree of homology between chromosome 07 in G. herbaceum and the A07 chromosome in G. hirsutum and G. barbadense ( Figure 4B ). Next, we observed the expression levels of MYB genes on the A07 chromosome in white and brown cotton at 0, 5, 10, 15, and 20 days post-anthesis (DPA). It was evident that only the Ghir_A07G002090 gene showed significant differences between white and brown cotton fibers. The differences were most pronounced at 5, 10, and 15 DPA, reaching approximately 2-3 times higher expression levels in white cotton fibers ( Supplementary Figure 18 ). Lastly, the synteny analysis revealed that the A07 chromosome remained evolutionarily conserved during the transition from diploid to tetraploid, showing no significant variation ( Supplementary Figure 20 ). This conservation is likely attributed to the importance of maintaining stability in the anthocyanin metabolism process, thereby preventing substantial changes.

Our protein-protein interaction analysis focused on the interactions between MYB proteins on the chromosome A07 and proteins involved in the anthocyanin and proanthocyanin pathways. For example, we observed that GhTT2-3A showed interactions with proteins related to anthocyanin synthesis, such as GhDFR_A06, GhC4H_A10 and GhDFR_D05. ( Figure 4C and Supplementary Table 6 ). This finding further suggests that MYB TF in the distal region of cotton chromosome A07 plays a crucial role in the anthocyanin metabolism process and likely influences the coloration of brown cotton fibers.