Deep mutational scanning of CYP2C19 reveals a substrate specificity-abundance tradeoff

Cytochrome P450s (CYPs) are a family of enzymes responsible for metabolizing nearly 80% of small molecule drugs. Variants in CYPs can substantially alter drug metabolism, which may result in improper dosing and severe adverse drug reactions. CYPs have low sequence conservation, making it difficult to anticipate whether variant effects measured in one CYP may extend to others based on sequence alone. Even closely related CYPs, like CYP2C9 and its closest homolog CYP2C19, have distinct phenotypic properties despite sharing 92% amino acid sequence identity. Thus, we used Variant Abundance by Massively Parallel sequencing (VAMP-seq) to measure the steady-state protein abundance, a proxy for protein stability, of 7,660 missense variants in CYP2C19 expressed in cultured human cells. Our results confirmed positions and structural features critical for CYP function and revealed how variants at positions conserved across all eukaryotic CYPs influence abundance. We jointly analyzed 4,670 variants whose abundance was measured in both CYP2C19 and CYP2C9, finding that the homologs have different variant abundances in substrate recognition sites within the hydrophobic core, and that substitutions in some regions reduced abundance in CYP2C19 but not CYP2C9. We also measured the abundance of all single and some multiple WT amino acid exchanges between CYP2C19 and CYP2C9. While most exchanges had no effect, substitutions in substrate recognition site 4 (SRS4) reduced abundance in CYP2C19. When nearby amino acids were exchanged in double and triple mutants, we found distinct interactions between the sites in CYP2C19 and CYP2C9, revealing a region that is partially responsible for the difference in thermodynamic stability between the two homologs. Since these positions are also important for determining substrate specificity, there may be an evolutionary tradeoff between stability and altered enzymatic function. Finally, we used our data to analyze 368 previously unannotated human variants, finding that 43% had decreased abundance. Thus, by comparing variant effects between two closely related and important human genes, we have uncovered regions underlying their functional differences and paved the way for a more complete understanding of one of the most versatile families of enzymes.


Introduction
Nearly 20,000 cytochrome P450 (CYP) heme monooxygenases have been identified across all domains of life 1 . CYPs catalyze a wide range of reactions with a diverse set of substrates, making them some of the most versatile enzymes in existence 2,3 . The 57 human CYP genes are grouped into 18 families with 43 subfamilies 4 , highlighting their genetic heterogeneity even within a single species. Despite their genetic and functional diversity, key structural and topological features of CYPs are highly conserved 5,6 . However, the relationship between CYP genetic variation, structure, and function is far from fully elucidated. For example, within CYP family 2, subfamily C (CYP2C), CYP2C19 (MIM: 124020, 609535) and CYP2C9 (MIM: 601130) are the most closely related subfamily members, sharing 92% amino acid sequence identity. Their protein structures have nearly identical organization, with the largest deviations between their Cα backbones in the substrate binding cavity being only~3 Å 7 . Yet, the two homologs are functionally distinct, with largely disparate sets of substrates 8,9 and divergent membrane interactions 10 . Moreover, CYP2C19's melting temperature is~11℃ higher than CYP2C9's 11 . Thus, even between these close homolog CYPs, the 43 diverged positions drive large functional differences.
Understanding the functional impact of variants across CYPs is particularly important because~12 of the 57 human CYPs contribute to metabolizing 70-80% of currently prescribed drugs that are processed by enzymes for elimination. Of those 70-80% CYP-metabolized drugs, CYP2C19 and CYP2C9 account for 20-30% 12 . Genetic variation in CYPs can substantially alter individual drug response leading to adverse drug reactions (ADRs), which are among the leading causes of morbidity and mortality 13,14 , and cost an estimated $30.1 billion annually 15 . To provide clinicians guidance for treating individuals with CYP variants, the Clinical Pharmacogenetics Implementation Consortium (CPIC) categorizes CYP genes into star (*) allele haplotypes according to enzymatic function: normal function, decreased function, no function, and increased function 16,17 . Genetic testing and employment of CPIC guidelines can prevent many ADRs. For example, up to 30% of the population may have a CYP2C19 variant with reduced function 18 which may result in impaired activation of the antiplatelet drug clopidogrel. Genotyping for CYP2C19 loss of function variants can avoid major adverse cardiovascular events [19][20][21] . However, only a very small number of CYP variants have established functional consequences, and it is unknown to what degree variant effects in one CYP can be applied to others.
Measuring CYP variant function individually is laborious and low throughput, but massively multiplexed methods can be used instead. Variant Abundance by Massively Parallel sequencing (VAMP-seq) measures steady-state protein abundance of thousands of variants in parallel 22 . In VAMP-seq, as in other similar methods [23][24][25][26] , steady-state protein abundance is used as a proxy for protein stability. Steady-state protein abundance refers to the final concentration of a protein when its rates of synthesis and decay are balanced 27 . To control for variation in protein synthesis, the VAMP-seq vector uses a single promoter to express two fluorescent reporters from a single mRNA transcript using an internal ribosome entry site (IRES). In each cell, the fluorescent signal of the target gene fused to enhanced green fluorescent protein (eGFP) is normalized to mCherry fluorescence meaning changes in fluorescence are due to degradation. The resulting measurements correlate with changes in thermodynamic stability, indicating that reduced steady-state abundance is likely due to loss of protein stability 22,24,28,29 .
Previously, we used the VAMP-seq assay 22 to measure the abundance of 6,370 of 9,780 possible missense variants in CYP2C9 30 . From the resulting variant effect map, we identified patterns of loss of abundance that revealed mutationally sensitive regions of the protein.
Additionally, we revealed hundreds of variants with reduced abundance in the human population in addition to providing variant effect measurements for thousands of variants not yet observed 30 .
Here, we used VAMP-seq to measure 7,660 of 9,780 possible missense variants of CYP2C19. We identified 4,698 variants that likely result in reduced protein abundance, with 1,122 of those exhibiting complete abundance loss equivalent to nonsense mutations. We first analyzed positions conserved across all eukaryotic CYPs, revealing that all but six of the 58 conserved positions were intolerant of substitutions. Four of the tolerant positions were catalytically important sites buried in the hydrophobic core where mutations are nearly always deleterious, suggesting that some sites critical for enzyme function may not impact abundance.
We jointly analyzed the CYP2C19 and CYP2C9 variant abundance dataset, and found 2,366 variants whose abundance differed between the two enzymes. Most differences were of small effect, though a fraction of the differences were large. While nearly all sites had at least one variant that differed, 83 of 489 (17%) sites were significantly different between the homologs. CYP2C9 had higher mutational tolerance in its hydrophobic core than CYP2C19, and variants in the structurally conserved K' helix are highly deleterious in CYP2C19, but tolerated in CYP2C9, even though all K' positions contain the same amino acid in both homologs. We analyzed WT amino acid exchanges between CYP2C19 and CYP2C9, revealing that sequence differences in a set of diverged positions partially drive the differences in thermodynamic stability between the homologs. These divergent positions are also important for substrate specificity, suggesting that reduced thermodynamic stability in CYP2C9 may have been evolutionarily tolerated in exchange for functional benefit 31 . Finally, we analyzed the effects of human CYP2C19 variants.
Our abundance scores are largely concordant with existing functional annotations indicating that, like for many other proteins, loss of abundance accounts for the majority of loss-of-function alleles. We provided abundance scores for 368 out of 408 (90.2%) previously unannotated missense variants in the gnomAD database. Thus, by conducting the first comparative analysis of closely related CYPs using large-scale variant effect data we provide fundamental insights into common CYP structural features that differentially impact abundance between CYP2C19 and CYP2C9. We also provide functional annotations for human CYP2C19 variants which could be used to improve genotype-guided dosing of drugs metabolized by CYP2C19.

Multiplexed measurement of CYP2C19 variant abundance
We used VAMP-seq to simultaneously measure the steady-state abundance of CYP2C19 variants in cultured human cells 22,30 (Fig. 1A). VAMP-seq relies on two fluorescent reporters: GFP fused to each CYP2C19 variant to read out abundance, and mCherry expressed via an internal ribosome entry site (IRES) as a transcriptional control. Because CYP2C19 is N-terminally inserted into the endoplasmic reticulum membrane, we fused GFP onto the C-terminus, as we did for a previous VAMP-seq experiment on CYP2C9 30 . Expression of the wild type (WT) CYP2C19 C-terminal GFP fusion led to strong fluorescent signal, and R433W, a known destabilizing CYP2C19 variant, had substantially lower signal indicating that the C-terminal GFP fusion construct was compatible with VAMP-seq (Fig. 1B).
We introduced a barcoded library of CYP2C19 variants into HEK293T cells using a recombinase-based landing pad, such that each cell expressed only one variant 32,33 . Cells were sorted into quartile bins based on the ratio of GFP:mCherry fluorescence. Each bin was deeply sequenced, variant-associated barcodes were counted, and abundance scores were calculated based on weighted average of barcode frequencies across bins (Fig. 1A). Abundance scores were highly correlated between seven replicate sorting experiments arising from three independent library recombinations (Supp. Fig. 1A-D; Pearson's R = 0.82 -0.98). Replicate scores were averaged, filtered (Supp . Fig 2A-D) and normalized such that the median nonsense variant had a score of 0 and WT had a score of 1 22,30 . Our final data set contained abundance scores for 8,480 of 10,290 (82%) possible variants, of which 7,660 were missense, 316 were nonsense, and 504 were synonymous.
Abundance scores of synonymous and nonsense variants were well separated, with the missense variant distribution spread between nonsense and synonymous variants (Fig. 1C).
Individually measured GFP:mCherry ratios for 10 variants spanning the range of abundance scores were highly correlated with VAMP-seq scores ( Fig. 1D; Pearson's R= 0.92). Our results are also highly consistent with a smaller scale VAMP-seq experiment encompassing 121 variants 34 (Supp. Fig. 3A, Pearson's R = 0.74). Thus, our VAMP-seq derived abundance scores faithfully reproduced variant abundance. Lastly, we classified variants according to their abundance score relative to the range of scores from nonsense and synonymous variants (Supp. Fig. 3B, Fig. 1C,E). The majority (58%, 4,620 variants) of missense variants decreased abundance (Fig. 1E).

Mutational tolerance at conserved CYP2C19 positions reflects function
We visualized the abundance scores as a variant effect map ( Fig. 2A) and projected position-averaged scores onto the CYP2C19 structure (Fig. 2B). Many of the low abundance variants occur within α-helices and β-sheets ( Fig. 2A-B), especially in amino acids on interior α-helix turns and in regions closer to the protein core (Fig. 2B).
While CYPs vary widely in sequence, key structural and functional features are highly In clusters 1, 2, and 3 nearly all substitutions, except those of the same biophysical type, reduced abundance (Fig. 2C, Supp. Fig. 4A). In cluster 4, substitutions caused moderate loss of abundance, with no consistent pattern across all positions. The sole exceptions were two of the three positions where glycine was the WT amino acid. These positions tolerated alanine and cysteine substitutions suggesting that amino acid size is an important factor. Cluster 5 contained M136, A297, E300, T301, K322 and I362, all of which were substantially more tolerant of mutations than the other conserved sites indicating that they are critical to CYP2C19 function but not abundance. The combined conservation and tolerance of M136 and K322 can be explained by the fact that these positions are located on the surface of the protein and that they likely bind to the critical cofactor cytochrome P450 reductase (CPR), as they do in the closely related CYP2C9 39,40 . However, amino acids A297, E300, T301, and I362 are buried in the hydrophobic core making their mutational tolerance more challenging to explain (Supp. Fig 4B).
Positions 297 and 362 influence substrate specificity, and >80% eukaryotic CYPs have hydrophobic amino acids at these positions 38 . Surprisingly, while substitutions are tolerated at these positions, some hydrophobic substitutions elicit moderate reductions in abundance. T301 is a critical threonine for oxygen activation and catalysis in 2C19 7,41-43 and contains hydrogen-bonding amino acids in >80% eukaryotic CYPs at this position, and most substitutions did not appreciably reduce abundance. Finally, E300 stabilizes a water network during proton delivery 43 , and substitutions other than aspartic acid were tolerated at this position.
Thus, substitutions at nearly all conserved positions caused reduced abundance.
However, positions 136, 297, 300, 301, 322, and 362, which participate in catalysis or cofactor binding, were largely tolerant of substitutions despite their location in the hydrophobic core of CYP2C19. We speculate that this tolerance is a consequence of the dynamic and flexible nature of CYP active sites, making these positions important for catalytic activity but not folding and stability 44 .
Comparing variant abundance effects between CYP2C19 and CYP2C9 reveals core-stabilizing regions with distinct mutational tolerance Next, we investigated variant effect patterns in CYP2C19 compared to its closest homolog, CYP2C9. CYP2C19 and CYP2C9 share 92% protein sequence identity and nearly identical crystal structures (Supp. Fig. 5, RMSD = 0.596 Å). However, they have important functional differences, notably their substrate profiles and membrane interactions 8,10,45,46 .
Moreover, the temperature at which they lose the ability to bind their heme cofactor, which reflects thermodynamic stability 47 , differs by 11°C 11 . Thus, small differences in sequence and structure translate into distinct functional and phenotypic characteristics.
To understand how these functional differences arise, we sought to estimate how much the abundance score is shifted between the CYP2C19 abundance data presented here and abundance data from a previous VAMP-seq experiment we conducted on CYP2C9 variants 30 .
The combined dataset contained 4,670 variants whose abundance was scored in both CYPs.
Most variants had similar effects in both homologs, though a subset of variants had large shifts in effects (Fig. 3A, Pearson's r = 0.77). While some shifts are likely due to actual biological differences, others are due to the noise inherent in any high-throughput experiment. To identify which shifts are most likely due to signal, we inferred shifts using a joint-modeling approach called multidms 48 . This open source package models two or more independent DMS experiments within the same predictive space, and utilizes appropriate forms of parameter regularization with the goal of identifying shifts in mutational effect We inferred shifts as the difference in abundance score in CYP2C19 relative to CYP2C9, such that positive shifts indicate a higher abundance score in CYP2C19 and vice versa. We tested regularization weights ranging from 0.0 to 1e-4, and selected 1e-5 as the optimal value for subsequent analysis (Supp. We calculated the mean of the shift values at each position to reveal the effect of regional and structural features (Fig. 3C). We identified positions with mean shift values that differed significantly from 0 using a randomization test (Supp. Fig. 9B). The region with the largest mean shift values was in the K' helix, which is part of a region that is both highly mobile and critical for packing of the hydrophobic core 5,49 (Fig. 3C). In this region, mean shift values were negative, meaning that substitutions were more deleterious in CYP2C19 than in CYP2C9 ( Fig. 3D, Supp. Fig. 10A).
Overall, CYP2C19 was more mutationally tolerant than CYP2C9 in the D, E, I, L, J, and J' helices ( Fig. 3C, Supp. Fig. 9A, Supp. Fig. 10A, Bi-ii), which form the majority of the hydrophobic core 5,49 . The sites that were most differentially tolerant in these helices were on portions of the helices that sit outside of the hydrophobic core (Supp. Fig. 10Bi-ii). Conversely, CYP2C9 was more mutationally tolerant than CYP2C19 at positions within the hydrophobic core near important sites for heme positioning and function (Supp. Fig. 10A). Many of these heme-associated positions reside within substrate recognition sites (SRSs) 50 , and CYP2C9 was more mutationally tolerant than CYP2C19 in SRSs relative to the other regions of the protein ( Fig. 3C, E-F). The mutational tolerances of positions in SRSs that were not heme-associated were similar between the homologs (Fig. 3F).
We also examined whether differences between CYP2C19 and CYP2C9 could be explained by sensitivity to variants of different biophysical types or by differences in the structures of the two homologs. However, we found that neither homolog is more sensitive to particular types of substitutions (Supp. Fig. 11A), and that shift values were unrelated to the distance between positions in the CYP2C19 and CYP2C9 crystal structures (Supp. Fig. 11B) Thus, comparison of variant effects between CYP2C19 and its closest homolog CYP2C9 revealed that CYP2C19's K' helix and, to a lesser extent, heme-associated positions in the hydrophobic core were more sensitive to mutation than CYP2C9, but that CYP2C19 was less sensitive than CYP2C9 to substitutions in other regions flanking the hydrophobic core.

Amino acid swaps reveal homolog-specific constraints on abundance at sites influencing substrate specificity
We investigated abundance shifts at all variants comparing CYP2C19 and CYP2C9.
However, the phenotypic differences in substrate recognition, membrane interaction, and thermodynamic stability between the two homologs must be driven by divergent sites. While most divergent sites are not localized to the catalytic site, some are critical for substrate specificity, regiospecificity, and stereospecificity [51][52][53][54][55][56] . In many cases, evolutionary pressures result in a protein's reduction of thermodynamic stability in exchange for new functionality 31 . We wondered whether we could link the differences in thermodynamic stability between the homologs to substrate specificity by measuring the protein abundance of the 43 divergent sites.
Thus, we investigated the abundance of the variants that partially convert CYP2C19 to CYP2C9 and vice versa.
We had abundance scores for 33 of the 43 2C19 variants that install the WT CYP2C9 amino acid, (e.g. CYP2C19→CYP2C9). For an exhaustive analysis, we individually measured GFP:mCherry fluorescence for each of the 10 CYP2C19→CYP2C9 variants not present in our abundance data (Fig. 4A). All but three CYP2C19→CYP2C9 substitutions were well tolerated.
R261Q and L295F caused modest loss of abundance. Position 295 is critical for the specificity of CYP2C19 for S-mephenytoin and for the specificity of CYP2C9 for diclonefac 45,57 whereas R261Q has not been studied in either homolog. V288E caused the largest loss of abundance of all of the swaps and was classified as "nonsense-like" (Fig. 4A). Position 288 alone does not have a known functional role. However, positions 241, 288, and 289 together have been suggested to play a role in abundance and substrate specificity 45,51,52,54,57 . In the CYP2C9 structure, K241 interacts electrostatically with E288 and hydrogen bonds with N289 to stabilize a region of the SRS4 in the I helix 51,56 . In CYP2C19, all three positions have different amino acids, E241, V288, and I289, and thus no electrostatic interaction between E241 and V288.
Thus, the loss of abundance caused by V288E in CYP2C19 was likely due to the introduction of an electrostatic clash between E241 and V288E (Fig. 4B). Since all three positions contribute to functional differences between CYP2C19 and CYP2C9, we sought to understand how positions 241, 288, and 289 might interact to influence abundance.
First, we individually measured the abundance of single and double mutants at positions When individually measured, E241K had no effect on CYP2C19 abundance and V288E profoundly reduced abundance, the same effects we measured using VAMP-seq (Fig. 4C).
Combining E241K and V288E partially restored CYP2C19 abundance. CYP2C9 K241E only modestly reduced abundance, even though this variant putatively results in an electrostatic clash similar to the one that dramatically reduced CYP2C19 abundance. CYP2C9 E288V had no effect on abundance, suggesting that the native K241-E288 electrostatic interaction likely does not contribute appreciably to thermodynamic stability 51 . Combining K241E and E288V fully restored CYP2C9 abundance (Fig. 4C). Thus, both homologs have a similar pattern, with installation of a second negative charge disrupting abundance. Elimination of one of the two negative charges even with variants from the other homolog restored abundance, although to differing degrees in each homolog.
Next, to incorporate position 289 into our analysis, we measured the abundance of the CYP2C19→CYP2C9 and CYP2C9→CYP2C19 241, 288, 289 triple mutants (Fig. 4C). The CYP2C19→CYP2C9 triple mutant had a modestly reduced abundance relative to CYP2C19 WT, largely restoring the low abundance of the E241K, V288E double mutant. The CYP2C9→CYP2C19 triple mutant had an abundance equivalent to CYP2C9 WT and to each of the two double mutants. Thus, while the interaction between these three positions is complex, it seems likely that sequence changes in this region of the protein contribute to the increased thermodynamic stability of CYP2C19.
CYP2C19 variants can increase, decrease, or eliminate an individual's ability to metabolize many important drugs, and knowing variant function can help avoid severe and expensive adverse events [13][14][15]58,59 . For example, the anti-platelet drug clopidogrel is activated by CYP2C19. Thus, individuals with deleterious CYP2C19 variants experience reduced or non-existent benefit from clopidogrel, requiring higher doses or alternative drugs. Genetic testing for CYP2C19 variants prior to clopidogrel treatment is important for avoiding major adverse cardiovascular events 19,20 . PharmVar is a repository for pharmacogene allelic variation and functional information, including CYP2C19. Alleles in PharmVar are known as "star alleles," and annotated using star notation 17 . For example CYP2C19*5 refers to R433W. Despite decades of study, 10 of the 39 CYP2C19 star alleles are of uncertain function. Thus, we analyzed the functional effects of CYP2C19 alleles in PharmVar. All four PharmVar "normal function" alleles had WT-like abundance scores (Fig. 5A). Of the eight "decreased" and "no function" alleles, six were low abundance. The remaining two decreased/no function alleles, *6 and *9, were WT-like abundance. PharmVar lists the *6 allele (R132Q) as "no function" with "definitive" evidence; however we measured a WT-like abundance score of 1.06 (95% CI 1.09 -1.03). This strongly suggests that the *6 allele's loss of function results from disrupted enzymatic activity rather than loss of abundance. Consistent with this interpretation, the *6 allele is intact enough to bind its heme cofactor, but has a decreased ability to metabolize substrates and disrupted electron flow from cytochrome P450 reductase 60,61 . Allele *9 (R144H) had an abundance score of 0.914 (95% CI 0.956-0.872). Consistent with these results, *9 has WT-like affinity for cytochrome P450 reductase 62 , suggesting that it is at least partially folded. However, our abundance results conflict with another, smaller scale VAMP-seq experiment in which CYP2C19*9 was identified as "decreased" abundance 34 with less than~50% of WT abundance.
We therefore individually validated this variant and reaffirmed its WT-like abundance in our hands (Supp. Fig. 12). In light of *9's moderately reduced activity against mephenytoin, ability to bind cytochrome P450 reductase normally and WT-like abundance in our assay, we suggest that, like *6, *9 has normal abundance but decreased catalytic activity. Overall, six of the eight known loss of function alleles had reduced abundance, and all normal function alleles had WT-like abundance. While variants with WT-like abundance could have low or no function, reflecting the many ways function can be compromised, low abundance variants were always low or no function alleles. Thus, abundance is a powerful method for identifying loss of function variants. Of the ten PharmVar alleles with uncertain function, four were present in our VAMP-seq library, and we found that *30 and *23 had decreased abundance strongly suggesting that they would disrupt drug metabolism (Fig. 5A).
As sequencing and genetic testing are more widely deployed, rare variants with unknown clinical consequences are being identified at an exponentially increasing rate 63 .
Reflecting this reality, the annotated alleles in the PharmVar database are only a fraction of all of the CYP2C19 variants discovered so far. There are 408 unique CYP2C19 missense variants in the exome database gnomAD, 390 of which have no CPIC annotation or functional information.
We annotated 368 (90.2%) of the variants in gnomAD (Fig 5B) and identified 131 (35.6%) variants with "decreased" abundance and 29 (7.88%) with "nonsense-like" or "possibly nonsense-like" abundance relative to WT, strongly suggesting that these variants have decreased or no function. We annotated 210 (57.0%) variants as "WT-like" or "possibly WT-like," indicating that these variants may have normal function. However, assessment of enzymatic activity would be needed to definitively determine if these "WT-like" or "possibly WT-like" variants have normal function since variants can eliminate activity without affecting protein stability. These results are broadly consistent with a study that genotyped 2.29 million participants for CYP2C19*2, *3, and *17 alleles. The study discovered that *2 was present in 15.2%, *3 in 0.3%, and *17 in 20.4% of individuals, and nearly 60% had at least one of these star alleles 64 . Thus, CYP2C19 variants with reduced abundance appear common in the population.

Discussion
The CYP family tree spans all animal kingdoms and comprises an exceptionally versatile set of enzymes. Understanding the phenotypic consequences of natural variation in human CYPs is particularly important since they catalyze the metabolism of most drugs currently in use.
However, even closely related CYPs, like CYP2C19 and CYP2C9, are functionally distinct, and the underlying causes of these distinctions are largely unknown. We used VAMP-seq to measure the abundance of 7,660 CYP2C19 missense variants. In addition to confirming positions known to be critical for CYPs structure and function, we revealed that variants at four conserved positions in the hydrophobic core do not impact CYP2C19 abundance. By jointly analyzing 4,670 shared CYP2C19 and CYP2C9 abundance scores, we discovered regions where the two homologs have different mutational tolerances. CYP2C9 has a more tolerant hydrophobic core, whereas CYP2C19 is more tolerant in regions surrounding the core. We measured the abundance of WT amino acid swaps between CYP2C19 and CYP2C9, discovering a region likely responsible for at least some of the thermodynamic stability difference between the homologs. Finally, our abundance scores identify known reduced activity CYP2C19 variants with high fidelity, and indicates that two star alleles of unknown function, *30 and *23, likely have reduced abundance. We also evaluated 368 of the 408 human CYP2C19 variants with no prior annotation. Notably, 43% of these variants are low abundance, warranting follow-up studies to measure their impacts drug metabolism. impact on substrate specificity but also in their specialized role, as they primarily influence specific functions rather than overall protein abundance.
We also jointly analyzed CYP2C19 and CYP2C9 30 abundance scans using multidms 48 to find variants with different abundances. Variants in the K'-helix reduced abundance in CYP2C19 but were tolerated in CYP2C9, suggesting a markedly different mutational tolerance of K' in the enzymes despite having identical WT amino acids. Moreover, CYP2C9 had a more tolerant hydrophobic core than CYP2C19, especially in SRSs that contain heme-associated positions.
CYP2C9's higher mutational tolerance in its core may indicate more flexibility. We speculate that, since the flexibility of CYP active sites is correlated with its promiscuity 44,67 , this may allow CYP2C9 to bind more substrates 9 .
Variants capable of imparting novel function, like those that alter substrate specificity, often reduce thermodynamic stability 31  Likewise, CYP2C9 K241E introduces the same opposing negative charge adjacent to E288 and also reduces abundance. This pattern of loss of abundance suggests that CYP2C9 may have evolved from CYP2C19. This is because CYP2C9 is the only enzyme in the family that has a negatively charged amino acid at 288, meaning that the ancestral sequence had valine at position 288 56 . Thus, the ancestral CYP2C9 likely acquired E241K or I289N first, both of which partially ameliorate the loss of abundance induced by V288E. We did not find that any combination of swaps at these three positions could fully restore CYP2C19 V288E abundance.
One possibility is that swaps at other sites, which by themselves do not affect CYP2C19 abundance, could fully rescue V288E. However, the reduced abundance of the CYP2C19 E241K-V288E-I289N variant is in line with CYP2C9's reduced thermodynamic stability.
CYP2C9's new substrate binding capabilities apparently made this loss of thermodynamic stability evolutionarily tolerable. Importantly, our VAMP-seq derived abundance data has limitations. First, the abundance of each variant arises from the balance between protein synthesis and degradation driven by cellular protein quality control systems. In VAMP-seq, all variants are expressed transgenically from an inducible promoter, and an mCherry reporter expressed via an IRES is used to control from cell-to-cell variation in expression. Thus, VAMP-seq derived abundance scores largely reflect changes in degradation which are, in turn, generally driven by stability-related changes in protein folding 22,24,28,29 . However, abundance changes can result from other mechanisms such as changes to degron sequences or protein localization. Second, we expressed the CYP2C19 cDNA from an inducible promoter, meaning we cannot detect variants that induce splicing defects or affect transcriptional regulation. Third, variants can affect function without affecting abundance. For example, variants may disrupt a critical substrate binding position or prohibit binding to critical cofactors like cytochrome P450 reductase (CPR) or cytochrome b5. Therefore, while variants that we identified with reduced abundance are likely to alter drug metabolism, variants with WT-like abundance may not necessarily have normal function. Finally, the VAMP-seq assay depends on fluorescent reporters and fluorescence activated cell sorting. As a result, subtle changes in abundance are difficult to discern.
In the future, we envision intersecting mutational scans from important CYPs in other subfamilies such as CYP2D6 and CYP3A4. Since multidms is capable of jointly analyzing more than two scans, it will be a powerful tool to compare mutational effects across additional CYPs, helping us understand the extent to which variant effects are conserved across the family. In addition to improving personalized drug dosing, such comprehensive profiling could expand our overall understanding of CYP function.

General reagents
Unless otherwise noted, all chemicals were obtained from (MilliporeSigma) and all the enzymes were obtained from New England Biolabs. All cell culture reagents were purchased from Thermo Fisher unless otherwise noted. All plasmids and oligonucleotides used in this study are listed in Table S3.

Growth media and culturing techniques
HEK293T cells (ATCC CRL-3216) and the derived landing pad cell line were cultured in Dulbecco's modified Eagle medium supplemented with 10% fetal bovine serum, 100U/mL penicillin, and 0.1 mg/mL streptomycin. Landing pad expression was induced with doxycycline at a final media concentration of 2.5 μg/mL. Cells were passaged by by detachment with trypsin 0.5%. All cell lines were tested negative for mycoplasma on a monthly basis.

Library mutagenesis
The CYP2C19 library was constructed using inverse PCR-based site-directed saturation mutagenesis 68 Table S2 for details). To generate barcoded libraries, the variant library was first digested with SacII and AflII at 37℃ for 1h, followed by heat inactivation at 65℃ for 20 minutes. We ordered barcode oligos with 18bp random sequences from IDT, resuspended them at 100 μM, and annealed them by To obtain more accurate library counts, we sequenced the libraries with Illumina sequencing. The forward and reverse reads were merged using Pear 69 , and barcode counts were estimated using Bartender 70 . Barcodes with fewer than 10 reads were filtered out, resulting in~200,000 unique barcodes for an average of 21x coverage.

PacBio sequencing for barcode-variant mapping
PacBio sequencing libraries were generated with SMRTbell Express Template Prep Kit 3.0 (Pacific Biosciences) according to manufacturer's instructions. The barcoded variant sequences were excised using restriction enzymes NheI-HF and HindIII-HF and purified with AMPure PB beads (Pacific Biosciences 100-265-900) at a 1:1 ratio of beads to DNA. Following end-repair, A-tail attachment, and ligation, the assembled product was extracted using a BluePippin instrument (Sage Science, BLU0001) using a 0.75% agarose precast cassette (Sage Science, BLF7510). Library purity and size was confirmed by 4200 TapeStation (Agilent, G2991BA) before sequencing. Samples were submitted to University of Washington PacBio Sequencing Services and sequenced on one SMRT cell in a Sequel II v2.0 run using a 15 hour movie.
We filtered long reads for a minimum of 3 passes. We then analyzed the circular consensus reads (CCSs) using PacRAT to identify and link the gene variants with the barcode region 71 . The filtered barcode-variant library contained 12,559 unique nucleotide sequences tagged by 176,372 unique barcodes (see Supplemental Table S3 for details).

FACS-based deep mutational scan (VAMP-seq)
HEK293T with a Bxb1 serine recombinase landing pad with an inducible Caspase 9 cassette (HEK293T-LLP-iCasp9) 33  Transfected HEK293T cells were sorted using a BD AriaIII sorter. Cells were gated for live, recombined singlets. In recombined cells, the ratio of GFP:mCherry fluorescence was calculated and plotted as a histogram. The histogram was split into four quartiles. Each quartile was sorted into separate 5 mL tubes. Cells from each bin were grown out for 1-2 days to ensure enough DNA for sequencing. Three biological replicates from separate transfections were collected for the FACS-based deep mutational scan. instructions. For each sample, 1 μL sample was mixed with 3 μL of Sample Buffer, thoroughly mixed, and run on a D1000 ScreenTape (Agilent Technologies) using an internal electronic ladder. The bands were quantified using the TapeStation analysis software. The samples were then pooled in equal amounts, loaded onto a 1% agarose gel with SYBR Safe, and then the gel was extracted using a freeze and squeeze column (Bio-Rad). Finally, the quantification of the pooled sample was done with the Qubit™ 1X dsDNA Assay Kit broad range (Q33266).

Library sequence analysis
Barcode sequences were trimmed and filtered for a minimum base quality of Q20 using the FASTX-toolkit. These barcodes were then used to generate a FASTQ file input for Enrich2 to count variants. Variants with insertions, deletions, or multiple amino acid substitutions were excluded. Barcode counts were then collapsed to variant counts, retaining variants with a total frequency greater than 4 x 10^-5 across all bins (Supp. Fig. 2). For each replicate, an abundance score was calculated using a weighted average of variant frequency across bins (w 1 = 0.25, w 2 = 0.5, w 3 = 0.75, w 4 = 1) 22 . Scores were normalized to synonymous and nonsense distributions, excluding the top 20% of nonsense scores. Missense variant abundance scores ranged from -0.09 to 1.5.
Abundance classes were determined as in previous studies 22,30 . To discern between "WT-like" and "decreased" scores, we used a synonymous score threshold. This threshold was set at the 5th percentile of synonymous scores (0.856). Variants were classified as "WT-like" if their lower confidence interval exceeded the threshold, or as "possibly WT-like" if only their score surpassed the threshold. Additionally, an upper threshold at the 95th percentile of synonymous scores (1.14) was used to differentiate between "WT-like" and "increased" scores.
To distinguish between "decreased" and "nonsense-like" scores, we used a threshold at the 95th percentile of nonsense scores (0.265). Variants were categorized as "nonsense-like" if both their score and upper confidence interval were below the nonsense threshold, or as "possible nonsense-like" if only their score fell below the threshold. Out of a total of 8,480 variants, 316 were nonsense, 504 were synonymous, and 7,660 were missense variants. The missense variants were categorized into the following abundance classes: 2,590 WT-like, 612 possibly WT-like, 3,146 decreased, 437 possibly decreased, 340 possibly nonsense-like, and 708 nonsense-like.

VAMP-seq internal validation with individual variants
We cloned 11 variants using IVA cloning 72 site directed mutagenesis into the VAMP-seq recombination vector (attB-CYP2C19-eGFP-IRES-mCherry) via primers listed in Table S3 (HB049 through HB073, GB143, and GB144).  Fig. 13). To use multidms, we calculated abundance scores at the level of individual barcodes such that the model had to infer expected mutation effects of the CYP2C9 variants (β), while jointly identifying if (and by how much) the respective mutation effects may have shifted in the CYP2C19 background (Δ). The input data consisted of barcode-level abundance scores from two biological replicate DMS experiments for both homologs. We then grouped these four raw datasets into two independent training datasets, each of which consisted of one replicate from each CYP2C19 and CYP2C9. We then separately fit multidms models to the training sets and found the resulting fits to have well-correlated parameter sets (Supp. Fig. 6-7).
The models were trained to predict abundance scores in CYP2C9 and shifts in scores in CYP2C19 relative to CYP2C9, minimizing the difference in predicted and experimentally measured scores, as quantified by a Huber loss function. A key feature of the multidms modeling approach is the application of a lasso (L1) penalty to the shift parameters during the fitting procedure, effectively driving them to zero unless they are strongly supported by the data.
To determine a reasonable penalty coefficient (λ) for the lasso, we compared the results of seven model fits -each using different coefficients ranging from λ=0, to λ=5e-4. This λ sweep was performed in duplicate using each of the distinct training sets, for a total of 14 model fits evaluated. Supp. Fig. 8 shows summary statistics of the model fits. This figure emphasizes the accuracy-simplicity tradeoff for any given value of λ. For example, as λ increases, we observe higher overall training set loss, but also increased correlation between replicate parameters values. We selected a penalty coefficient of λ = 1e-5 as a reasonable value to balance model accuracy with shift sparsity. The shifts in the main text are then reported as the averaged values between the replicate model fits, each using the selected penalty coefficient. See https://github.com/matsengrp/CYP-multidms for the code used for this analysis.
We identified positions whose mean shift value was significantly different from the distribution of all shift values using a randomization test. To generate a null distribution, we randomly sampled 10 shift values, the average number of abundance scores per position, and calculated the mean of the shifts. This procedure was repeated 100,000 times. We calculated p-values for each position by counting the number of randomly generated mean shifts more extreme than the position mean and dividing it by 100,000, the total number of randomly generated shifts. P-values were then adjusted for repeated hypothesis testing using the Benjamini-Hochberg method 73 with a false discovery rate of 5%. Positions with p-values less than 0.05 were considered significant.

Additional files:
Table S1: CYP2C19 library fluorescence activated cell sorting

Data and code availability:
The accession number for the sequencing data reported in this paper is NCBI GEO: GSE244489. Code and processed variant scores generated during this study are available at GitHub: https://github.com/FowlerLab/cyp2c19_2c9. multidms analyses and data are available at GitHub: https://github.com/matsengrp/CYP-multidms