Interindividual Variability in Lymphocyte Stimulation and Transcriptomic Response Predicts Mycophenolic Acid Sensitivity in Healthy Volunteers.

Mycophenolic acid (MPA) is an immunosuppressant commonly used to prevent renal transplant rejection and treat glomerulonephritis. MPA inhibits IMPDH2 within stimulated lymphocytes, reducing guanosine synthesis. Despite the widespread use of MPA, interindividual variability in response remains with rates of allograft rejection up to 15% and approximately half of individuals fail to achieve complete remission to lupus nephritis. We sought to identify contributors to interindividual variability in MPA response, hypothesizing that the HPRT1 salvage guanosine synthesis contributes to variability. MPA sensitivity was measured in 40 healthy individuals using an ex vivo lymphocyte viability assay. Measurement of candidate gene expression (n ± 40) and single-cell RNA-sequencing (n ± 6) in lymphocytes was performed at baseline, poststimulation, and post-MPA treatment. After stimulation, HPRT1 expression was 2.1-fold higher in resistant individuals compared with sensitive individuals (P ± 0.049). Knockdown of HPRT1 increased MPA sensitivity (12%; P ± 0.003), consistent with higher expression levels in resistant individuals. Sensitive individuals had higher IMPDH2 expression and 132% greater stimulation. In lymphocyte subpopulations, differentially expressed genes between sensitive and resistant individuals included KLF2 and LTB. Knockdown of KLF2 and LTB aligned with the predicted direction of effect on proliferation. In sensitive individuals, more frequent receptor-ligand interactions were observed after stimulation (P ± 0.0004), but fewer interactions remained after MPA treatment (P ± 0.0014). These data identify a polygenic transcriptomic signature in lymphocyte subpopulations predictive of MPA response. The degree of lymphocyte stimulation, HPRT1, KLF2, and LTB expression may serve as markers of MPA efficacy.

Mycophenolic acid (MPA) is a first‐line immunosuppressive agent used to prevent allogeneic transplant rejection and is prescribed to > 93% of all renal transplant recipients. Although originally approved for prevention of renal transplant rejection, it is now used for a spectrum of solid‐organ transplant rejection prophylaxis. Due to its pharmacodynamic effects and favorable toxicity profile, MPA has also become first‐line therapy for several autoimmune diseases, including lupus nephritis. , Even after the introduction of modern immunosuppressive agents, such as MPA, 1‐year rejection rates range from 8–15%. , Furthermore, significant interindividual variability in MPA response remains in the treatment of lupus nephritis, as rates of complete response are still only 45–60%.

WHAT IS THE CURRENT KNOWLEDGE ON THE TOPIC?

☑ Mycophenolic acid (MPA) inhibits the enzyme IMPDH2, selectively reducing lymphocyte proliferation in the absence of the alternative salvage purine synthesis pathway. Despite the widespread use of MPA, there is a clinical need to reduce solid‐organ transplant rejection and improve rates of complete therapeutic response in autoimmunity.

WHAT QUESTION DID THIS STUDY ADDRESS?

☑ We explore the transcriptomic contributors to interindividual variability in MPA response. We test the hypothesis that the salvage purine synthesis pathway contributes to this variability.

HOW DOES THIS STUDY ADD TO OUR KNOWLEDGE?

☑ This study identifies a polygenic transcriptomic signature in lymphocyte subpopulations predictive of MPA response. It supports a contribution from the salvage pathway and HPRT1 gene expression to MPA response.

HOW MIGHT THIS CHANGE CLINICAL PHARMACOLOGY OR TRANSLATIONAL SCIENCE?

☑ The predictive markers uncovered may identify individuals at risk for MPA resistance, prompting clinicians to target higher trough concentrations or select alternative immunosuppression strategies.

Mycophenolic mofetil undergoes rapid presystemic bioactivation to mycophenolic acid by carboxylesterases following oral administration and is primarily excreted in the urine after glucuronidation. , MPA inhibits a rate limiting enzyme in the de novo purine synthesis pathway, inosine‐5′‐monophosphate dehydrogenase (IMPDH), to deplete the guanosine pool within lymphocytes and reduce proliferation (Figure  ). , There are two IMPDH isoforms. IMPDH1 is expressed in many cell types, but IMPDH2 is expressed predominantly in stimulated lymphocytes. MPA is a more potent inhibitor of IMPDH2, therefore, exhibiting greater specificity for inhibition of stimulated lymphocyte proliferation. , In addition, other cell types maintain an active salvage purine synthesis pathway catalyzed by hypoxanthine‐guanine phosphoribosyl transferase (HGPRT). HGPRT allows for the production of guanosine independent of IMPDH. Canonically, the absence of this pathway in stimulated lymphocytes is thought to contribute to the specificity of MPA, as well as its favorable side effect profile. , Therapeutic drug monitoring of MPA has been suggested and genetic variants within genes involved in the uptake and metabolism of MPA (SLCO1B1/3, ABCC2, and UGTs) as well as its IMPDH targets have been shown to be associated with MPA pharmacokinetics and response. , , , , , , , However, in this study, we examine the noncanonical salvage pathway as a contributor to MPA resistance.

Figure 1

De novo and salvage purine synthesis pathway. A schematic of the purine synthesis pathway highlighting the de novo and salvage pathways. Mycophenolic acid is a noncompetitive, reversible inhibitor of inosine‐5′‐monophosphate dehydrogenase (IMPDH) within the de novo pathway. Hypoxanthine‐guanine phosphoribosyl transferase (HGPRT) is the enzyme responsible for the alternative salvage purine synthesis pathway. PRPP, phosphoribosyl pyrophosphate synthase.

In this study, we conducted both hypothesis‐driven and unbiased assessments of the contributors to MPA sensitivity in lymphocytes. We hypothesize that part of the interindividual variability in MPA response is attributed to the HGPRT salvage pathway as a compensatory mechanism for guanosine synthesis. We tested this hypothesis in healthy volunteers by measuring purine synthesis gene expression in the de novo and salvage pathways and assaying MPA dose response in an ex vivo lymphocyte assay. Subsequently, we used single‐cell sequencing of individuals exhibiting either high resistance or sensitivity to MPA in order to better ascertain additional transcriptomic contributors to the variability in MPA response. The relationship of candidate differentially expressed genes identified from these analyses to MPA sensitivity was queried using siRNA knockdown.

Methods

Healthy volunteers

We enrolled and obtained informed consent from 40 adult men and women with no medical history of autoimmune diseases or transplant (Table  ). Blood samples were collected in EDTA tubes (BD Vacutainer, Franklin Lakes, NJ) and immediately processed for the lymphocyte viability assay. Individuals were selected for single‐cell sequencing and siRNA knockdown based on having the highest or lowest lymphocyte viability ratios and their availability for follow‐up. This study was approved by the Indiana University School of Medicine Institutional Review Board (Approval #1603179330).

Lymphocyte viability assay

Blood samples were diluted in 3 parts 1× DPBS (Corning Cellgro, Manassas, VA). Ficoll‐Paque PLUS (GE Healthcare, Pittsburgh, PA) was added carefully to the diluted blood and centrifuged at 600 g for 30 minutes. The buffy coat layer was removed and centrifuged at 800 g for 10 minutes. The supernatant was discarded and the cell pellet was resuspended in 10 mL RPMI media (Corning Cellgro) supplemented with 10% fetal bovine serum containing bovine serum albumin (HyClone FBS; GE Healthcare) and 1% penicillin/streptomycin (MP Biomedicals, Solon, OH) with 50 µL of Lectin from Phaseolus vulgaris (PHA‐P) dissolved in phosphate buffered saline (5 µg/mL final concentration; Sigma Aldrich, St. Louis, MO). Cells were incubated at 37°C 5% CO2 for 24 hours. After 24 hours, cells suspended in media were removed and treated with or without antibodies to CD3/CD28 (25 µL per 1 mL of media; ImmunoCult Human CD3/CD28 T Cell Activator, STEMCELL Technologies, Cambridge, MA) and plated at 1,000 cells per well (100 µL total volume) in white‐bottom 96‐well plates for 48 hours. Stimulated and unstimulated cells were treated with 0, 0.01, 0.1, 0.5, 1, 2.5, and 5 µg/mL final concentration of mycophenolic acid (DOT Scientific, Burton, MI) initially dissolved in DMSO (1 µg/µL) in triplicates for 24 hours. The Cell‐Titer Glo Luminescent Cell Viability Assay (Promega, Madison, WI) was performed according to the manufacturer's protocol. Triplicate luciferase values were averaged and ratios were calculated using each individual's baseline value. The degree of stimulation was calculated as the lymphocyte viability ratio between stimulated and unstimulated lymphocytes in the absence of MPA.

siRNA knockdown

We performed siRNA knockdown for HPRT1, LTB, CCL4, and KLF2 by electroporation, 6 hours prior to the MPA treatment of our lymphocyte viability assay workflow. Lymphocytes were treated with 300 nM pooled HPRT1 Stealth siRNAs (HPRT1HSS105004, HPRT1HSS105005, and HPRT1HSS105006), Silencer Select siRNAs CCL4 (s12574), KLF2 (s20269, s20270, and s20271), LTB (s8311, s8312, and s194597), or scramble control (Negative Control Median GC Duplex; Invitrogen, Carlsbad, CA) using the Amaxa P3 Primary Cell 4D‐Nucleofector X Kit S (Lonza, Walkersville, MD). The E0‐115 method was run on the 4D Nucleofector Core Unit (Lonza). Lymphocytes were maintained in T25 flasks prior to electroporation and plated into 96‐well and 6‐well plates after a 6‐hour electroporation recovery period.

Measurement of gene expression

Gene expression of IMPDH2, HPRT1, and GAPDH was measured at baseline (after Ficoll Isolation), post‐anti‐CD3/CD28 stimulation with and without MPA treatment. Gene expression was also measured to confirm HPRT1, CCL4, LTB, and KLF2 siRNA knockdown. In conjunction with the lymphocyte viability assay procedure, stimulated lymphocytes were plated in 6‐well plates (2 mL per well). RNA was isolated using the miRNeasy Mini Kit (Qiagen, Hilden, Germany) and cDNA was made using the QuantiTect Reverse Transcription Kit (Qiagen) according to the manufacturer's protocol using the 2720 Thermal Cycler (Applied Biosystems, Waltham, MA). HPRT1 (Hs02800695_m1), IMPDH2 (Hs00168418_m1), CCL4 (Hs99999148_m1), and LTB (Hs00242739_m1) gene expression were measured using TaqMan Gene Expression Assays (Life Technologies, Foster City, CA) with GAPDH (Hs02786624_g1) as the control gene with technical triplicates. We measured KLF2 (Forward Primer: 5′‐TGCGGCAAGACCTACACCAA‐3′, Reverse Primer: 5′‐AAATACCAGTCACAGTTTGGGAGG‐3′) using previously published primers with GAPDH (Forward Primer: 5′‐TGCACCACCAACTGCTTAGC‐3′, Reverse Primer: 5′‐GGCATGGACTGTGGTCATGAG‐3′) as the control gene using SYBR Green (iTaq Universal SYBR Green Supermix; Bio‐Rad, Hercules, CA) with technical triplicates according to the manufacturer's protocol. All gene expression was measured on the ViiA 7 Real‐Time PCR System (Applied Biosystems).

Single‐cell RNA sequencing

The Center for Molecular Genetics at Indiana University School of Medicine performed single cell 3′ RNA sequencing at baseline, poststimulation with and without MPA treatment using the Chromium single cell system version 2 (10× Genomics, San Francisco, CA) and the NovaSeq6000 sequencer (Illumina, San Diego, CA). Sequencing was performed in three sensitive and three resistant individuals based on lymphocyte viability ratios. CellRanger 2.1.0 was utilized to process the raw sequence data generated. Briefly, CellRanger used bcl2fastq to demultiplex raw base sequence calls generated from the sequencer into sample‐specific FASTQ files. The FASTQ files were then aligned to the human reference genome GRCh38 with RNAseq aligner STAR. The aligned reads were traced back to individual cells and the gene expression level of individual genes were quantified based on the number of unique molecular indices detected in each cell. The filtered gene‐cell barcode matrices generated by CellRanger were used for further analysis.

Quality control checks led to the inclusion of cells between 400 and 3,000 genes expressed, < 15,000 counts of RNA and between 1% and 8% of mitochondrial gene content for the baseline samples. For stimulated samples, cells that expressed > 1,800 genes, < 100,000 counts of RNA, and < 8% mitochondrial gene content were included. In total, 120,194 cells were retained for downstream analysis. We performed the standard preprocessing, feature selection, dimension reduction, identification of anchors between samples, and integration for reference assembly according to the Seurat version 3 anchoring method. , We chose a dimensionality of 80 principal components for the integration of our 18 samples and a resolution of 0.8 for clustering analysis. We manually classified clusters based on common gene expression markers identified in the literature for different lymphocyte subpopulations. Single‐cell sequencing data is available on the Gene Expression Omnibus database (Accession #GSE141026).

Statistical analyses

K‐means clustering using the Hartigan–Wong method was used to classify individuals as resistant, moderate, or sensitive to MPA using the lymphocyte viability data and age. This classification was used for all downstream analyses. Kruskal–Wallis followed by Dunn's post hoc test with a Benjamini–Hochberg correction was performed on lymphocyte viability and gene expression data (d CT values). Gene expression data was normalized to GAPDH using the ΔΔCT method. A paired t‐test was performed between siRNA treatment of each gene and the scramble control. Spearman's correlation and asymptotic P values were calculated between lymphocyte viability and gene expression features. For single‐cell sequencing data, a χ2 test was conducted to compare lymphocyte subpopulation proportions. Differential gene expression was identified using a Wilcoxon Rank Sum test with Bonferroni adjusted P values. A P value < 0.05 was considered statistically significant. The differentially expressed genes were investigated for enriched pathways in Gene Ontology, Kegg, and Reactome databases using ClusterProfiler and ReactomePA. , All previously mentioned data analyses were performed using R 3.5.2. Enriched receptor/ligand interactions were identified using CellPhoneDB using genes expressed in > 10% of the cells in a cluster and 1,000 statistical iterations in Python version 3.7.

Results

Interindividual variability in ex vivo lymphocyte proliferation among healthy volunteers

We measured changes in lymphocyte viability after MPA treatment in 40 healthy volunteers, as described in Figure  . The ex vivo assay measures the viability of lymphocytes after MPA treatment with concentrations ranging from 0 to 5 µg/mL. This information is used to calculate the degree of stimulation and the residual lymphocyte proliferation after MPA treatment. The assay was highly reproducible in the same individual over the course of 17 months (Figure  ). The lymphocyte viability ratios at each MPA concentration and age were used for k‐means clustering to classify individuals as resistant, moderate, or sensitive to MPA treatment (Figure  , Table  ). It is of note that no African Americans were classified as resistant to MPA treatment. At each MPA concentration, significant differences in lymphocyte viability were identified between the groups. The 1 µg/mL dose‐response vector contributed most to an individual's k‐means clustering. At this clinically relevant concentration of 1 µg/mL MPA, lymphocyte viability ranged from 41% (most sensitive) to 67% (most resistant).

Figure 2

Predictors of ex vivo lymphocyte proliferation response to stimulation and mycophenolic acid (MPA) treatment. (a) Ex vivo lymphocyte viability assay workflow. The degree of stimulation is the ratio between stimulated lymphocytes (A) and unstimulated lymphocytes (C) in the absence of MPA. Percent lymphocyte viability is the ratio between MPA‐treated lymphocytes (B) and the absence of MPA (A). The three treatments used for downstream analyses are in blue. Baseline, poststimulation (0 µg/mL MPA) and poststimulation with MPA (1 µg/mL MPA). (b,c) Lymphocyte viability data among 40 healthy individuals. (b) Resistant, moderate, and sensitive individuals classified using k‐means clustering. X‐axis: final MPA concentration (µg/mL), Y‐axis: lymphocyte viability ratio normalized to 0 µg/mL MPA. Error bars are ±SD. Significance based on Kruskal–Wallis test. (c) X‐axis: k‐means clustering MPA response classification. Y‐axis: degree of stimulation is the ratio between stimulated and unstimulated lymphocytes at 0 µg/mL MPA ± SD. Significance based on Dunn's Kruskal–Wallis test with Benjamini–Hochberg correction. (d) IMPDH2 and HPRT1 gene expression across treatments. X‐axis: k‐means clustering MPA response classification. Y‐axis: 40 – delta CT (threshold cycle values). Higher value means greater expression. Error bars are ±SD. Significance based on Dunn's Kruskal–Wallis test with Benjamini–Hochberg correction using d CT values. (e) Spearman correlation matrix among age, lymphocyte viability, and gene expression features. Only significant correlations are shown. (f) Lymphocyte viability curve after HPRT1 vs. scramble siRNA knockdown (n = 5). X‐axis: MPA concentration (µg/mL). Y‐axis: lymphocyte viability ratio normalized to 0 µg/mL MPA. Error bars are ± SEM. Significance based on a paired t‐test. *P < 0.05, **P < 0.01, ***P < 0.001, ****P < 0.0001.

Lymphocyte stimulation altered among stratified healthy volunteers

We define the degree of stimulation as the ratio between stimulated and unstimulated lymphocytes in the absence of MPA. Sensitive individuals had 48% (P = 0.02) and 132% (P = 0.0002) increase in stimulation when compared with moderate and resistant individuals, respectively (Figure  ). The degree of stimulation was inversely correlated with lymphocyte viability ratios after treatment with 1 µg/mL (r = −0.43; P = 0.006), 2.5 µg/mL (r = −0.60; P = 4.2e−5), and 5 µg/mL MPA (r = −0.62; P = 2.3e−5; Figure  ).

Altered gene expression of purine synthesis genes among stratified healthy volunteers

Gene expression measured by quantitative polymerase chain reaction revealed higher IMPDH2 levels in sensitive individuals compared with resistant individuals at baseline (5.2‐fold; P = 0.012), post‐anti‐CD3/CD28 stimulation (2.8‐fold; P = 0.015) and post‐MPA treatment (2.6‐fold; P = 0.001; Figure  ). Baseline (r = 0.45; P = 0.004) and post‐MPA treated (r = 0.55; P = 0.0002) IMPDH2 gene expression levels had a significant positive correlation with the degree of stimulation (Figure  ).

In contrast, higher HPRT1 levels were detected in resistant individuals after stimulation (2.1‐fold; P = 0.049; Figure  ). This difference in HPRT1 expression is not observed at baseline or after MPA treatment. Stimulated lymphocytes’ HPRT1 levels were positively correlated with the proliferation ratios after treatment with 5 µg/mL MPA (r = 0.38; P = 0.01; Figure  ).

HPRT1 knockdown increases sensitivity to MPA

To determine if HGPRT contributes to MPA dose response, we performed siRNA knockdown of HPRT1 prior to MPA treatment in a resistant individual. Compared with the scramble siRNA control, the siHPRT1 treated lymphocytes became more sensitive to MPA treatment with a 12% (P = 0.003) decrease in viability at 1 µg/mL MPA (Figure  ). On average, we achieved 90% (±4.4 SEM) HPRT1 knockdown in these lymphocytes.

Single‐cell sequencing and clustering analysis among resistant and sensitive individuals

Because HPRT1 knockdown did not fully explain the variability in MPA interindividual response, we conducted single‐cell sequencing to further elucidate the transcriptomic contributors to MPA sensitivity. Three sensitive (1 female, 2 male) and three resistant (2 female, 1 male) individuals were sequenced under three conditions: baseline, poststimulation, and post‐MPA treatment at 1 µg/mL. Unsupervised clustering led to the identification of 20 distinct cell types, which were annotated based on known expression markers (Figure  ). , , , Subclusters of CD8+ T lymphocytes, as well as Th1 and Th2 helper T cells, were defined by the level of expression of the most highly expressed gene in the cluster. The resulting UMAP and t‐SNE depict the 20 clusters identified with separation among various T cell, B cell, and monocyte populations (Figure  ).

Figure 3

Integration, clustering, and classification of 120,194 peripheral blood mononuclear cells across different treatments. (a) Dot plot highlighting selected genes used to classify unsupervised clustering of cells into known subpopulations. (b) Uniform Manifold Approximation Projection (UMAP) and (c) t‐Stochastic Neighbors Embedding (t‐SNE) depiction of 20 subpopulations. (d) Proportion of composite lymphocyte populations between sensitive and resistant individuals. X‐axis: baseline, anti‐CD3/CD28 stimulated (stimulation), and stimulation with 1 µg/mL mycophenolic acid (MPA; stimulation w/MPA) treated lymphocytes. Y‐axis: log2 ratio between proportions of cells classified within each composite in resistant and sensitive individuals. A χ2 test was performed to determine whether the distribution of the proportion of the composite lymphocyte populations were equal between sensitive and resistant individuals. *P < 0.001, **P < 0.0001. CD4+ T cell composite does not include TH0 cells or regulatory T cells.

Variability in the proportion of lymphocyte populations between sensitive and resistant individuals

Cell subpopulations were merged into one of seven clusters for downstream analyses: B cells, naïve CD4+ T cells, regulatory T cells (Tregs), CD4+ T cells (except naïve or Treg), CD8+ T cells, natural killer cells, and monocytes/dendritic cells. Monocytes and dendritic cells were only analyzed at baseline because too few cells remained after stimulation with or without MPA to appropriately analyze them. These cells were subject to negative selection during the stimulation protocol due to their physical properties and because their proliferation is not enhanced by antibodies to CD3 and CD28. Significant differences in cell type proportion were observed between sensitive and resistant individuals (Figure  ). After stimulation, resistant individuals had higher proportions of naïve CD4+ T cells, differentiated CD4+ T cells, Tregs, and B cells. In contrast, sensitive individuals had a higher proportion of cytotoxic T cells and natural killer cells. Similar direction of effects were observed under all three conditions. The proportion and direction of effect of cell types across all 20 populations were comparable to the merged cluster data (Figure  ).

Variability in the gene expression of lymphocyte populations between sensitive and resistant individuals

Differential expression was assessed between the sensitive and resistant individuals in six merged clusters (Figure  , Table  ) and across all individual clusters (Figures  , Table  ). At baseline, increased expression of the HBA1, HBA2, HBB, and KLF2 genes was identified in resistant individuals and LYZ expression was higher in sensitive individuals across multiple populations (Figure  , Figure  ). Poststimulation, IL2, CCL3, and CCL4 were upregulated in resistant individuals whereas LTB and many mitochondrial genes were upregulated in sensitive individuals (Figure  , Figure  ). Poststimulation and after MPA treatment, CCL4 and IL2 were upregulated in resistant individuals and GZMA and ISG15 were upregulated in sensitive individuals (Figure  , Figure  ). Aligning with our hypothesis‐driven exploration of HPRT1 expression, certain cell‐type clusters revealed increased expression of HPRT1 in resistant individuals, including TH0 cells after stimulation (the largest cluster in the resistant population) as well as B cells and CD4+ T cells after MPA treatment (Table  ).

Figure 4

Differential gene expression between sensitive and resistant individuals based on composite lymphocyte populations. Volcano plots of differential genes expressed among composite clusters between resistant and sensitive individuals at baseline (a), poststimulation without mycophenolic acid (MPA) treatment (b) and with MPA treatment (c). X‐axis: log2 fold change. Y‐axis: ‐log10 of adjusted P value. Statistical significance based on a Wilcoxon Rank Sum test with Bonferroni adjusted P values. A positive fold change means higher expression in sensitive individuals and a negative fold change means higher expression in resistant individuals. CD4+ T cell composite does not include TH0 cells or regulatory T cells.

We identified a number of enriched pathways in each cluster after stimulation (Figure  ) and after stimulation with MPA treatment (Figure  ). After stimulation, sensitive individuals were more enriched for receptor signaling and antigen processing pathways, whereas resistant individuals were enriched for protein translation pathways. After MPA treatment, protein translation was upregulated in sensitive individuals and cell cycle pathways were upregulated in resistant individuals. For example, the TNFR2 noncanonical NF‐kB pathway involving LTB, was upregulated in CD8+ T cells in the resistant group (P = 1.21E‐7). After MPA treatment, CD4+ T cells in the resistant group exhibited upregulation of the G1/S transition pathway (P = 1.8E−22). All pathway enrichment is provided in Table  .

Figure 5

Pathway enrichment among composite lymphocyte populations between sensitive and resistant individuals. Pathways enriched after stimulation (a–f) and after MPA treatment (g–l). X‐axis: descending order of enriched pathways by P value. Y‐axis: ‐log10 P value for enriched pathways using KEGG, Gene Ontology, and ReactomePA. (a/g) B cells. (b/h) TH0 cells. (c/i) CD4+ T cells. (d/j) T Regulatory cells. (e/k) CD8+ T cells. (f/l) Natural killer cells. Pathways upregulated in resistant individuals (pink/circles). Pathways upregulated in sensitive individuals (blue/triangles).

To validate the differential expression findings, we performed siRNA knockdown of three genes. CCL4 and KLF2 were predicted to increase MPA sensitivity and decrease lymphocyte viability after MPA treatment. LTB was predicted to decrease MPA sensitivity and increase lymphocyte viability after MPA treatment. CCL4 knockdown resulted in a 5% (P = 0.02) decrease in lymphocyte viability at 0.1 µg/mL MPA concentration but trended in the opposite direction at higher MPA concentrations (Figure  ). For KLF2, the expected decrease in lymphocyte viability was observed at 2.5 µg/mL (12%; P = 0.0008) and 5 µg/mL MPA concentrations (14%; P = 0.003; Figure  ). For siLTB treatment, the expected increase in lymphocyte viability was observed 1 µg/mL MPA (9%; P = 0.022; Figure  ). We confirmed a knockdown of 77% (±12 SEM) for CCL4, 40% (±10 SEM) for KLF2, and 72% (±8 SEM) for LTB compared with scramble control.

Figure 6

Validation of select differentially expressed genes within the single‐cell sequencing data. Lymphocyte viability after (a) CCL4, (b) KLF2, or (c) LTB (solid/pink line) vs. scramble (dashed/blue line) siRNA knockdown. X‐axis: mycophenolic acid (MPA) concentration (µg/mL). Y‐axis: lymphocyte viability ratio normalized to 0 µg/mL MPA. Significance based on a paired t‐test. Error bars are ± SEM. CCL4 (n = 5), KLF2 (n = 3), and LTB (n = 4). *P < 0.05, **P < 0.01, ***P < 0.001.

Differential receptor‐ligand communication among lymphocyte populations is found in sensitive and resistant individuals

Using CellPhoneDB to predict receptor ligand communication between populations of lymphocytes, there were no significant differences in cross talk at baseline. Among resistant individuals, 201 interactions were detected compared with 217 among sensitive individuals (Figure  ). However, after stimulation, there was significantly greater cross talk (i.e., proportion of interactions specific to resistant or sensitive individuals), among sensitive individuals as compared with the resistant individuals (P = 0.0004). The opposite direction of effect was observed after MPA treatment, where more receptor ligand interactions were maintained in the resistant group (P = 0.0014). Venn diagrams to show overlap in interactions between sensitive and resistant individuals within treatment groups are shown in Figure  and the list of interactions are provided in Table  .

Figure 7

Differential ligand‐receptor interactions among composite lymphocyte populations between sensitive and resistant individuals. (a) Circlize plots to depict receptor‐ligand communication between populations. Edges reflect the specific lymphocyte population. The axis around the edges reflect the number of interactions involving the population. (b) Venn diagrams depicting the shared and unique receptor‐ligand interactions between resistant and sensitive individuals. Statistical significance based on a χ2 test. CD4+ T cell composite does not include TH0 cells or Tregs. MPA, post‐mycophenolic acid treatment; R, resistant; S, sensitive; Stim, poststimulation.

Discussion

Through this multifaceted investigation, we identified transcriptomic markers of interindividual variability in MPA response. These markers help to elucidate the complex and polygenic nature of MPA's pharmacodynamic profile. In a hypothesis‐driven approach, we found that both HPRT1 expression and the degree of lymphocyte stimulation contribute to MPA response. Although the HGPRT salvage pathways is not thought to contribute to MPA response in lymphocytes, , we have shown that resistant individuals have higher HPRT1 gene expression after stimulation as well as increased MPA sensitivity after siHPRT1 knockdown. Through single cell sequencing, we narrowed the source of this higher expression to naïve TH0 cells, the largest and most disproportionate cluster in resistant individuals as compared with sensitive individuals. Further, after MPA treatment, both B cells and CD4+ T cells expressed higher HPRT1 in resistant individuals.

It is known that MPA has more potent cytostatic effects on stimulated lymphocytes due to their increase in IMPDH2. , Our study confirmed that IMPDH2 levels are higher in the sensitive individuals after stimulation and that these individuals are stimulated 132% more compared with the resistant individuals, aligning with the predicted pharmacodynamics of MPA. Of note, increased IMPDH2 after MPA treatment was unexpected, but may be attributed to a reflexive rise in expression after enzyme inhibition because sensitive individuals’ lymphocytes still showed a reduction in viability. In addition, 100% of African Americans and 75% of Hispanics were classified as moderate or sensitive to MPA in our population. Although a small sample size, this aligns with data from the Aspreva Lupus Management Study showing that African Americans and Hispanic groups had a higher response rate to mycophenolic mofetil than Asians and Caucasians for the treatment of lupus nephritis.

Despite the observed effects of HPRT1 and IMPDH2, these genes did not fully account for all of the variability observed between sensitive and resistant individuals. We observed up to a 26% absolute difference in lymphocyte viability between sensitive and resistant individuals. However, siRNA knockdown of HPRT1 led to a modest increase in lymphocyte viability by 12% at the 1 µg/mL treatment concentration. In an unbiased exploration, single‐cell sequencing identified additional candidate genes and pathways that may contribute to MPA efficacy across the lymphocyte populations. Again, the contribution of each individual gene was small, even when assayed by siRNA knockdown, with many candidate genes having predicted effects on MPA sensitivity that oppose each other. MPA response is polygenic and balancing the multitude of opposing factors to fully recapitulate drug sensitivity can be understood more directly in aggregate, through pathway and receptor‐ligand interaction enrichment. Indeed, we identified significant differences in enriched pathways, with various lymphocyte subpopulations in the sensitive individuals having more receptor‐ligand interactions and T cell receptor signaling pathways upregulated after stimulation and resistant individuals had more receptor‐ligand interactions and cell cycle pathways upregulated after MPA treatment.

A notable limitation of our study is the ex vivo nature of the lymphocyte viability assay conducted in healthy volunteers. However, this limitation is counterbalanced by the increased sensitivity to detect pharmacodynamic markers of MPA efficacy, independent of pharmacokinetic ones like concentration, exposure duration, or metabolism. The ex vivo assay represents a primary culture as cells are not frozen, passaged, or immortalized. The dose‐response curve for MPA reduction in lymphocyte viability in the assay aligns with clinically relevant concentrations. Typical trough concentrations in humans range from 1 µg/mL to above 3 µg/mL. , , Because several participants reached a plateau of viability above 50%, an alternate determinant of sensitivity was utilized beyond a half‐maximal viability concentration. Based on k‐means clustering, the strongest vector contributing to interindividual variability was the lymphocyte viability ratio at 1 µg/mL, which precedes the plateau phase of viability and is within the range of clinically relevant concentrations. As such, this concentration was selected for downstream transcriptomic single cell analysis. Additionally, previous literature has shown 100% inhibition of lymphocyte proliferation using the [3H]‐tritiated thymidine incorporation assay that allows the detection of newly synthesized strands of DNA. , Our study design measures cell viability through the amount of ATP present at the time of the luciferase assay. Because MPA is cytostatic and we do not measure cell viability at the 72 hour time point, we were able to detect lymphocytes that were present prior to MPA treatment; however, we believe our results are consistent with complete inhibition at the 5 µg/mL MPA concentration based on these previous studies and the observed plateau at the higher concentrations. In addition, we did not explore the impact of MPA treatment on IMPDH1 gene expression or genetic variation on MPA sensitivity due to small sample size. However, our gene expression data may account for genetic variations related to HPRT1 and IMPDH2 gene expression.

The network of genes uncovered that impact lymphocyte stimulation and MPA response is vast. The direction of effect and effect size of these genes vary across the subpopulations. This is to be expected since Tregs, B cells, and CD4+ T cells maintain different roles in proliferation. , The siRNA validation experiments provide additional confidence for the mechanistic contribution of a few genes to interindividual variability, albeit in the collective pool of all lymphocytes. However, it is not feasible to test all of the predicted interactions observed in our dataset. The differential receptor‐ligand interactions and their related pathways reveal a complexity of the MPA response phenotype, which cannot be captured by the knockdown of a select few genes. It is more likely that these genes, including HPRT1, could serve as representative markers of MPA response, prompting clinicians to target higher trough concentrations of MPA or select alternative immunosuppressive agents when warranted.

Our study was conducted in healthy volunteers. As this technology is adapted to patients with lupus nephritis or a renal transplant, an entirely different cell proliferation signature may be uncovered in those with autoimmune disease or chronic MPA exposure. Future directions may include examination of MPA response using HPRT1 expression or lymphocyte stimulation as biomarkers in these populations. Comparison of our single‐cell data in healthy volunteers to individuals with autoimmunity or immunosuppression may help uncover key pathways that explain the underlying biology of lymphocyte proliferation and MPA resistance. Such information will facilitate the translation of the uncovered markers to clinical care and allow us to better predict how an individual will respond to MPA therapy.

Funding

This work was supported by funding from the National Institutes of Health‐National Institute of Diabetes and Digestive and Kidney Diseases (NIH‐NIDDK) K08‐DK107864 (Eadon).

Conflict of Interest

The authors declared no conflict of interest.

Author Contributions

K.S.C., P.C.D., and M.T.E. wrote the manuscript. K.S.C., M.D.D., P.C.D., and M.T.E. designed the research. K.S.C., Y.‐H.C., D.J., N.A.K., and C.W. performed the research. K.S.C., R.M.F., H.G., P.C.D., and M.T.E. analyzed the data.

Figure S1

Click here for additional data file.

Figure S2

Click here for additional data file.

Figure S3

Click here for additional data file.

Figure S4

Click here for additional data file.

Figure S5

Click here for additional data file.

Table S1

Click here for additional data file.

Table S2

Click here for additional data file.

Table S3

Click here for additional data file.

Table S4

Click here for additional data file.