Plant CPR sequences from 20 plant species were obtained from the web-based resource for Arabidopsis P450, cytochrome b5, NADPH-cytochrome P450 reductases, and family 1 glycosyltransferases¹, NCBI², and each plant genome database. Accession numbers of CPR genes and genome databases used in this study are listed in Supplementary Table 1. Multiple sequences were aligned using ClustalW and were used for tree construction using MEGA7 with the neighbor-joining method. The position of protein helices was predicted using an in silico transmembrane helix prediction tool³ (Supplementary Figure 1). The conserved membrane anchor, FMN-, FAD-, CYP-, and NADPH-binding domains were predicted as described by Qu et al. (2015). Motif analysis of the FMN domain of legume CPR classes I and II was performed using the MEME Suite online software⁴ using a minimum of 37 sequences for each CPR class from 24 different legume species (Supplementary Figure 2). The structure model of MtCPRs was constructed using 5gxu.1.A (crystal structure of Arabidopsis ATR2) as a template with >71% identity in the SWISS-MODEL (Arnold et al., 2006). PyMol (DeLano, 2020) was used to visualize the amino acids at the active site of the enzyme.

β-amyrin, erythrodiol, oleanolic acid, α-amyrin, uvaol, ursolic acid, and lupeol were purchased from Extrasynthese (Genay, France). Betulin, methyl jasmonate, and salicylic acid were purchased from Sigma-Aldrich (St. Louis, MO, United States). Betulinic acid and gibberellin were purchased from Tokyo Chemical Industry (Tokyo, Japan). Sophoradiol, 24-hydroxy-β-amyrin, 24-carboxy-β-amyrin, 11-oxo-β-amyrin, 30-hydroxy-β-amyrin, and 11-deoxoglycyrrhetinic acid were kindly gifted by Dr. Kiyoshi Ohyama (Tokyo Institute of Technology).

Methods of G. uralensis sample preparations, RNA extraction, and library construction for de novo assembly have been described in Supplementary Methods.

Co-expression analysis of M. truncatula and L. japonicus was performed using an online transcriptomic database from the gene expression atlas web servers https://lipm-browsers.toulouse.inra.fr/pub/expressionAtlas/app/mtgeav3 and https://lotus.au.dk/expat/ (Miyakojima MG20 v 3.0), respectively. Transcriptomic data for M. truncatula, L. japonicus, and G. uralensis were obtained from 274, 74, and 37 different samples, respectively. Pearson’s correlation coefficients (PCCs) of different CPRs and CYPs were calculated using Excel (Microsoft Corp., Redmond, WA, United States).

The position of protein helix inside, transmembrane, and outside of the endoplasmic reticulum (ER) were predicted using in silico transmembrane helix prediction (see foot note 3) (Supplementary Figure 1). In this study, the N-terminal domain are represented by the protein located inside ER together with the transmembrane helix. To swap the N-terminal domain of MtCPR1 and MtCPR2, the first 46 amino acids of MtCPR1 are swapped with the first 60 amino acids of MtCPR2 using NEBuilder^® HiFi DNA Assembly cloning to produce N-terminal sequence of MtCPR1 fused with truncated-N-terminal MtCPR2 (M1N-M2C) and N-terminal sequence of MtCPR2 fused with truncated-N-terminal MtCPR1 (M2N-M1C). Two different fragments of each designed chimeric CPRs (M2N-M1C and M1N-M2C) were amplified using pENTR-MtCPR1 and pENTR-MtCPR2 as template and primers listed in Supplementary Table 2. A minimum of 20-bp overlapping nucleotides of each CPR mutant and pENTR™ as vector backbone were designed for fusion cloning.

Cytochrome P450 reductases class I and II genes from M. truncatula, L. japonicus, and G. uralensis were amplified from the respective cDNA using PCR with PrimeSTAR^® Max DNA Polymerase (Takara Bio, Kyoto, Japan) and named MtCPR1, MtCPR2, LjCPR1, LjCPR2, GuCPR1, and GuCPR2. The primers used in this study are listed in Supplementary Table 2. Wild-type and chimeric CPR genes were cloned into the pENTR™ vector using the D-TOPO^® (Thermo Fisher Scientific, Waltham, MA, United States) or NEBuilder^® HiFi DNA Assembly (Ipswich, MA, United States) cloning method, to produce entry clones of CPR genes. Yeast expression clones of CPR, using pAG415-GAL-ccdB (plasmid number 14145, Addgene, Watertown, MA, United States) as the destination vector, were constructed using LR reaction with LR Clonase™ II Enzyme mix (Thermo Fisher Scientific), to produce pAG415-CPR. pAG415GAL-ccdB was a gift from Susan Lindquist. Yeast expression clones of MtCYP716A12, MtCYP72A63, LjCYP716A51, LjCYP93E1, LjCYP72A61, and GuCYP88D6 were constructed using LR reaction with LR Clonase™ II Enzyme mix (Thermo Fisher Scientific) into pYES-DEST52 (Invitrogen, Carlsbad, CA, United States) and a Gateway™-compatible version of pESC-HIS generated previously in our laboratory (Yasumoto et al., 2016) as the destination vectors, to produce pYES-DEST52-CYP and pESC-HIS-CYP (Seki et al., 2008, 2011; Fukushima et al., 2011; Suzuki et al., 2019).

Two different yeast strains were used in this research, INVSc1 (MATa his3D1 leu2 trp1-289 ura3-52; Thermo Fisher Scientific) and PSIII strain. PSIII strain was generated by transforming pESC-TRP[P_GAL10/tHMG1-T2A-upc2-1][P_GAL1/LjbAS]) as described in Supplementary Methods 2 into PSI strain (BY4742/P_erg7:P_MET3-ERG7/trp1:P_ACT1-Gal4dbd-ER-VP16) constructed in a previous study (Srisawat et al., 2020). INVSc1 strain carrying pYES3-ADH-aAS, pYES3-ADH-OSC1, or pYES3-ADH-LUS (Fukushima et al., 2011) and PSIII strain were then co-transformed with pAG415-CPR, pYES-DEST52-CYP, and pESC-HIS-CYP consecutively using Frozen-EX Yeast Transformation II™ (Zymo Research, Irvine, CA, United States).

In case of INVSc1, yeast sample with cell density at an OD₆₀₀ value in the range of 2.2–2.5 (late logarithmic phase) were cultured in 5 mL appropriate synthetic defined medium (Clontech, Palo Alto, CA) containing 2% glucose at 30°C and 220 rpm for 24 h. Yeast cells were collected by centrifugation, resuspended in 5 ml of appropriate synthetic defined medium (Clontech) containing 2% galactose and then incubated at 30°C and 220 rpm for 48 h. In case of the PSIII strain, yeast strains with cell density at an OD₆₀₀ value in the range of 1.4–1.6 (mid-logarithmic phase) were cultured in 5 mL appropriate synthetic defined medium (Clontech) containing 2% glucose, 100 nM β-estradiol, and 1 mM L-methionine and then incubated at 30°C and 220 rpm for 48 h. For the feeding assay, INVSc1 yeast was cultured as described previously and supplemented with 10 μM of erythrodiol at the same time as galactose addition, and yeast were then cultured at 30°C and 220 rpm for 60 h.

Before extraction, 20 ppm of ursolic acid or uvaol were added to 5 mL yeast culture as internal standards. Triterpenoid metabolite extraction was performed as previously reported (Fanani et al., 2019). The dried extract was dissolved in 500 μL of chloroform:methanol 1:1 (v/v). Upon Gas chromatography–mass spectrometry (GC-MS) analysis, 100 μL of the sample was evaporated to dryness and trimethylated with 50 μL of N-methyl-N-(trimethylsilyl) trifluoro acetamide (Sigma-Aldrich) for 30 min at 80°C.

GC-MS was performed as previously reported (Fanani et al., 2019). For peak identification, authentic standards were used to confirm the peak in the sample and compare the retention time and mass fragmentation patterns. The relative amounts of triterpenoids in each strain were compared by calculating the area of each peak from the extracted chromatogram based on m/z values of 203 (Kim et al., 2018). Each sample of the different CYP-CPR combinations was analyzed using three biological replicates.

The significance of difference between two yeast strains harboring different CPR classes from same species was determined by single-factor analysis of variance (ANOVA). One-way ANOVA was used to test the difference of triterpenoid levels among different yeast strains and the significance for the means were separated by Tukey’s test. p-values less than 0.05 (P < 0.05) were considered to be significant for this study. All the statistical analyses were performed with SPSS 20.0 for Windows (IBM, Chicago, IL, United States).

The phylogenetic analysis results from 20 different plant species showed clear branching of CPR class I and II sequences (Figure 2). Each dicot plant species was found to have a minimum of one CPR sequence in each class. Monocots, gymnosperms, and mosses possess several copies of CPR-II genes, but no CPRs belonging to CPR class I have been found to date (Jensen and Møller, 2010). Amino acid sequences of plant CPRs from different plant families were compared to create a table of similarity matrices (Supplementary Table 3). Same CPR classes of the same plant family were found to have 75–95% similarity, while same CPR classes of different plant families were found to have 64–75% similarity. On the other hand, different CPR classes of the same or different plant families were found to have 55–65% similarity. Motif analysis was performed on different plant CPR classes (Supplementary Figure 2).

Multiple sequence alignments from several plant CPRs showed that CPR-I has approximately 10–20 amino acids shorter N-terminal membrane sequences than CPR-II (Figure 3). In this study, we focused on the structural analysis of the FMN domain of legume CPRs using CPR from a model legume, MtCPR. Based on motif analysis of amino acids from CPR classes I and II from 24 legume species, there were some differences in the amino acid residues between CPR classes I and II in the FMN domain, which are conserved in each CPR class (Supplementary Figure 2). By multiple sequence alignment with hCPR, we highlighted eight amino acid residues of MtCPR1 and MtCPR2, which are co-localized with the position of hCPR mutations and negatively charged amino acid residues in the FMN domain. These amino acids have been reported to play a role in CPR:CYP-specific interactions (Campelo et al., 2018; Esteves et al., 2020; Figure 3). The location of each highlighted amino acid in the protein structure was observed through structural analysis of MtCPRs (Figure 4).

We selected two model legumes, Medicago truncatula and Lotus japonicus, and one of the high-value medicinal legumes, Glycyrrhiza uralensis, to analyze the gene expression levels of different CPR classes in legumes. Available transcriptomic data of M. truncatula and L. japonicus were retrieved from the respective transcriptome database, while transcriptome data of G. uralensis were obtained in this study (Supplementary Tables 4.1–4.3). RNA-seq analysis from 37 data-sets of G. uralensis resulted in 351,138,706 reads and were assembled into 226,599 unigenes (Supplementary Table 5). Co-expression analysis of CPR-I and CPR-II genes of M. truncatula, L. japonicus, and G. uralensis showed distinct expression patterns (Figure 5). CPR-Is were found to be constitutively expressed with lower and more stable expression levels (Figure 5). In contrast, CPR-IIs were generally found to have higher expression levels than CPR-Is, which varied depending on tissues and treatments (Figure 5). Overall, CPR-I showed higher expression levels than that of CPR-II in roots or nodules, while CPR-II showed higher expression levels than of CPR-I in chemically or biologically treated samples (Supplementary Table 6). Global analysis of the genes from each legume species was performed by calculating the PCC of each CPR class with each of the transcripts found in the assembly (Supplementary Table 7). CPR-I showed strong correlation (PCC > 0.70) with genes involved in basal metabolism or molecular functions, such as ATP-binding/-carrier proteins (Supplementary Table 7). In contrast, CPR-II correlated strongly with receptor genes such as the ubiquitin, glutamate, and jasmonate receptors, which are known to play a role in plant defense responses (Guo et al., 2018; Qiu et al., 2020; Doroodian and Hua, 2021) (Supplementary Table 7). Interestingly, despite its critical role in providing electrons for CYPs, the global transcriptomic analysis of genes correlated with CPRs showed no correlation with CYPs, at a PCC > 0.4 (Supplementary Table 7).

In order to analyze the correlation strength between the different CPR classes and CYPs in each species, we calculated the PCC values of different CPR classes in combination with CYP families 716A, 72A, and 93E from the three legume species, in addition to CYP88D6 from G. uralensis. MtCYP716A12, LjCYP716A51, and GuCYP716A179 catalyze the same reaction of three steps of oxidation at the C-28 position of α-amyrin, β-amyrin, and lupeol to produce the final product of ursolic acid, oleanolic acid, and betulinic acid, respectively, (Fukushima et al., 2011; Tamura et al., 2017; Suzuki et al., 2019; Figure 1). Both MtCYP72A63 and GuCYP72A154 oxidize the C-30 position of β-amyrin to produce 30-OH β-amyrin, but only MtCYP72A63 catalyzes the further oxidation to produce 30-COOH β-amyrin (Seki et al., 2011). LjCYP72A61 catalyzes the one-step oxidation at the C-22 position of β-amyrin to produce sophoradiol (Ebizuka et al., 2011). MtCYP93E2, LjCYP93E1, and GuCYP93E3 catalyze the same reaction of three steps of oxidation on the C-24 position of β-amyrin to produce the final product of 24-COOH β-amyrin (Seki et al., 2008; Fukushima et al., 2011; Suzuki et al., 2019). GuCYP88D6 catalyzes the two-step oxidation at the C-11 position of β-amyrin to produce 11-oxo β-amyrin (Seki et al., 2008).

The PCC table shows that different CPR classes had different correlation values with different CYPs (Table 1). The probeset ID used for PCC calculation is presented in Supplementary Table 8. CYP716A12 had a higher positive correlation value with MtCPR1, while CYP72A63 had a higher positive value with MtCPR2. However, from the data-set, no significant difference was observed in the correlation of MtCPRs with CYP93E2 (Table 1A). From the L. japonicus database MG20 v3.0, LjCPR1 showed higher correlation values with CYP716A51, CYP72A61, and CYP93E1, as compared to LjCPR2 (Table 1B). Based on the transcriptome data of L. japonicus in the Gifu version 2.0, LjCPR1 also showed higher positive correlation value with CYP716A51 and CYP93E1, as compared to LjCPR2. However, CYP72A61 did not have strong correlation with any of the LjCPRs, but showed very high correlation value (>0.5) with CYP93E1 (Supplementary Table 9). GuCPR2 showed higher correlation with CYP716A179, CYP72A154, and CYP93E3, as compared to GuCPR1. In contrast, CYP88D6 correlated negatively with GuCPR1 (Table 1C).

To analyze the effect of different CPR classes on heterologous triterpenoid production, different CYPs and CPR classes from M. truncatula and L. japonicus were co-expressed in transgenic S. cerevisiae INVSc1, a commonly used strain for recombinant protein expression. By comparing the relative value of each metabolite using a semi-quantitative method with an internal standard, we compared the accumulation of oxidized triterpenoids resulting from the heterologous expression of different CPR-CYP pairs in transgenic yeast (Figure 6). The amount of triterpenoid backbone in the control samples was considered to be 100%. In yeast strains expressing MtCYP716A12 and LjCYP716A51, the intermediate compounds oleanolic and ursolic aldehyde could be detected based on the mass spectrum profile (Supplementary Figure 3) compared to previously reported mass spectra of oleanolic and ursolic aldehyde authentic standards (Misra et al., 2017). However, betulinic aldehyde could not be detected in the yeast extract and was co-eluted with betulinic acid in the yeast extract; (Suzuki et al., 2018) therefore, the amount of betulinic acid shown in Figure 7 corresponds to the sum of both compounds. C-30 and C-24 aldehydes were excluded because of the instability of the compound, lack of authentic standards, and no previous reports on their mass spectrum profile.

The results show that upon using INVSc1 β-amyrin-producing yeast, MtCYP716A12 and LjCYP716A51 showed higher relative amount of oleanolic acid when paired with MtCPR1, as compared to when they were paired with MtCPR2, with a comparable amount of erythrodiol and oleanolic aldehyde (Figures 6A,B). Pairing MtCYP716A12 with LjCPR1 showed higher conversion of β-amyrin into erythrodiol compared to pairing it with LjCPR2. However, pairing MtCYP716A12 with LjCPR2 showed a higher conversion rate of erythrodiol into oleanolic aldehyde and oleanolic acid than with LjCPR1 (Figure 6A). Conversely, pairing LjCYP716A51 with LjCPR2 resulted in a higher conversion rate of β-amyrin into erythrodiol than with LjCPR1, while the conversion rate of erythrodiol into oleanolic acid was higher when pairing LjCYP716A51 with LjCPR1 than LjCPR2 (Figure 6B). These trends were also similar in the yeast feeding assay results of INVSc1 strain supplemented with erythrodiol as a substrate, although the difference was not as significant as in the in vivo assay (Supplementary Figure 4). In contrast, pairing MtCYP72A63 with CPR-II resulted in significantly higher conversion of β-amyrin into 30-OH β-amyrin and 30-COOH β-amyrin levels (Figure 6C) than pairing it with CPR-I. Similar to MtCYP72A63, pairing LjCYP72A61 with CPR-IIs resulted in significantly higher conversion of β-amyrin into sophoradiol (Figure 6D). In case of the legume-specific CYP subfamily, CYP93E1 paired with CPR-II showed higher conversion of β-amyrin into 24-COOH β-amyrin, as compared to CYP93E1 paired with CPR-I (Figure 6E).

CYP716As can also catalyze oxidations on different triterpenoid backbone. Therefore, to confirm that the effect of the CPR on the conversion rate can also be applied to different triterpenoid backbone, we then co-expressed MtCYP716A12 in α-amyrin and lupeol-producing strains (Figure 7). The result showed that MtCYP716A12 paired with MtCPR1 and LjCPR1 showed higher conversion of α-amyrin into ursolic acid and uvaol, respectively (Figure 7A). A similar trend was observed in the lupeol-producing strain in which MtCYP716A12 paired with MtCPR1 showed a higher conversion rate of betulin into betulinic acid compared to that with MtCPR2, and MtCYP716A12 paired with LjCPR1 showed a higher conversion rate of lupeol to betulin and betulinic acid compared to LjCPR2 (Figure 7B). All the chromatograms of yeast in vivo are shown in Supplementary Figures 5–14.

This yeast in vivo assay might provide insight into CPR-CYP preferences, which are better justified by their protein expression levels. Quantifying the CPR protein itself would not be sufficient because triterpenoids are the products of CYP. However, protein quantification of CYP remains challenging. We propose this study to analyze the effect of CPR class on heterologous production in yeast, but not on the catalytic activity of CPR or CYP. We admit that this is a limitation of yeast in vivo studies. Therefore, further in-depth in vitro analysis of the effect of CPR class on CYP activity might be needed in future experiments to provide more justifiable results.

Based on the multiple sequence alignment result, the main difference of CPR class I and II is located in the N-terminal domain. CPR class I has shorter amino acid sequence than CPR class II (Figure 3). To analyze whether this N-terminal domain is responsible for the differences in the conversion rate of β-amyrin when CPR was co-expressed with different CYPs in transgenic yeast (Figure 6), we investigated the effect of this domain by swapping the N-terminal membrane sequence of MtCPR1 and MtCPR2. MtCPRs were used as the representative since M. truncatula is a model legume. We speculated if the N-terminal domain of CPR is indeed responsible for the CPR activity, swapping the N-terminal domain of CPR class I and II will affect the triterpenoids profile.

The location of transmembrane helix of MtCPR1 and MtCPR2 was predicted in silico (Supplementary Figure 1). Amino acids of MtCPR1 number 1–26 are located inside the ER, number 27–46 are the transmembrane helix, and number 47–692 are located outside the ER (Supplementary Figure 1A). Amino acids of MtCPR2 number 1–40 are located inside the ER, number 41–60 are the transmembrane helix, and number 61–701 are located outside the ER (Supplementary Figure 1B). In this study, the first 46 amino acids of MtCPR1 located inside the ER together with the transmembrane helix were swapped with the first 60 amino acids of MtCPR2. The resulting chimeric MtCPRs were N-terminal membrane sequence of MtCPR1 fused with truncated-N-terminal MtCPR2 (M1N-M2C) and N-terminal membrane sequence of MtCPR2 fused with truncated-N-terminal MtCPR1 (M2N-M1C) (Figure 8A).

Yeast in vivo assay was then performed by co-expressing the chimeric MtCPRs with CYP716As and CYP72As to observe its effect on triterpenoids conversion rate (Figures 8B–E). The results showed that by swapping the N-terminal domain of CPR class I with that of class II, the product from CYP716A12 and CYP716A51 were significantly reduced (Figures 8B,C). The production of erythrodiol by CYP716A51 co-expressed with mutant CPRs was even lower than the background level observed in native yeast CPR (Figure 8C). A similar effect was also seen in CYP72A61, where the chimeric CPR caused a significant decrease in conversion of β-amyrin into sophoradiol (Figure 8D). Interestingly, only CYP72A63 was not affected by the chimeric MtCPRs. Moreover, the results showed that CYP72A63 paired with the N-terminal membrane sequence of MtCPR2 fused with truncated N-terminal MtCPR1 (M2N-M1C) showed higher conversion of β-amyrin into 30-OH β-amyrin and 30-COOH β-amyrin (Figure 8E).

To further analyze the effect of different CPR classes on improving heterologous production, different CPR classes were co-expressed with different CYP families in the PSIII strain. PSIII is an engineered yeast strain that has been optimized for triterpenoid biosynthesis by means of some modifications in the upstream mevalonate pathway. This strain is modified from the PSII strain constructed in previous research, which can produce a 100 × times higher triterpene backbone as a substrate than its parental strain (Srisawat et al., 2020). In contrast to the INVSc1 results, in the PSIII strain, almost all the triterpenoid sapogenin production profiles had a lower ratio of β-amyrin, and the trends of CYP-CPR preferences observed in the INVSc1 strain seem to be less obvious in this strain (Figure 9). Nevertheless, some notable changes were observed in the PSIII strain. MtCYP716A12 paired with LjCPR1 that previously showed different conversion ratios of triterpenoids with LjCPR2 in the INVSc1 strain, showed similar conversion ratio profiles in the PSIII strain (Figure 9A). The ratio of oleanolic acid to erythrodiol was lower in PSIII strain co-expressing MtCYP716A12 (Figure 9A). The PSIII strain co-expressing MtCYP72A63 resulted in much lower 30-COOH β-amyrin to 30-OH β-amyrin ratio, as compared to the INVSc1 strain (Figure 9C). Co-expressing all CPRs with LjCYP72A61 and MtCYP716A51 resulted in similar conversion ratio in PSIII strain (Figure 9D). However, similar to the INVSc1 strain, co-expressing MtCYP716A12 with MtCPR1 in the PSIII strain showed higher oleanolic aldehyde and oleanolic acid conversion rates than MtCPR2. In addition, co-expressing MtCYP72A63 with LjCPR2 in the PSIII strain showed higher 30-OH β-amyrin and 30-COOH β-amyrin conversion rates compared to that with LjCPR1.

Co-expression of different CYP and CPR classes from model legumes, Medicago truncatula and Lotus japonicus, have suggested that CYP shows a preference for certain CPR classes. Therefore, to further study the effect of CYP-CPR pairs, we added CPRs from a non-model legume species, Glycyrrhiza uralensis, which produces glycyrrhizin, a highly valuable metabolite found only in Glycyrrhiza sp. Furthermore, unlike CYP716As and CYP72As, which are commonly found in plant species, GuCYP88D6 uniquely presents in Glycyrrhiza sp. and is very specific in function. Therefore, it is interesting to facilitate contrasting results. GuCYP88D6 catalyzes the two-step oxidation at the C-11 position of β-amyrin to produce 11-oxo β-amyrin, the precursor for glycyrrhizin (Seki et al., 2008). The PCC value of GuCYP88D6 with GuCPRs also indicates that GuCYP88D6 might have different preferences for GuCPR classes I and II, which shows a negative correlation with GuCPR1 (Table 1C). Corresponding to this hypothesis, the results showed that co-expressing GuCYP88D6 with class II GuCPR in the PSIII strain showed the highest conversion of β-amyrin into 11-oxo β-amyrin, as compared to the other CPRs tested in this study (Figure 10).