Co-regulatory network analysis of the main secondary metabolite (SM) biosynthesis in Crocus sativus L.

Saffron (Crocus sativus L.) is being embraced as the most important medicinal plant and the commercial source of saffron spice. Despite the beneficial economic and medicinal properties of saffron, the regulatory mechanism of the correlation of TFs and genes related to the biosynthesis of the apocarotenoids pathway is less obvious. Realizing these regulatory hierarchies of gene expression networks related to secondary metabolites production events is the main challenge owing to the complex and extensive interactions between the genetic behaviors. Recently, high throughput expression data have been highly feasible for constructing co-regulation networks to reveal the regulated processes and identifying novel candidate hub genes in response to complex processes of the biosynthesis of secondary metabolites. Herein, we performed Weighted Gene Co-expression Network Analysis (WGCNA), a systems biology method, to identify 11 regulated modules and hub TFs related to secondary metabolites. Three specialized modules were found in the apocarotenoids pathway. Several hub TFs were identified in notable modules, including MADS, C2H2, ERF, bZIP, HD-ZIP, and zinc finger protein MYB and HB, which were potentially associated with apocarotenoid biosynthesis. Furthermore, the expression levels of six hub TFs and six co-regulated genes of apocarotenoids were validated with RT-qPCR. The results confirmed that hub TFs specially MADS, C2H2, and ERF had a high correlation (P < 0.05) and a positive effect on genes under their control in apocarotenoid biosynthesis (CCD2, GLT2, and ADH) among different C. sativus ecotypes in which the metabolite contents were assayed. Promoter analysis of the co-expressed genes of the modules involved in apocarotenoids biosynthesis pathway suggested that not only are the genes co-expressed, but also share common regulatory motifs specially related to hub TFs of each module and that they may describe their common regulation. The result can be used to engineer valuable secondary metabolites of C. sativus by manipulating the hub regulatory TFs.

www.nature.com/scientificreports/co-regulatory network of those genes was constructed.Before carrying out the network analysis, t-SNE analysis was applied to all data sets and assessed patterns and trends of clustering among samples.This analysis explained the local similarities among the data set.t-SNE revealed that the majority of the homologous stigma transcriptomic projects of different ecotypes were grouped into the same category (Supplementary Fig. S1).We used the normalization methods on gene expression data to remove the samples and features with missing information and adjust for sources of noise.We clustered reference batch based on t-SNE analysis.Overall, four batches were found.Removing batches and reducing the unwanted variations improved the identification of outliers and error correction in this study.The results demonstrated that combat-based normalization appeared to give consistent results for building a network.The weighted gene co-expression network analysis (WGCNA) approach is applied for mining specific modules and highly correlated genes related to the secondary metabolites.In order to discover hub genes and co-regulatory networks which play main roles in various secondary metabolites of C. sativus stigmas, 19 samples were analyzed by WGCNA.The optimal soft threshold was determined at 5 to construct a scale-free network (Supplementary Fig. S2).The adjacency matrix was computed using the function adjacency, and the topological overlap matrix was established based on dissimilarity between genes by measuring the distTOM = 1-TOM function.Clustered genes were divided into modules, and modules were obtained based on the dynamic cutting tree branches algorithm (Supplementary Fig. S3).Module eigengenes were preform to merge modules whose dissimilarity was below the 0.2 cutoffs (Fig. 1).A total of 11 modules related to the metabolic biosynthesis pathway were gained (Supplementary Fig. S4).The module size ranged from 37 purple to 137 turquoise modules.In order to explore the notable modules that play the regulatory role in the biosynthesis of secondary metabolic pathways in C. sativus, KEGG enrichment analysis (www.kegg.jp/ kegg/ kegg1.html) 44 was performed in the genes on the interest modules, and considerable enrichment pathways were screened.
The results demonstrated that in the stigma of C. sativus: the "Apocarotenoid, Carotenoid biosynthesis, MEP Pathway, Phenylpropanoid biosynthesis, ABC transports" were significantly enriched in brown module."MVA Pathway, Saponin biosynthesis, monoterpene, Apocarotenoid" were enriched in the blue module.In the turquoise module, the most enrichment pathways were "Flavone and flavanol biosynthesis"."Flavonoid biosynthesis" and "Apocarotenoid" were enriched in pink and green modules, respectively.Interestingly, numerous genes related to Apocarotenoid biosynthesis were the most significant in brown, blue, and green modules.In addition, the list of the enrichment pathways in the over-representative modules is presented in Table 1.

Hub TF recognition and visualization
In order to identify hub TFs associated with metabolites biosynthesis, hub gene analysis was performed for genes in the modules.Hub gene analysis recognized transcription factors including MADS, C2H2, ERF, bZIP, HD-ZIP, and zinc finger protein in the brown module with the most degree and MCC.bZIP and MYB in the green module were also recognized as hub TFs.The blue module contained HD-ZIP, and HB GATA (Fig. 2 and Supplementary Table S2), and the turquoise module included MIKC and HB hub TFs.MYB-related and NF-YB were also highly co-expressed in the pink module.The results of the PPI network indicated the existence of strong relationships and significant connections among most hub genes.The constructed PPI networks based on the hub TFs of three modules related to apocarotenoid biosynthesis indicated ideal and significant connectivity, which emphasized the effectiveness of our approach to organizing functional modules that comprised a set of proteins having similar functions.Hence, these modules might influence the regulation of secondary metabolites biosynthesis.Therefore, warrants further validation for our findings (Supplementary Fig. S5).

Promoter analysis
Genes with similar regulatory motifs could be co-regulated.One of the best approaches for the description of co-regulated genes is the determination of common cis-acting elements located on their promoter.The commonly shared motifs were derived as promoter analysis by MEME.According to transcription factor binding sites (TFBSs) screening, multiple unique TFBSs of hub TFs in the promoter regions of co-expressed genes in modules related to apocarotenoids were detected (Table 2, Supplementary Table S5).The cis-elements related to

Validation of the relationship between hub TFs and related genes with RT-qPCR
To verify the chief conclusion drawn from the co-regulation network results for stigma samples, the relative expression levels of hub TFs and candidate genes involved in apocarotenoid synthesis including MADS, C2H2, ERF, CCD2, ADH11367, GLT2 in the brown module, HB and ADH3F1 in the blue module and bZIP, MYB, UGT12 and UGT73D1 in the green module between high-content metabolite ecotypes ("Ghaen", "Arjenak" and "Zaveh") and low-content metabolite ecotypes ("Shahr-e Kord", "Torbat" and "Hamedan") were determined (Fig. 4).In total, the RT-qPCR results were consistent with the results of the coregulation network analysis.Several research studies have demonstrated a close relationship between metabolite content and mRNA expression level 24,25,[48][49][50][51][52] ).The RT-qPCR analysis has shown that the relative expression of regulatory TFs and genes tended to give higher up-regulation levels in the high-content ecotypes group than in the low-content ecotypes.
As shown in Supplementary Fig. S6, the correlation between hub TFs and genes was positively significant.The relative expression of UGT73D1 which is involved in the synthesis of picrocrocin 53 , in "Ghaen" as a high-content metabolite ecotype was fourfold higher than "Shahr-e Kord" and the correlation between hub TF MYB and UGT73D1 was more than 0.9.The gene expression of GLT2, which is involved in the conversion of crocetin to crocin 19 increased 8.1 fold in "Ghaen" compared with "Shahr-e Kord" and there was also a significant correlation between hub TFs (MASS, C2H2, ERF) and GLT2 gene (R 2 > 0.86).Herein, the partial least squares (PLS) regression result illustrated that MADS and C2H2 positively affect related-genes in the apocarotenoid pathway.The PLS result showed that there is also a strong correlation between the C2H2 and GLT2, and ERF and ADH367 (Fig. 5).It is recommended that the potential participation of a regulatory network exists between hub TFs and related genes in the metabolite biosynthesis pathway among different ecotypes.We also performed PLS analysis to determine the effect of regulon hub genes on their targets in nine transcriptome sequencing data samples from three developmental stages (RED, -2 DAY, and 0 DAY) of stigmas in C. sativus which their apocarotenoid contents were also assessed (low and high amounts of apocarotenoids, respectively) 54 .The results showed the strong effect of hub TFs such as MADS, C2H2, ERF, and HD-ZIP on related apocarotenoid genes especially CCD2 and GLT2 (Supplementary Fig. S7).Interestingly, the result of our analysis was nearly consistent with the above results which confirmed the role of the identified hub TFs in the regulation of apocarotenoid biosynthesis in C. sativus.

Discussion
C. sativus is a chief source of crucial medicinal metabolites such as apocarotenoids crocin, crocetin, picrocrocin and safranal.Most of the studies in this plant, just have examined metabolite contents assay and gene expression without recognizing relationships among genes and considering regulatory factors including transcription factors.Whereas evaluating the interactions among all regulatory factors can obtain an immense insight into the regulatory mechanisms involved in the secondary metabolite pathways.In this study, we used RNA-Seq data of C. sativus stigma to establish a co-regulatory network for genes involved in the secondary metabolite biosynthesis pathway by using WGCNA analysis.The achievement of integrated regulatory networks can help to understand a more comprehensive knowledge of molecular mechanisms involved in secondary metabolite synthesis.We identified 11 modules related to the metabolite pathways in C. sativus and three noticeable modules, namely brown, blue and green modules, which were related to the apocarotenoids synthesis pathway.It is notable that the high connectivity among genes in the same modules can be attributed to the close expression of some genes in the module.In order to discover these potential genes that may be related to the regulation of the biosynthesis of secondary metabolites, different approaches were performed, including assigning modules, functional enrichment analysis of modules, and hub-hub genes recognition.Our findings presented evidence that the regulated hub TF genes may be associated with the coordinated expression of genes involved in the biosynthesis of secondary metabolites.The outcomes also designated that one module can concurrently be related to multiple components, and one component can be regulated by multiple modules.
Hub TFs are considered highly correlated genes for a given module in a biological interaction network.Hub genes were significantly connected in terms of the PPI network.Therefore, TFs with a high degree of connectivity within secondary metabolite modules were also associated at the protein-protein interaction level.In the present study, transcription factors MADS, C2H2, ERF, HB, MYB, and bZIP were recognized as hub genes in modules related to apocarotenoid biosynthesis as showed co-regulated expression with their associated genes in the network.Moreover, the high correlation and positive effect were demonstrated among the genes involved in apocarotenoid production and the hub regulatory TFs in both high-and low-content metabolite ecotypes.PLS analysis showed a considerable effect of C2H2 and MADS on the related genes of apocarotenoid, especially CCD2 and GLT2.Previous research reported that MADS, MYB, MYB-related, WRKY, C2C2-YABBY, and bHLH contain the necessary regulatory machinery for apocarotenoid biosynthesis 40,55 .The TFs including zinc-finger motifs have been previously recognized due to their potential biological functions related to the regulation of apocarotenoid biosynthesis 56 .The screen of MADS, ERF, C2H2, MYB, and zinc-finger protein transcription factors found that these straightly regulate carotenoid genes to definitely and coordinately modulate carotenoid metabolite in plants such as apple fruit, 'Benin Shogun' , 'Yanfu 3' fruit flesh [56][57][58][59] .The roles of three of these mentioned hub TFs have been determined and confirmed in previous research in C. sativus.For example, the study conducted by Malik et al. 60 highlighted the role of CstHD, whose expression pattern corresponds to apocarotenoid accumulation in C. sativus stigmas.Overexpression of CstHD led to an increase in apocarotenoid content by upregulating biosynthetic pathway genes.Transient expression of CsHD which is identified as hub TFs in our study of co-regulation network, raised apocarotenoid content in metabolism.Furthermore, zinc-finger transcription was identified as a potential regulator of apocarotenoid biosynthesis, with AN20/AN1 being nuclear localized and activating reporter gene transcription in yeast, suggesting its involvement in regulating apocarotenoid biosynthesis 56 .MYB, the main regulatory factor, and genes involved in the apocarotenoid pathway were found in the same modules.
Transient overexpression experiments involving MYB genes in C. sativus have further confirmed their regulated role in apocarotenoid metabolism.Specifically, these TFs were revealed to regulate apocarotenoid biosynthesis by interacting with the promoters of genes involved in the secondary metabolite pathway 61 .Our findings revealed that more effective and crucial hub TFs such as MADS, C2H2, ERF, HB, and bZIP could be important regulatory genes for regulating the apocarotenoid biosynthesis pathways.
Promoter analysis is a powerful method to describe the presence of common cis-acting elements among genes of the co-expression modules.MYB, WRKY and MADS TFs commonly regulate different genes in the same module which are involved in the metabolite biosynthesis pathway by recruiting to the same cis-acting elements that are placed on their promoters [62][63][64][65] .Presence of a few and the same cis-regulatory elements in the promoter region of genes in the same module highlights this point, that genes that participate in the same pathway have a co-regulatory network 66,67 .The results confirmed that "Ghaen", "Arjenak", and "Zaveh" not only had the highest expression in RNA transcript level of the hub TFs and related genes but also produced the highest amount of the apocarotenoids among the studied ecotypes.The expression level of the apocarotenoid genes in three ecotypes, "Shahr-e Kord", "Torbat" and "Hamedan" was notably less than the other ecotypes.A study conducted by Tan et .al2019 indicated that red and 0 DAY stigmas of the flowering stage, which illustrated the low and high amounts of apocarotenoids, respectively, showed low and high expression of MADS, MYB, zing finger TFs, and UGT2 gene, respectively.A study of the effect of abscisic acid treatment in saffron demonstrated that crocin and safranal contents and the gene expression of CsGT2 and CsLYC were increased coordinately 68 .It seems that a positive correlation would be between metabolite contents and expression profiles in medical plants 24,25,[48][49][50][51][52] .A thorough understanding of their mechanism will be the aim of future research.In addition, five genes involved in ABCC Transporters (MATE1, extrusionprotein1b sub-familyCmember4b, Cmember4c, and sub-familyCmember4c) were identified in the brown module.Previous studies were shown ABCC transporters mediated the vacuolar accumulation of crocins in saffron stigmas and were co-expressed with the gene encoding the first dedicated enzyme in the crocin biosynthetic pathway, CsCCD2 69 .This gene had highly coordinated expression and correlation with some of the TFs such as MADS, C2H2, ERF, bZIP, HD-ZIP, and zinc finger protein in the brown module.
Glycosyltransferases play various roles in cellular metabolism by modifying the activities of regulatory metabolites.UGTs catalyze the connections between glucose molecules and specific receptors by glycosidic bonds.Crocetin glucosyltransferase 2 (GLT2) is a glucosyltransferase that plays important role in the biosynthesis of crocins.Herein, GLT2 was more highly expressed in high-content metabolite than in low-content metabolite ecotypes.This gene indicated the highest correlation with associated TFs relative to other genes involved in apocarotenoid biosynthesis based on the test for statistical significance of the correlation coefficient.Several other studies also showed that GLT2 was expressed at high levels in stigmas which facilitate sequential glucosylation of crocetin to crocin in C. sativus 19,70,71 .The UGT85 family as stress-regulated UDP-glucosyltransferase-encoding genes have been co-expression with genes related to metabolite biosynthesis in the brown module.The UGT85 family that comprises members related to stress responses and cell cycle regulation in C. sativus 72 .In this study, the turquoise module contained the genes of the flavonols pathway that were highly co-expressed.The expression of main flavonols genes including CsGT45 (UGT75P1), and UGT703B1 is highly correlated with the presence of kaempferol 73 .The results illustrated that a close relationship has existed between the expression of CsZCD and stress-regulated genes in the relevant modules based on the increased coordinated expression levels of these genes in stress conditions 74 .So, these findings can predict stress influence on apocarotenoid synthesis by modifying the gene expression and oxidative and peroxidative reactions.Further explorations are clearly required to uncover the relationship between stress conditions and the biosynthesis of apocarotenoids.

Conclusion
C. sativus is the major source of crucial medicinal apocarotenoid metabolites.Herein, we performed a co-regulatory network analysis of the transcriptome projects of C. sativus.We identified 11 modules related to C. sativus secondary metabolites which three specialized modules covered functions in the apocarotenoid biosynthesis pathway.The result revealed that hub TFs such as MADS, C2H2, ERF, bZIP, MYB, and HB were involved in regulating apocarotenoid biosynthesis.The relative expression of genes involved in apocarotenoid production and associated TFs revealed a strong correlation among various ecotypes, which represented different apocarotenoid

RNA-Seq data analysis
Before constructing an integrated regulatory network, data preparation was performed.All processes were run on the Linux operating system (Ubuntu 20.04 LTS).The quality assessment of raw RNA-Seq data was performed by using the FastQC tool (version 0.11.9) 75.To obtain high-quality reads and remove ambiguous base (N) and adaptor sequences, raw reads were pre-processed by fastp tool (version 0.20.1) 76.
The clean reads were used for all the downstream analyses.De Bruijn graph-based assembler Trinity package with 32 k-mer was used for de novo assembly 77 .The mapping was performed by using Bowtie 2 (version 2.4.1).To annotate assembled transcripts, we used BLASTx with a cut-off of E-value ≤ 1e − 5 against the NR database, the GO database, and the Swiss-Prot protein database 78 , and to identify transcription factors (TFs) in the C. sativus transcriptome data, the PlantTFDB v5.0 database was used 79 .Also, genes and TFs identification along with considering the relevant literature review was performed 40,55,80 .Salmon (version 1.9.0) 81 was used to quantify the expression of each transcript, and a gene expression matrix was generated based on Transcripts Per Million (TPM) as well as applied the Log2 transformation method to diminution residual variability.

Batch effects detection and normalization
Batch effects and unwanted variation in RNA-seq experiments are triggered by various factors.Evaluation and removal of batch effects from data could enhance prediction performance in gene expression data.It is essential to normalize the larger merged dataset.We performed tSNE (Rtsne package version 0.16) to recognize the pattern of clustering among different samples.ComBat (SVA package version 3.44.0)function in surrogate variable analysis (SVA) was applied to estimate batch.We used Mean-only adjustment, the popular ComBat model, which adjusts for batch effects and allows the operator to specify the batch name or number to be used as the reference batch 82 .

Weighted coregulation network analysis and module detection
Regulatory networks and modules of genes were constructed by the WGCNA package in R language 30 .In order to ensure the reliability of the constructed network, the outlier samples were checked and omitted.For this aim, adjacency matrices of expression matrix were considered and sample network connectivity was standardized based on the distances.GoodSamplesGenes function in the WGCNA package was applied to exclude the unqualified genes.The argument maxPOutliers was considered as 0.05 to specify outliers.
To construct scale-free network, a suitable soft-thresholding power was computed based on scale-free topology.The scale-free network is defined with power-law degree distribution and includes several nodes with few interactions and few nodes with high interactions, which are named hub nodes.The approximate scale-free soft threshold was chosen as β = 5.The power 5 exhibited the lowest possible power level at which the scale-free topology fit index curve plateaued, achieving a high value of 0.8, while also attaining a high mean connectivity value.based on this fact, the mean of connectivity falls as power goes up.For module detection, step by step modules function in the WGCNA package was performed with the following main parameters: power = 5, cor-Type = "bicor", networkType = "Signed Hybrid", TOMType = "signed", minModuleSize = 30, deepSplit = 2.The adjacency matrix was converted into a topological overlap matrix (TOM), and the corresponding dissimilarity was measured.The identification of modules was carried out by average linkage hierarchical clustering by the topological overlap matrix of samples.The dynamic tree cut algorithm was used to detect clusters in the hierarchical dendrogram, then module eigengenes were applied to merge highly similar modules.

Detection of hub genes
The assignment of intramodular connectivity or association between a gene within a specific module and other genes in that module was defined by module membership (kME).MEDissThres was considered as 0.20 and kME threshold was considered > 0.8.Genes with high module membership (kME > 0.8) were good representatives of the overall expression profile in the module and genes with high module membership tend to be "hub" genes in the module (high intramodule connectivity) in our study.To prevent the phenomenon of fuzzy clustering, we considered high KME values.The applied threshold ensures the selection of genes with strong module membership, indicating their importance in the network's structure and function.Also, the genes with the highest association were known as hub genes.CytoHubba plugin of Cytoscape software (version 3.9.1) 83was applied to obtain hub genes in each module using Maximal Clique Centrality (MCC) and degree methods for topological analysis.STRING online tool (Version: 12.0) (https:// string-db.org/) 84 was used to determine the possible protein interactions network for the hub genes.The background file was chosen monocot genome annotation information

Plant cultivation
The corms of saffron were collected from different geographical locations of traditional saffron production areas in Iran (Fig. 6 and Supplementary Table S4).The taxonomy of the C. sativus plant was confirmed by Dr Hamzeh-Ali Shirmardi, from Chaharmahal and Bakhtiari Agricultural and Natural Resources Research and Education Center, Shahrekord, Iran (shirmardi1355@gmail.com).The plant corms and the voucher specimen (Her-barium No.1400402/381) were stored in the collection of medicinal plants of the Department of Agricultural Biotechnology at the Tarbiat Modares University.Different ecotypes of C. sativus were grown in a climate-controlled greenhouse under a 10 h photoperiod in Shahr-e Kord, Iran, with coordinates 32.3282° N, 50.8769°E, at the temperature of 14-18 °C and decent loamy soil (20% sand, 30% silt, and 30% clay) with pH 7.8-8.3from September to December 2021.In order to minimize fungal diseases, triggered chiefly by Fusarium spp, 1% fungicide water solution of copper oxychloride was used to treat the corms for 30-60 s.The biological replicates of saffron stigmas were randomly harvested in the early morning hours during the flowering stage (on November 11th, 2021, Supplementary Fig. S8).At this stage, some of the morphological traits such as fresh stigma weight, dried stigma weight, stigma length, corm weight, horizontal diameter, flower fresh weight, and day -to-flower were measured in the studied ecotypes, the result presented in Supplementary Table S3.The fresh stigmas were frozen in liquid nitrogen and preserved at -80 °C for metabolite analysis and RNA extraction.

High performance liquid chromatography (HPLC) Analysis
The stigma tissues were freeze-dried with a freeze dryer machine (CHRiST model ALPHA 1-4) and then were stored at 4 •C until further analysis.Chemical standards of picrocrocin (CAS number 138-55-6, Sigma-Aldrich,  S4. (> 90% purity) were used.10 mg of stigma tissues were grounded and then homogenized in 10 ml ethanol at 80%.Then samples were sonicated by an ultrasonicator (Stima sonic Company-Iran) for over 15 min in dark conditions at 25 °C.Next, the samples were centrifuged for 5 min at 8000 rpm and filtered following the supernatant was used for chromatography by HPLC/DAD to determine saffron components 88 .HPLC equipment was performed by KNAUER System equipped with a binary well chrome K1001, a multiple wavelength UV-Vis (DAD)-2800 model KNAUER Eurospher equipped with a 20 μl injection loop and RP C18 (4.6 × 250 mm × 5 µm) the column was used for analytical separation.A combination of acetonitrile (solvent A) and water (solvent B) was applied as a mobile phase with a flow rate of 1.0 ml/min at room temperature.The gradient method was used at 90% for 5 min then decreased to 20% in 20 min and was kept for more than 5 min.The acquisition of the chromatograms was adjusted at different maximum absorption including 308 nm, and 440 nm to determine picrocrocin, and crocins, respectively.The identification of each compound was achieved by spiking its retention time with standards under equal conditions.The quantification analysis was assessed by evaluating the integrated peak area and concentration levels, which was considered using the calibration curve by plotting the peak area against the different concentration levels of the particular standard compounds.Fourpoint regression curves were obtained in the concentrations of 2.5,10, 50, and 200 ppm for the crocin standard and 25,50, 100, and 200 ppm for the picrocrocin standard.The correlation coefficients (R 2 ) were measured for all authentic standards of compounds, which were about 0.99.

Validation of selected genes of the modules with qRT-PCR
Based on network analysis results, we selected several identified hub TFs and associated genes from the noticeable apocarotenoid modules for validation by RT-qPCR. 100 mg red stigma tissue in all bulked stigma samples (15 stigmas for each replicate) was used for RNA extraction.Total mRNA extraction was performed using DENAzist Column RNA Isolation Kit (DENAzist Asia Co., S-1010, Mashhad, Iran) and treated with DNase I (Thermo Fisher Scientific, Waltham, MA, USA) to remove remaining genomic DNA contamination according to the manufacturer's instructions.The RNA quality and integrity were checked using 1.2% agarose gel, and the concentration was examined using Nanodrop (BioTek, EPOCH).Approximately 1 μg of DNase I-treated total RNA was subjected to the cDNA synthesis using Add Scrip cDNA synthesis kit (add-bio, Korea, Lot:22,701) based on the manufacturer's instructions.Primers were designed by PerlPrimer (version 1.1.21)software 89 and Oligo-Analyzer (Version 3.1) web-based program (https:// www.idtdna.com/ pages/ tools/ oligo analy zer) subsequently, were checked using in silico NCBI Primer-BLAST tool (https:// www.ncbi.nlm.nih.gov/ tools/ primer-blast/) to

Figure 1 .
Figure 1.Average linkage hierarchical clustering by the Topological Overlap Matrix and adjacency-based dissimilarity applied to recognize modules.the original modules along with the assigned merged modules are shown with different colors in the below the dendrogram.

Figure 3 .
Figure 3. HPLC analysis of the accumulation of crocins and picrocrocin in different ecotypes of C. sativus.The bubble graph shows the contents of crocins (A) and picrocrocin (B).The size of each bubble indicates the mean of three biological replicates.

Figure 5 .
Figure 5. Partial Least Squares regression (PLS) was carried out for determining complex relationships between TFs and dependent genes.The results indicate a strong correlation and more effect between the hub TFs (C2H2 and MADS) and related genes (GLT2 and CCD2) and between ERF and ADH11367.PLS analysis in low-(A), high-(B), content metabolite and total (C) ecotypes.

Table 1 .
Weighted Gene Co-Expression Network Analysis (WGCNA) for significantly expressed genes in saffron stigma and pathway enrichment analysis of main modules in Crocus sativus.The brown, blue and green modules are the specific modules of apocarotenoids.P values are computed from Fisher's exact test (α < 0.05).

Figure 2. Co-regulated network of hub TFs and structure genes related to secondary metabolites biosynthesis pathway in C. stativus. Each TF and the related edge are demonstrated with a specific color. Yellow circles indicate apocarotenoid genes. The notable apocarotenoid modules including blue module
(A), brown module (B), and green module (C) are shown, respectively.

Table 2 . Consensus motifs identified in promoter region of genes of modules involved in apocarotenoids biosynthesis from C.sativus.
Vol:.(1234567890) Scientific Reports | (2024) 14:15839 | https://doi.org/10.1038/s41598-024-65870-z . Overall, the result of the co-regulatory network analysis of C. sativus provided insights into the specific regulatory mechanisms of apocarotenoid biosynthesis at the gene expression level and further provided some candidate hub TFs for regulon engineering.

Table 3 .
The list of primers for hub TFs and related genes used in quantitative Real-Time PCR expression analysis.