To investigate the features of the CYP75 gene family in Orchidaceae, 13 orchids with completed whole-genome sequencing were selected. They include C. goeringii (Chung et al., 2021; Sun et al., 2021), C. sinense (Yang et al., 2021), C. ensifolium (Ai et al., 2021), Gastrodia elata (Yuan et al., 2018; Xu et al., 2021), D. catenatum (Zhang et al., 2016; Niu et al., 2021), D. chrysotoxum (Zhang et al., 2021b), D. huoshanense (Han et al., 2020), Phalaenopsis aphrodite (Chao et al., 2018), P. equestris (Cai et al., 2015) of Epidendroideae, Platanthera guangdongensis (Li et al., 2022), Pl. zijinensis (Li et al., 2022) of Orchidoideae, Vanilla planifolia (Hasing et al., 2020) of Vanilloideae, Apostasia shenzhenica (Zhang et al., 2017) of Apostasioideae. The genome download urls for 13 orchids can be found in Table S1 . And four CYP75A (AUB13331.1 from Horgeum vulgare; ABI95365.1 and AAZ79451.1 from D. hybrid; AEB96145.1 from D. moniliforme) and four CYP75B (BAJ93256.1 from H. vulgare; AT5G07990.1 from A. thaliana; XP015613041.1 from Oryza sativa; AF155332.1 from Petunia hybrida) proteins were downloaded from the NCBI (https://www.ncbi.nlm.nih.gov/). All the protein sequences can be found in the Table S2 .

It is complicated to identify CYPs due to their extremely high sequence conservation (Mizutani and Ohta, 2010; Babu et al., 2013; Kumar et al., 2014; Reddy et al., 2014). Therefore, eight CYP75 proteins were used as queries to perform a blast search (built-in TBtools; Chen et al., 2018) against 13 orchid genomes with an E-value threshold of e^-20 and a requirement for amino acid sequence similarity greater than 50% (Falginella et al., 2010). Pfam verification was performed on all protein sequences in the blast results. The CYP domains (PF00067) built on the hidden Markov model were downloaded from Sanger center (http://pfam.sanger.ac.uk/). The identification of all alternative orchid CYPs was conducted using the Hmmsearch program (built-in Tbtools; Chen et al., 2018). The subsequent screening study included all protein sequences containing CYP domains. To facilitate identification, 14 A. thaliana CYP proteins (AtCYP71A12, AtCYP71A13, AtCYP71B15, AtCYP73A5, AtCYP76C1, AtCYP77B1, AtCYP77A6, AtCYP78A10, AtCYP82G1, AtCYP84A1, AtCYP89A2, AtCYP98A3, AtCYP701A3, AtCYP706A1) belonging to subfamilies other than CYP75 were also downloaded from the TAIR (https://www.arabidopsis.org/). Then, ML phylogenetic trees were constructed for further screening of CYP75 genes based on all CYP genes of each orchid and other species. The phylogenetic analysis was conducted using the maximum likelihood (ML) approach, and the ML tree was constructed using the RAxML on the CIPRES Science Gateway web server (RAxML-HPC2 on XSEDE; Miller et al., 2015) with 1,000 bootstrap iterations. The phylogenetic tree of CYPs for each orchid is detailed in the Figure S1 . The target genes were clustered with eight CYP75 query genes. Ultimately, we aligned all candidate CYP75 protein sequences of 13 orchids and further eliminated the incomplete gene annotation sequences ( Figure 1 ). The completed protein sequences of orchid CYP75s can be found in Table S2 . The physicochemical properties of CYP75 proteins were predicted by ExPASy database (https://www.expasy.org/) (Artimo et al., 2012). Subcellular localization was predicted by Plant-mPloc (http://www.csbio.sjtu.edu.cn/bioinf/plant-multi/#) (Chou and Shen, 2010).

We performed multiple sequence alignment of CYP75 proteins from 13 orchids and other species using MEGA 7.0 software (Kumar et al., 2016). The alignment sequences selected with the ClustalW program, Gap Opening and Gap Extend, are 15 and 6.66, respectively; the DNA Weight Matrix selection is the IUB; other values keep the default. The phylogenetic analysis was conducted using the maximum likelihood (ML) approach, and the ML tree was constructed using the RAxML on the CIPRES Science Gateway web server (RAxML-HPC2 on XSEDE; Miller et al., 2015) with 1,000 bootstrap iterations.The best model (JTT+I+G4) was determined by using the Akaike Information Criterion (AIC) and the Bayesian Information Criterion (BIC) based on modeltest-ng-0.1.3 (Darriba et al., 2019). Keep other settings as default. The output phylogenetic tree file was polished using Evolview (http://www.evolgenius.info/evolview/) (He et al., 2016).

To illustrate the variations of motifs among orchid CYP75 proteins, the MEME motif search tool was employed to detect the conserved motifs in all CYP75 proteins (Bailey et al., 2009). The MEME parameters were optimized as follows: the maximum number of motifs found was set to 20, and the optimal motif width ranged from six to 50 bases. GSDS6 (http://gsds.gao-lab.org/) (Hu et al., 2015) was used for analyzing structure. The protein motifs and gene structures of CYP75s were visualized with TBtools (Chen et al., 2018).

Utilizing chromosome-level genome assemblies of C. goeringii, C. ensifolium, and D. chrysotoxum, genomic FASTA files were merged pairwise to generate a database for BLASTp queries. The merged BLAST files and modified GFF3 files for each species were analyzed using MCscanX (Wang et al., 2012) to identify collinear blocks of CYP75 genes between C. goeringii and C. ensifolium, C. goeringii and D. chrysotoxum. The dual synteny plotter tool of MCscanX (JCVI kit) was utilized for visualization of the results pertaining to collinearity.

We selected six orchids with chromosome-level genomes to analyze the distribution of CYP75 genes on their chromosomes. According to the genome data and annotation data of C. goeringii, C. ensifolium, D. huoshanense, D. chrysotoxum, Pl. guangdongensis and Pl. zijinensis, the chromosomal localization of five orchid CYP75 genes were visualized through TBtools (Chen et al., 2018).

A total of 2,000bp upstream and downstream of CgCYP75s was extracted via TBtools, respectively (Chen et al., 2018). The online software PlantCARE (http://bioinformatics.psb.ugent.be/webtools/plantcare/html/; Lescot, 2002) was used to identify and annotate the cis-acting elements found in the upstream and downstream regions. Cis-acting element number and response function were visualized using TBtools (Chen et al., 2018).

For transcriptome analysis, RSEM (Li and Dewey, 2011) was utilized to quantify transcripts and compute the fragment per kilobase of transcript per million mapped reads (FPKM) score for every gene. Expression heatmaps were produced using TBtools with the FPKM matrix (Chen et al., 2018). To confirm the expression patterns of the CYP75s, sepals, petals, labellums (lips), and gynostemiums were sampled from yellow-green C. goeringii (‘YG’) and purple-red C. goeringii (‘PR’), respectively, these were grown at Fujian Agriculture and Forestry University for quantitative real-time PCR (qRT-PCR) experiment. Each tissue type was sampled in three replicates. Total RNA of these tissues was extracted using the FastPure Plant Total RNA Isolation Kit (Vazyme Biotech Co., Ltd., Nanjing, China). First-strand DNA was synthesized with TransScript® All-in-One First-Strand cDNA Synthesis SuperMix for quantitative PCR (TransGen Biotech, Beijing, China). Premier 5 software was used to design primers for candidate genes and internal reference genes for qRT-PCR (Zhao et al., 2022). Gene-specific primers for two chosen genes and their corresponding internal control genes are presented in Table S3 . The qRT-PCR was performed to verify the specific expression of CgCYP75A1 and CgCYP75B1 in the floral organs of ‘YG’ and ‘PR’ C. goeringii. All experiments were conducted in triplicate, with each run consisting of three technical replicates. The relative expression of genes was determined using the 2^-ΔΔCT method (Wang et al., 2021; Zhao et al., 2023).

EggNOG-mapper v2 (http://eggnog-mapper.embl.de/) was used to perform a search against the eggNOG5.0 database for gene ontology (GO) functional annotation (Huerta-Cepas et al., 2019). The prediction of orthology was conducted through sequence alignment, while applying bit-score or E-value filtering to improve the quality of orthology assignments. Functional classification was attained by associating the GO annotation terms with the proteins involved in established biological processes. The results of the GO analysis are available in the Table S7 and visualized using Tbtools (Chen et al., 2018).

A total of 72 CYP75s were identified from 13 orchids, with the number of CYP75s in each orchid ranging from two to ten ( Table 1 ). Our study reveals that CYP75B subfamily has a significantly higher number (53/72) compared to CYP75A subfamily (19/72). Each of the 13 orchids contains one−four CYP75A genes and between two−six CYP75B genes. Among them, V. planifolia has the largest number of CYP75 genes, with a total of ten (four CYP75A and six CYP75B), while Pl. guangdongensis has the smallest number with only two CYP75B genes and lacks CYP75A gene.

The 72 CYP75 protein sequences range from 235−668 amino acids, with a mean of 493. The molecular weight ranges from 26.10−76.21 kDa, with a mean of 54.80 kDa. Around 88.89% (64/72) of the CYP75 proteins have high isoelectric points (pI>7), with an average of 7.57. The average instability index (II) is 41.39, and 39 CYP75 proteins are below this index, indicating good protein stability. The average aliphatic index (AI) for the 72 CYP75 proteins is 99.40, indicating high thermal stability. Moreover, the calculated mean hydrophilic index (GRAVY) of CYP75 proteins in all orchids is negative, indicating a high degree of hydrophilicity. All CYP75 proteins are localized within the endoplasmic reticulum (ER), as evidenced by subcellular localization results. This localization pattern is consistent with the majority of CYP proteins, which are known to primarily function in the ER (Neve and Ingelman-Sundberg, 2010).

Gene ontology analysis was performed to delineate gene functional classifications of orchid CYP75s and investigate the important biological processes they might be involved in. As a result, GO terms “response to stimulus,” “response to auxin”, “response to organic substance”, “response to endogenous stimulus”, “response to chemical”, “response to hormone”, “secondary metabolic process”, “secondary metabolite biosynthetic process”, “obsolete oxidation-reduction process”, “membrane” and “oxidoreductase activity” constituted the greatest number of genes for GO ontologies “Biological Process”, “Cellular Component,” and “Molecular Function”, respectively ( Figure S2 ; Table S7 ). The data suggests that the gene ontology of orchid CYP75s is significantly enriched in the “Biological Process”, which is strongly linked to the production of plant metabolites. Additionally, it is heavily concentrated within response elements that relate to the plant’s reaction to external environmental factors.

A phylogenetic tree was constructed to analyze the evolutionary patterns of orchid CYP75s ( Figure 2 ). 72 CYP75 proteins were used, and eight query protein sequences were used. The maximum likelihood (ML) method was used for the phylogenetic tree, which was constructed with RAxML on the CIPRES Science Gateway web server (RAxML-HPC2 on XSEDE; Miller et al., 2015). Bootstrap values were 1,000 replicates with the JTT+I+G4 model. The phylogenetic tree indicated that CYP75s belonged to two categories: CYP75A and CYP75B, which is consistent with the previous studies (Tanaka and Brugliera, 2013). The results clearly show that the number of members in the CYP75B subfamily has far exceeded than that of the CYP75A subfamily in the course of evolution. The branch relationships of the eight known CYP75 protein sequences are found to be correct.

Motifs of CYP75 proteins in 13 orchids were examined using the online analysis tool MEME, and 20 motifs were set as upper bound ( Figure 3A ). There are four conserved motifs commonly found in cytochrome P450 enzymes: PERF motif, K-helix region, I-helix region and the consensus sequences of the heme-binding region ( Figures 1 , 3B ), also referred to as the “P450 signature” (Crooks et al., 2004; Qi et al., 2017). As depicted in Figure 1 , all CYP75 proteins found in orchids possess the four common CYP motifs referred to earlier. A total of 20 motifs were detected in the CYP75 proteins using the MEME software (Bailey et al., 2009). The number of CYP75 motifs ranges from nine to 18. The results show that most CYP75 proteins have identical sequence beginning with motif 18 and followed by motif 5, 13, 2, 14, 10, 6, 8, 16, 12, 3, 1, 17, 11, 4, 7, 15, ending with motif 9. Furthermore, all protein sequences of CYP75 exhibit the highly conserved motif 12, 3, 1, 17, 11, 4, 7, 15, and 9. ( Figure 3A ). Among the 20 motifs, motif 1 corresponds to K-helix region, motif 4 to heme-binding region, motif 11 to PERF motif, and motif 12 to I-helix region ( Figure 3B ). Furthermore, the unique conserved motifs in each subfamily are also found: motifs 19 are specific to the CYP75A subfamily, while motifs 17 are exclusively present in the CYP75B subfamily. This finding will greatly facilitate the identification of subfamilies within the CYP75 gene family. These 72 CYP75s share many common motifs, indicating high conservation. Most of the CYP75 protein sequences are conserved with differences mainly in the N terminus and intermediate regions. However, there are some CYP75 protein sequences containing a small number of motifs but also contain the typical motifs of CYP protein.

To further explore the characteristics of CYP75 genes in orchids, intron-exon structure is analyzed as shown by Figure S3 . The results show that the orchid CYP75 family is composed of one−eight exons and one−seven introns, and in DchCYP75B3 and DhuCYP75B4, the exons are split into many small fragments by introns. Most CYP75 protein sequences have relatively long intron regions, while all CYP75A genes contain only one intron, which is a unique feature of the CYP75A subfamily.

As shown in Figure 4A , both CYP75A genes of C. goeringii are located on chromosome 02, while the remaining CgCYP75B genes are scattered on chromosomes 01, 08, 11, and 14. Among them, chromosomes 02, 08 and 11 all have two CgCYP75s located at the same site. The four CYP75s of C. ensifolium are scattered on chromosomes 01, 02, 09 and 11, respectively ( Figure 4B ). The CYP75 genes of D. huoshanense are evenly distributed across five chromosomes, specifically chromosomes 01, 06, 10, 16, and 18. The two DhuCYP75A genes are located together on chromosome 06 ( Figure 4C ). In D. chrysotoxum, only DchCYP75B4 and DchCYP75B5 are co-located on chromosome 07, while the remaining DchCYP75 are separately distributed on various chromosomes. Among them, DchCYP75A1 has not been assembled onto the chromosome, but is located on scaffold 787 ( Figure 4D ). For both Platanthera species, two PguCYP75Bs are co-located on chromosome 02 and 08 of Pl. guangdongensis, while three PziCYP75s are individually located on chromosome 01, 10 and 13 of Pl. zijinensis ( Figures 4E, F ).

To investigate the evolution of CYP75 genes in orchids, the collinear relationship among CYP75s in C. goeringii, C. ensifolium, and D. chrysotoxum was analyzed. Our collinear analysis revealed a one-to-one correspondence among all CYP75 genes in the three orchids, indicating limited reshuffling of CYP75 orthologs and significant genomic rearrangements following the divergence of Dendrobium and Cymbidium lineages ( Figure 5 ). Furthermore, we also examined that CYP75 gene tandem duplication occurred on the chromosomes of both C. goeringii and D. chrysotoxum, which was consistent with the results of chromosome localization ( Figure 4 ). Our results suggest that a small-scale tandem duplication may have led to the expansion of CYP75 gene family in orchids.

To explore the regulatory roles of CYP75s, we retrieved the 2,000 bp upstream and downstream regions of CYP75 genes in C. goeringii to identify potential cis-elements. We identified a total of 1,257 cis-acting elements, including 38 types and ten responsive functions ( Figure 6 ; Table S4 ). Among these elements, TATA-box made up the most common elements (46.38%), followed by CAAT-box (22.83%) ( Table S5 ). Cis-element functions included phytohormone responsiveness for gibberellin, auxin, methyl jasmonate (MeJA), salicylic acid, and abscisic acid (ABA); stress responses, such as anoxic, anaerobic, low-temperature, and defense; and growth and development elements, such as light response and MYB binding site ( Figure 6 ). Each CgCYP75 gene contained multiple types of elements with light responsiveness as the most occurring element function ( Figure 6 ), supporting that light is one of the most important environmental factors affecting flavonoid biosynthesis in plants (Zoratti et al., 2014). The second and third most abundant types of elements identified were MeJA-responsive and ABA-responsive elements ( Table S4 ). The results suggest that these elements may play a role in modulating these two phytohormones. Of particular note are the MYB binding site elements, which are key transcription factors involved in pigment synthesis. The significant proportion of these elements further emphasizes the critical role of the CYP75 genes in the synthesis of plant pigments.

The orchid transcriptomic data from nine orchids out of 13 were visualized as heatmaps via TBtools software as shown in Figure 7 . Expression analysis was conducted using transcriptome data from nine orchids, including various plant organs, such as sepals, petals, labellums (also known as lips), the gynostemium, inflorescence, stems, tubers, pollinia, and complete flowers. The expression profile showed that CYP75B genes were expressed broadly in flower organs, while the expression of CYP75A genes was more confined to C. ensifolium, D. catenatum and G. elata, particularly in the gynostemiums ( Figure 7 ). In C. goeringii, CYP75A genes showed little expression in all tissues of two varieties, while CgCYP75B1 exhibited an exclusive expression in four organs of ‘PR’ C. goeringii and gynostemium of ‘GY’ C. goeringii ( Figure 7A ). Furthermore, CeCYP75B1 and CeCYP75B2 exhibited elevated expression levels across all four floral organs ( Figure 7B ). In D. chrysotoxum, the expression of the DchCYP75B4 gene was significantly higher on the lip than other parts, due to the predominant distribution of its anthocyanin in the red macula of the lip ( Figure 7C ). A similar situation existed in D. catenatum, DcaCYP75B1 and DcaCYP75B3 were highly expressed in the lip with purple-red spots ( Figure 7D ). PeqCYP75B1 was expressed prominently in various tissues of P. equestris, with higher expression in the darker colored lip region. On the other hand, PeqCYP75B5 showed significant expression levels in the sepals. ( Figure 7E ). AshCYP75B3 of A. shenzhenica was expressed in all parts except for the pollinium, with the highest expression in the inflorescence ( Figure 7F ). In both of the Platanthera species, PziCYP75B1 and PguCYP75B4 exhibited the highest expression levels in fleshy underground tubers ( Figure 7G ), possibly indicating their involvement in non-biological stress response (Iwashina, 2003; Pourcel et al., 2007; Zhang et al., 2020). Similarly, G. elata was a mycoheterotrophy orchid without anthocyanins in its tissues, and the expression of GelCYP75A and GelCYP75B were comparable in various tissues ( Figure 7H ). Their functions are not related to anthocyanin synthesis but are possibly associated with physiological activities required for adaptation to its unique habitat (Iwashina, 2003; Pourcel et al., 2007; Zhang et al., 2020).

To further investigate the specific roles of CYP75 gene expression in two C. goeringii varieties, their sepals, petals, lips and gynostemiums were analyzed by qRT-PCR ( Figure 8 ). In two C. goeringii varieties ( Figure 8A ), the CgCYP75B1 showed extremely high expression in ‘PR’ C. goeringii, especially in petal which matched the transcriptome data, but were barely detected in ‘GY’ C. goeringii, further verifying that CYP75B genes have an anthocyanin-specific expression in the flower organs. As for CgCYP75A1, its expression level was relatively low in four types of tissues of the two C. goeringii plants ( Figure 8B ). Interestingly, CgCYP75A1 and CgCYP75B1 showed a certain degree of stable expression in the gynostemiums of both C. goeringii plants ( Figure 8B ). The function of these CYP75 genes may not be related to anthocyanin synthesis, but rather to attracting pollinators to promote reproduction (Samanta et al., 2011).