The stigmas of hickory were collected from the experimental base in Linglong Mountain, Lin’an, China (119°38′51″E, 30°12′39″N, elevation: 119 m). A total of six hickory trees with similar growth and good vigor were selected, and female flowers were collected every 2–3 days from late April, according to the developmental state. Also, 4 sampling time points were selected, noted as S1, S2, S3, and S4, respectively. The pistillate flowers were taken back to the laboratory immediately after collection. After the sepals were removed, the stigma and ovary parts were separated, placed in liquid nitrogen for freezing, and then stored at −80°C in the refrigerator for pigment determination and RNA extraction for transcriptome sequencing. Each stage (S1 to S4) had three biological replicates.

The extraction and analysis of the total content of chlorophylls and carotenoids were determined according to the method described (Lichtenthaler, 1987; Xue et al., 2020). First, 20 mg of freeze-dried hickory stigma powder was added to 1,000 μl 95%(v/v) ethanol/water solution. After 24 h, it was centrifuged at 5,000 × g for 10 min, measured the absorbance at 652, 665, and 470 nm of the supernatants using an enzyme-labeled instrument (Tecan Spark^®, Swiss). The following equations were used to calculate total chlorophyll and carotenoid contents:

Chl a (μg/g) = (16.29A₆₆₅–8.54A₆₅₂)*10;

Chl b (μg/g) = (30.66A₆₅₂–13.58A₆₆₅)*10;

Total Chl (μg/g) = Chl a + Chlb;

Total carotenoids (μg/g) = (1000A₄₇₀–1.63Chl a–104.96Chl b)÷221*10

The extraction and analysis of the total content of anthocyanin were performed using the pH-differential method (Wrolstad and Giusti, 2001; Xue et al., 2020). First, 4 stages of the stigmas were placed in a lyophilizer (Christ Alpha 2–4 LD plus, Germany) for 48 h, and then grind thoroughly to powder. Next, 25 milligrams of freeze-dried hickory stigma powder were added to 250 μL extraction buffer (Trifluoroacetic acid: formic acid: water: ethanol = 1:2:27:70, v/v), and the mixture was stood 24 h at 4°C in the dark. Following centrifugation at 12,000 rpm for 10 min, the supernatant was transferred to a fresh tube. Next, 100 μL supernatant was mixed with 900 μL.025 mol potassium chloride buffer (pH = 1.0), and another 100 μL supernatant was mixed with 900 μL.4 mol sodium acetate buffer (pH = 4.5). The mixture was stood for 1.5 h (complete sedimentation) at room temperature. The absorption of the mixture was recorded at 510 and 700 nm using the enzyme-labeled instrument with 3 repetitions. Finally, total anthocyanin content was calculated based on the molar absorbance of a cyanidin-3-glucoside standard by the equation:

total anthocyanins (μg/g) = (A × MW × DF × 10000)/(ε × L), whereA = (A₅₁₀−A₇₀₀)_pH1.0−(A₅₁₀−A₇₀₀)_pH4.5, ε = 26900 L⋅mol⁻¹cm⁻¹(extinction coefficient of cyanidin-3-glucoside at 510 nm), L = path length, MW = molecular weight of cyanidin-3-glucoside, DF = dilution factor.

In this experiment, 300 ul anhydrous ethanol was added to a 30 mg sample, sonicated for 3 h, and then macerated overnight. The extract was centrifuged at 8000 rpm for 10 min, and the supernatant was prepared at a concentration of 100 mg/ml (Burits and Bucar, 2000).

The DPPH radical scavenging capacity of hickory female stigma was determined by referring to the method with minor modifications (Burits and Bucar, 2000). The extracts were diluted to 0.25, 0.5, 1, 1.5, and 2 mg/ml with anhydrous ethanol and prepared for use later. Took 500 ul of the diluted extract in a test tube, added 500 ul of.2 mM DPPH free radical ethanol solution, reacted at room temperature, avoided light for 30 min, and then measured the absorbance at 517 nm using an enzyme marker. The absorbance Ac of 500 ul of DPPH free radical ethanol solution mixed with 500 ul of ethanol was measured simultaneously. The IC50 (concentration of samples scavenging 50% of DPPH free radicals) was used to evaluate the magnitude of antioxidant activity of hickory female flower stigma extracts from four periods. The lower the IC50 value, the higher the antioxidant activity.

Clearance rate (%) = (A_c – A_s)/A_c × 100%

The ABTS radical scavenging capacity of hickory female flower stigma was determined according to the method (Petretto et al., 2016). First, ABTS radical stock solution was prepared by mixing 7 mmol/L ABTS solution with 2.45 mmol/L potassium persulfate at a 1:1 ratio by volume and left overnight. The prepared ABTS radical solution was diluted with anhydrous ethanol to produce an absorbance value of 0.7 ± 0.02 at 734 nm. Then, 50 μL of the extract dilutions (0.25, 0.5, 1, 1.5, and 2 mg/mL) were mixed with 1 ml of ABTS radical solution, and the absorbance values were measured at 734 nm after a 30-min reaction at room temperature and protected from light.

Clearance rate (%) = (A_c – A_s)/A_c × 100%

Ac is the absorbance value of anhydrous ethanol mixed with ABTS radical solution, and As is the absorbance value of the extract mixed with the solution with ABTS radical.

Samples were immediately frozen in liquid nitrogen and stored at –80°C for RNA-sequence analysis. Total RNA was extracted from female stigma at different developmental time points using CTAB-PBIOZOL reagent and ethanol precipitation. After extracting, the concentration, 28S/18S and RIN/RQN of total RNA were detected by Agilent 2100, and OD260/280 and OD260/230 were quantified by the NanoDrop. The construction of cDNA libraries and transcriptome sequencing were completed by BGI Technology (Shenzhen, China). First, total RNA was disposed of by mRNA enrichment using Oligo (dT)-attached magnetic beads. The mRNA was further broken using the fragment buffer after ligating sequence adapters. PCR amplification and cyclizing were used to get a single-strand cDNA library. The DNBSEQ platform performed the sequencing after the library was certified. All of the experiments were performed with three replicates. Raw reads were obtained from the sequencing, and the low-quality reads with adaptors were processed to obtain clean reads.

Raw reads obtained were filtered with SOAPnuke (version 1.4) (Li et al., 2008) to filter out low-quality reads, contaminated junctions, and high levels of unknown base N. The clean reads were obtained and stored in FASTQ format, then aligned to the reference genome (hickory genome), using Hierarchical Indexing for Spliced Alignment of Transcripts (HISAT2 version 2.1) (Kim et al., 2015), followed by Bowtie2 (version 2.2.5) to align the clean reads to the gene set (Langmead and Salzberg, 2012). RSEM was used to calculate gene expression levels for individual samples (v1.2.8) (Li and Dewey, 2011). Once the new transcripts are obtained, we added the new transcripts with potential protein-coding to the reference gene sequence to form a complete reference sequence and then calculate gene expression levels. Differential gene expression analysis was performed by DESeq2 (Love et al., 2014) with a Q-value < 0.05. Finally, multiple samples are tested for differentially expressed genes as required, and in-depth clustering and functional enrichment analyses of differentially expressed genes (DEGs) are performed.

The assembled unigenes were functionally annotated using the Kyoto Encyclopedia of Genes and Genomes (KEGG)¹ and Gene Ontology (GO) databases. The Blast2GO program provided GO annotation. WEGO then generated GO classification maps. The biological processes and unigenes annotations of the pathway were analyzed using KEGG. The KEGG comprehensive database resource was compared with our data for analysis. DEGs were identified using the NOISeq method (Tarazona et al., 2015) based on the gene expression level of each sample with a fold change of at least 2, Q-value < 0.01, then were analyzed using RNAseq (3 biological replicates per group), characterized by gene ontology enrichment analysis.

To characterize patterns across samples and identify interesting and highly covariant sets of genes, we constructed a gene co-expression network based on symbolic hybrid network types using the weighted correlation network analysis (WGCNA) package (Langfelder and Horvath, 2008; Linn et al., 2017). The WGCNA package was also used to construct gene co-expression networks from identified DEGs to establish the transcriptional regulatory structure of anthocyanin biosynthesis, carotenoid biosynthesis, chlorophyll metabolism, and ROS metabolism. The network was visualized using Cytoscape software (Smoot et al., 2011).

Totals of 21, 18, and 21 unigenes related to chlorophyll, carotenoid, and anthocyanin metabolism were selected for qRT-PCR analysis. Total RNA was isolated from the stigma collected at four different developmental time points as previously described, and then transcribed 1000 ng of total RNA into cDNA for qRT-PCR using PrimeScript ™ RT reagent Kit with gDNA Eraser (removing genomic DNA) (Perfect Real Time) (TaKaRa, Japan) according to the instructions. The templet cDNA was diluted into 100 ng/μL (± 3 ng/μL). NCBI primer blast was used to design PCR primers to quantify gene expression involved in pigments formation, shown in Supplementary Table 1. The reaction system for RT-qPCR is as follow: ddH₂O 3.8 μL, TB Green 5 μL, primerF.2 μL, primerR.2 μL, cDNA0.8 μL, with 3 repetitions. All qRT-PCR assays were performed using TB Green^® Premix Ex Taq™ (Tli RNaseH Plus) (TaKaRa, Japan) in a CFX96 Touch™ Real-Time PCR Detection System (BIO-RAD, United States) with the following reaction conditions: 95°C for the 30 s and 40 cycles of amplification (95°C for 5 s, 60°C for 30 s). The relative expression levels of target genes were calculated using the 2^–ΔΔCt method against the internal control (Livak and Schmittgen, 2001), and the CcUBC9-5 gene was used as a control to normalize the relative expression levels of target genes. Experiments were performed with three independent biological replicates and three technical replicates.

The protein physicochemical properties, including protein length, isoelectric points (PI), molecular weight (MW), and grandaverage of hydropathicity (GRAVE), were predicted using ProtParam². The ProtScale³, SignalP3⁴, and Plant-mPLoc⁵ were used to predict protein hydrophobic/hydrophilic, signal peptides, and subcellular localization, respectively.

Statistical analysis was carried out with one-way ANOVA in SPSS17 software (ANCOVA; SPSS17, SPSS Inc., Chicago IL, United States). Significant differences among treatments were obtained according to P-values determined by Tukey’s honestly significant difference (Tukey’s HSD) with HSD (P < 0.05). All the experiments were repeated 3 times, and the results were presented as mean ± standard deviation (SD).

During the development of the hickory flower, the stigma color changed obviously from green to dark-purple (S1: green, S2: light-purple red, S3: dark-purple red, S4: dark purple) (Figure 1A). To clarify the pigment components and their dynamic change, total chlorophyll, carotenoid, and anthocyanin at four stages of stigma development were further measured. It was shown that total chlorophyll content was relatively high overall and showed a continuous slight downtrend with stigma development (Figure 1B), while total carotenoid content was increased from S1 to S3 and reached the highest value of 36.4 (ug/g, while decreased gradually from S3 to S4 (Figure 1C). Total anthocyanin content gradually increased during the four periods of stigma development, sharped reaching the highest content (131.4 ug/g, with cyanidin-3-glucoside used as the standard) at the S4 (Figure 1D). These results suggest that the color change of the stigma of hickory flowers is due to the joint action of these pigments. The green color at the early stage is mainly caused by chlorophyll. In the middle and late stages, the accumulation of carotenoid and anthocyanin gradually makes the stigma appear red and purple.

It is well-known that anthocyanins have a strong antioxidant capacity (Raven, 2000; Severo et al., 2015). The sharper changes of anthocyanins contents than the other two pigments from S3 to S4 inspire us to measure the ROS scavenging capacities. To explore the antioxidant capacity of stigma, removal-capacities of DPPH and ABTS radicals were further measured. It was found that the scavenging ability of DPPH radicals gradually increased with the process of stigma development and could reach 84.14% at the S4 period (Figure 1E). Similarly, the change of ABTS free radical scavenging capacity was close to that of DPPH, with 78.46% scavenging capacity at S4 (Figure 1F). It was also observed that DPPH and ABTS radical scavenging rates gradually increased with the increasing concentration of the hickory stigma extracts (Figures 1E,F). DPPH and ABTS IC50 could also be used to describe the oxidative scavenging capacity, the extracts from the stigma samples at S4 were more active in scavenging DPPH radicals, and their IC50s were significantly (P < 0.05) lower than those of the other three periods (S1, S2, and S3) (Figure 1G). The oxidative scavenging capacity obtained by the ABTS IC50 showed similar results. The highest clearance activity was observed during the S4 period and the lowest during the S1 period (Figure 1G).

To obtain a global transcriptome profile of developing stigma, twelve stigma samples at S1, S2, S3, and S4 were sequenced using the BGIseq500 platform. After removing adapter pollution, poly-N sequences, and low-quality reads, totals of 75.72 Gb clean bases were acquired, including an average of 45.24 Mb raw sequencing reads and 42.06 Mb clean reads, in each sample. The Q30 base percentage was greater than 90.78% and the average ratio of clean reads to raw reads was 92.97%. Using HISAT and Bowtie2, 93.43–94.53% of the clean reads could be mapped to the genome of hickory, 56.9–62.31% could be mapped to the genome of hickory uniquely. 67.88–71.17% (total), 52.83–54.87% (unique), and the clean reads were mapped to the hickory genes per sample, so that 31,027 genes were identified (Supplementary Table 2). Principal components analysis (PCA) showed that the sample clusters with high similarity converged, and the 4 sampling time points were separated without an outlier sample (Supplementary Figure 1). Based on the expression levels of all genes in every two samples, the Pearson correlation coefficient (PCC) was calculated and the value was at least 0.94, suggesting that the repeatability of the sample is very good (Figure 2A). Moreover, Hierarchical cluster analysis (HCA) was accompanied and displayed in the form of heatmaps plotted by a heatmap package in R software (Figure 2B), showing that significant differences in the expression levels of genes occurred between the sample groups. Using Mfuzz, twelve time-related gene clusters were obtained based on the similar expression patterns of some genes during stigma development (Figure 3A and Supplementary Table 3). KEGG analysis showed that some genes in cluster 5 were associated with “photosynthesis-antenna proteins,” “porphyrin and chlorophyll metabolism,” “carbon fixation in a photosynthetic organism,” “photosynthesis,” and “carotenoids metabolism,” with gradual downtrends from S1 to S4 (Figure 3B and Supplementary Table 4). These results indicated that photosynthesis mainly dominated in the pre-developmental period. Furthermore, the “Flavonoid biosynthesis” pathway was enriched in cluster 10, with a smooth trend in the pre-period, and increased significantly from S3 to S4 (Figure 3B and Supplementary Table 4). Other flavonoids, such as flavone, flavonol, and isoflavonoid, were found to be enriched in cluster 11, with a drastic reduction tendency from S1 to S2 but a slight change from S2 to S4 (Figure 3B and Supplementary Table 4). Notably, it was also found that oxidative phosphorylation-related genes were clustered in cluster 12, which showed a trend of increasing first from S1to S3 and then decreasing from S3 to S4 (Figure 3B and Supplementary Table 4).

Based on the selected threshold values of fold change with at least 2 (Q-value < 0.01, P-value of the cutoff standard < 0.05), a total of 10,432 DEGs were screened out in all the pair-wise comparison groups. In these DEGs, 2,604 (S2 vs. S1), 2,780 (S3 vs. S2) and 4,955 (S4 vs. S3), 4722(S3 vs. S1), and 6915 (S4 vs. S1) DEGs were counted out, including 1,801, 1026, 3011, 2383 and 3795 upregulated genes and 803, 1,754, 1,944, 2339 and 3120 downregulated genes, respectively (Figure 4A). In addition, Venn diagrams were mapped to figure out the relationship between the DEGs in the above comparison groups, indicating that a total of 185 DEGs existed in all five comparison groups (Figure 4B). These results suggested that there were dynamic changes in identified genes during stigma development and the active transcription occurred during the late period. KEGG analysis indicated that all identified DEGs were enriched in 15 pathways (Supplementary Table 5). It was also found that these DEGs enriched in the “Flavonoid biosynthesis” pathway (S4 vs. S1, S4 vs. S3, ko00941), “Photosynthesis-antenna proteins” (S3 vs. S1, S3 vs. S2, and S4 vs. S1, ko00196) pathways, “Photosynthesis” (S4 vs. S1, ko00195) pathway, “anthocyanin biosynthesis pathway” (S3 vs. S1, and S4 vs. S1, ko00942), “carotenoid biosynthesis pathways” (S4 vs. S1, ko00906) and ABC transporters (S4 vs. S3, ko02010) (Figure 4C and Supplementary Table 5). Furthermore, volcanic maps in S3 vs. S1, S3 vs. S2, S4 vs. S1, and S4 vs. S3 were drawn, intuitively showing the up-regulated genes (blue dots), down-regulated genes (red genes), and non-differential genes (green dots) (Figures 4D-G). These results indicated that the genetic manifestations of these pathways were active at the middle and late stages of female stigma development, offering a clear direction for elucidating the internal molecule mechanisms of the color formation and dynamic changes.

The above results indicated that chlorophyll, carotenoid, and flavonoid, including anthocyanins, may be involved in color formation. According to the whole genomic information of the hickory and KEGG database, the metabolism pathways of chlorophyll, carotenoid, and flavonoid, including anthocyanins, were drawn and genes encoding enzymes and transcription factors were annotated in these pathways (Figures 5–7 and Supplementary Table 6). Flower color-related genes with critical roles were further analyzed and 91 genes were screened out. Among these genes, 27 genes are required for chlorophyll metabolism, 28 genes were related to carotenoid metabolism and 36 genes were associated with anthocyanins metabolism. Among these genes, 21, 19, 22, and 18 from S1, S2, S3, and S4, respectively, had FPKM values of over 100. In the chlorophyll metabolic pathways, some genes, such as GSA, HEMB, HEMD, HEMF1/2, CLH1/2, FCC, and NYE1, were upregulated, while HEMA, HEMC, HEME, HEMG, CHLD, CHLH, CHLI, CRD, DVR, POR, CAO, CHLG, NYC1, and PAO, were downregulated (Figure 5 and Supplementary Table 6). In the carotenoid biosynthesis pathway, some genes, including DXS, MCT, CMK, MDS, HDS, GGPPS, and the downregulated genes DXR, PSY, PDS, ZISO, ZDS, LCY-E, LCY-B, LUT5, LUT1, ATVDE, and ZEP, were shown to be upregulated (Figure 6 and Supplementary Table 6).

Based on the results showing that there were greater changes in anthocyanin level and strong antioxidant capacity, a detailed analysis of the expression patterns of anthocyanin-related genes was conducted according to the data in Supplementary Table 6. Through the analysis of the synthesis pathway, we found that the expression patterns of most genes, including PAL, C4H, 4CL1/5, 4CL2, 4CL3, CHS, F3H, F3′H, F3′5 ′H, DFR, ANS, UF3GTs, UGT, UGT73B1/B2/B3/B4/5, and UGT88A1, were basically consistent, showing upward trends from S3 to S4, contributing to rapid accumulation of anthocyanins at a late stage of stigma development. It was also observed that the FLS gene required for quercetin-Der, Kaempferol-Der, and myricetin-Der biosynthesis showed the opposite trends as the above-described genes, suggesting that more premise substances enter into the leucocyanidin, leucopelargonidin, and leucodelphin for anthocyanins accumulation. After the anthocyanins biosynthesis is finished, they will be transported into the vacuole, involving some transporter proteins. It could be seen that the expression trends of TT2 coding transporter proteins were similar to that of anthocyanins biosynthesis genes, while others, including TT19, MRPs, and AHA10, showed the opposite trends from S3 to S4. This is an indication that TT2 may be the main transporter for anthocyanins in hickory stigma. In addition, previous studies have reported that the anthocyanin metabolism pathway was regulated by a suite of transcription factors that include MYB, bHLH, and WD-repeat proteins (Walker et al., 1999; Nesi et al., 2000; Gonzalez et al., 2008; Jaakola, 2013). It was observed that some TFs, including TT2, TT1, and TTG1, showed upward trends from S3 to S4, while others like MYB75, MYB11, MYB12, MYC1, TT8, and GL3 were on the contrary. This may be attributed to the fact that these TFs usually form protein complexes to function. To regulate the expression of anthocyanin-related genes, these TFs need to be transcribed earlier. When the data were homogenized via Z-score normalization, the trends in gene expression became clearer (Supplementary Figure 3). Together, the expression trends of these pigment-related genes provided plausible explanations for the dynamic change of hickory stigma color with development.

To validate the results of RNA-seq, 20, 16, and 21 unigenes related to chlorophyll metabolism, carotenoid biosynthesis, and anthocyanin biosynthesis were selected for qRT-PCR analysis. As was shown in Figures 5–7, the results of the qRT-PCR were similar to the RNA-seq analysis, showing a high correlation, except for minor differences in some genes. These results further confirmed the reliability of the RNA-seq data in this study.

To gain insight into the patterns of genetic association between different samples and further study the high covariate gene sets associated with pigment and ROS metabolism, we used WGCNA to construct a co-expressed gene network based on 10,432 DEGs (Supplementary Table 7). After removing the low expression genes, nine different gene modules were obtained by HCA and were demonstrated in different colors. Among nine modules, four significant modules were found through analyzing the relationship between modules and traits (Figure 8A). The darkred-, brown2-, firebrick4- and dark turquoise-modules were significantly correlated with the trait at S1 (PCC = 0.92, P = 2*e^–5), S2 (PCC = 0.99, P = 8*e^–10), S3 (PCC = 0.97, P = 2*e^–7), and S4 (PCC = 1, P = 3*e^–12), respectively. Figure 8B further showed the module eigengene in the four modules, which indicates these four modules were strongly interrelated in the four stages (Figure 8B). This result illustrated that the DEGs in the darkred, brown2, firebrick4, and the dark turquoise module were more important to correlate with the 4 stages during hickory stigma development. In the 4 hub modules, 17 chlorophyll-related DEGs and 12 carotenoid-related DEGs were mostly contained in the darkred module, while most anthocyanin biosynthesis-related DEGs were distributed in the dark turquoise module, with 112 DEGs (Supplementary Table 8). GO analysis further showed the “photosynthesis” and “photosystem” were enriched in the darkred module, which was closely related to S1 (Figure 8C and Supplementary Table 9). Notably, “oxidoreductase activity, acting on single donors with incorporation of molecular oxygen, GO:0016701,” “oxidoreductase activity, acting on single donors with incorporation of molecular oxygen, incorporation of two atoms of oxygen, GO:0016702,” and “active transmembrane transporter activity, GO:0022804” related pathway were concentrated on the brown2 module (Figure 8C and Supplementary Table 9). These related DEGs in the firebrick4 module were associated with “hormone-binding, GO:0042562,” “ethylene receptor activity, GO: 0038199,” “ethylene binding, GO:0051740” and “NAD(P)H dehydrogenase (quinone) activity, GO:0003955” (Figure 8C and Supplementary Table 9). In the dark turquoise module, these DEGs were involved in “oxidoreductase activity, acting on diphenols and related substances as donors, oxygen as acceptor, GO:0016682,” “phenylpropanoid metabolic process, GO:0009698” and hydroquinone: oxygen oxidoreductase, GO:0052716” (Figure 8C and Supplementary Table 9). These results further suggested that pigment and ROS metabolism were very active during hickory stigma development.

According to pigment and ROS metabolism as well as related transcription factors in 4 key modules, a co-expression network was further constructed (Figure 8D). It was found that this network included 23 chlorophyll-related genes, 19 carotenoid-related genes, 141 anthocyanin-related genes, 46 ROS-related genes, and 623 transcription factors (TFs) (Supplementary Table 10). These transcription factors belonged to different families, including MYB (99), AP2(52), bHLH (52), NAC(30), WRKY(28), C₂H₂(21), HB-HD-ZIP(20), bZIP9(19), Trihelix(18), and others (284) (Supplementary Figure 2). Notably, 7 transcription factors (6 in the MYB family and 1 in the bHLH family) were shown to be related to anthocyanin metabolism. Moreover, the number of connections between genes was also counted, ranging from 49 to 447 (Figure 8E and Supplementary Table 10). In the chlorophyll metabolism pathway, HEME2 (CCA1503S0015) and CAO (CCA1377S0020) had relatively high connections with more than 400 connections. It was also found that the connections of only the PSY (CCA0887S0096) gene in carotenoid biosynthesis-related genes were over 400. In the anthocyanins metabolism pathway, genes with more than 400 connections were FLS1 (CCA0006S0007), 4CL8 (CCA0881S0026), F3′H (CCA0743S0031), F3′H (CCA0507S0050), and TT12 (TRANSPARENT TESTA 12-like, CCA1606S0005). In ROS metabolism, GSTU18 (glutathione S-transferase TAU 18, CCA0747S0058) was shown to have over 400 connections. Among genes encoding all the transcription factors, there were 21 genes with more than 400 connections, of which, the gene expression difference of WRKY75 (WRKY transcription factor 75, CCA0573S0068) was the largest during stigma development. Remarkably, MYB113 (MYB transcription factor 113, CCA0887S0030) was more connected to other genes in 7 anthocyanin-related transcription factors. Based on this, we further mapped the co-expression network of these two transcription factors and genes related to pigment and ROS metabolism, which more clearly showed the regulatory relationship of transcription factors on them (Figure 8E). Using ProtParam, the physicochemical properties of the two transcription factors were further analyzed, as shown in Table 1. It was demonstrated that the protein length, pI, MW, and GRAVY of WRKY75 and MYB113 were 152 and 264 aa, 9.64 and 9.1, 17.64 and 30.21 kD, and –1.156 and –0.536, respectively. Both proteins were hydrophilic with no signal peptides. Furthermore, the subcellular localization prediction showed that they were both located in the nucleus. Taken together, these results provide in-depth guidance for revealing the transcriptional regulation of pigment and ROS metabolism in the future.