Apremilast Pharmacogenomics in Russian Patients with Moderate-to-Severe and Severe Psoriasis

One of the target drugs for plaque psoriasis treatment is apremilast, which is a selective phosphodiesterase 4 (PDE4) inhibitor. In this study, 34 moderate-to-severe and severe plaque psoriasis patients from Russia were treated with apremilast for 26 weeks. This allowed us to observe the effectiveness of splitting patient cohorts based on clinical outcomes, which were assessed using the Psoriasis Area Severity Index (PASI). In total, 14 patients (41%) indicated having an advanced outcome with delta PASI 75 after treatment; 20 patients indicated having moderate or no effects. Genome variability was investigated using the Illumina Infinium Global Screening Array. Genome-wide analysis revealed apremilast therapy clinical outcome associations at three compact genome regions with undefined functions situated on chromosomes 2, 4, and 5, as well as on a single single-nucleotide polymorphism (SNP) on chromosome 23. Pre-selected SNP sets were associated with psoriasis vulgaris analysis, which was used to identify four SNP-associated targeted therapy efficiencies: IL1β (rs1143633), IL4 (IL13) (rs20541), IL23R (rs2201841), and TNFα (rs1800629) genes. Moreover, we showed that the use of the global polygenic risk score allowed for the prediction of onset psoriasis in Russians. Therefore, these results can serve as a starting point for creating a predictive model of apremilast therapy response in the targeted therapy of patients with psoriasis vulgaris.


Introduction
Psoriasis is a chronic inflammatory disease that affects 2-3% of the world's population [1,2]. The annual incidence in Russia is 65, and the prevalence is 234.8 per 100,000 people [3]. Plaque psoriasis (i.e., psoriasis vulgaris-PV) is the most common variant of psoriasis, characterized by erythematous scaly patches or plaques that commonly occur on extensor surfaces. However, it can also affect intertriginous areas, palms, soles, and nails. The condition is often associated with a wide range of co-morbidities such as obesity, psoriatic arthritis, inflammatory bowel disease, cardiovascular, and psychosocial conditions that impair quality of life [1].
PV has a strong genetic component, as was confirmed by twin and family studies in populations of European descent. Estimates of heritability ranged from 50% to 90% [4]. Linkage analysis and the first genome-wide association study (GWAS) identified several loci as risk factors for onset PV. Most associations were found in genomic regions related with the immune system, i.e., major histocompatibility complex, signal transduction, apoptosis, cytokines, and receptors [5][6][7][8][9][10][11]. Recent data from meta-analysis revealed more than 80 loci associated with PV, which struck a multifactor nature of PV [2].
Combining multiple loci with modest effects into a global genetic risk score (GRS) might improve the identification of persons who are at risk for common complex diseases. Chen et al. were the first to develop such evaluation criteria for PV predisposition [12]. The GRS captured considerably more risk than any single-nucleotide polymorphism (SNP) considered alone. The area under curve (AUC) was chosen as the basic evaluation criteria for GRS fitness. Several subsequent GRS studies were based on the increasing number of loci associated with PV. The AUC rose to 0.8225 using 14 SNPs in the Chinese population [13]. The AUC reached 0.789 using just 5 SNPs in the Polish population [14]. For populations of Caucasian origin, the AUC achieved 0.76 using 63 SNPs [13]. A recent study created a PV prediction GRS based on 38 highly associated SNP [15]. In this study, we used GRS to differentiate patients with severe and moderate-to severe PV from the general Russian population.
The choice of PV treatment is based on the severity of the disease, which is determined according to the extent of the diseases and involvement of sensitive areas such as the face, scalp, and genital areas. The severity is roughly divided into mild, moderate, and severe, but in reality, overlap exists. Generally speaking, for mild PV treatment, various topical preparations for moderate PV phototherapy and oral systemic medications, as well as for severe PV orals, systemic medications and newer biological agents are used [16].
The understanding of PV molecular pathogenesis leads to the creation of target therapy approaches [1]. One such target drug is apremilast, which is the only selective phosphodiesterase 4 (PDE4) inhibitor approved by the U.S. Food and Drug Administration (FDA) and European Commission. Further, it is the only selective PDE4 inhibitor approved in Russia to treat moderate-to-severe and severe plaque psoriasis [17]. PDE4 inhibition, specifically hydrolyze cyclic adenosine monophosphate, leads to the cAMP intracellular level elevation [18]. Thus, the pro-inflammatory response downregulates within the T helper 1, T helper 17, and interferon pathways, all of which are pivotal in PV pathogenesis [19]. This, in turn, downregulates the production of pro-inflammatory cytokines, such as TNF-α, IL-2, IL-8, IL-12, IL-23, and IFN-γ. As such, this suppresses the PV pathway. On the other hand, PDE4 inhibition was reported to suppress the production of the IL-10 anti-inflammatory mediator and to increase the production of IL-6, which has both pro-and anti-inflammatory characteristics [20]. The side effects of apremilast are caused by PDE4 inhibition, which affects other organs (i.e., the adipose tissue) that may lead to weight loss and gastrointestinal symptoms (i.e., nausea, vomiting, and diarrhea). This most commonly occurs in the first week of therapy [17,19].
Apremilast is a safe and efficacious drug that is well tolerated. It has been shown that 33% of patients receiving apremilast achieved 75% or greater reduction from baseline versus those assigned to the placebo arm at week 16, whereas 61% patients achieved 75% at week 32 versus those assigned to the placebo. Almost the same figures were obtained during the clinical characterization of patients enrolled in the present study [3,20]. Considering the drug's long-term application, the expense of apremilast therapy, and the differences in drug response, it is reasonable to find personalized approaches for apremilast administration.
As confirmed by many studies, genetic polymorphisms are important factors for individual variations in drug responses [21,22]. One pharmacogenomic approach is SNP association. Certain studies, however, have found different clinical outcomes [23]. The aim of the present study was to assess the influence of genomic variability on PV apremilast therapy outcomes. Since the nature of psoriasis is complicated and there is a wide spectra influence of cyclic adenosine monophosphate (cAMP) levels, we used whole genome SNP coverage analysis based on Illumina Infinium Global Screening Array-24 v2.0. To search genomic regions that contain SNPs associated with a multifactor disease, two common approaches were used: candidate SNPs and genome-wide association. In this study, we explored apremilast target therapy pharmacogenomics considering the differences between patient groups with distinct clinical outcomes for both genome-wide associations and pre-selected PV-associated SNP set.

Clinical Characteristics of Patients
Phenotype classification involved subjects to be directly inspected by a dermatology specialist. The recruitment of patients was based on the following inclusion criteria: individuals over the age of 18 with a clinical diagnosis of psoriasis vulgaris (including the presence of joint damage) in a moderate or severe form (PASI ≥ 10) and at least 6 months duration of the disease (Figure 1). The exclusion criteria were as follows: other forms of PV (e.g., erythrodermic, guttate, or pustular psoriasis); drug-induced PV (e.g., onset or current exacerbation of the disease due to beta-blockers, calcium channel blockers, or lithium treatment); previous use of any drugs that directly affected interleukins or their receptors within 90 days prior to the study; presence of skin diseases other than PV (such as eczema) that may interfere with clinical assessment; pregnancy or lactation (breastfeeding), where pregnancy is defined as a woman's condition after conception and before the termination of pregnancy, confirmed by a positive human chorionic gonadotropin blood test; somatic diseases (i.e., metabolic, hematological, neurological, endocrine, infectious diseases, and diseases of the liver, kidneys, lungs, heart, or gastrointestinal tract), which significantly reduced the patient's immunity; uncontrolled arterial hypertension; the presence of lymphoproliferative or any other malignant neoplasms [3]. The clinical assessment of PV severity was carried out via PASI (Psoriasis Area Severity Index), degree of induration, desquamation, erythema, duration of exacerbation, and the presence of complications (PA). The PASI score was calculated as follows: four anatomical regions (head, upper limbs, trunk, and lower limbs) were assessed according to several indicators (i.e., erythema, infiltration, and desquamation). The severity of skin damage was rated on a 5point scale: 0-absence, 1-mild, 2-moderate, 3-severe, 4-very severe. The percentage The PASI score was calculated as follows: four anatomical regions (head, upper limbs, trunk, and lower limbs) were assessed according to several indicators (i.e., erythema, infiltration, and desquamation). The severity of skin damage was rated on a 5-point scale: 0-absence, 1-mild, 2-moderate, 3-severe, 4-very severe. The percentage area affected by PV in a separate anatomical region was evaluated as a percentage of the total area of a given anatomical region. Afterward, it was assigned a numerical value in accordance with the degree of psoriatic lesion (from 0 to 6): (0)-0%, (1)-1-9%, (2)-10-29%, (3)-30-49%, (4)-50-69%, (5)-70-89%, (6)-90-100%. The PASI score for each body area was calculated by multiplying the sum of the severity by the area affected; the result was multiplied by the coefficient corresponding to the given anatomical region (head-0.1; trunk-0.3; upper limbs-0.2; lower limbs-0.4). Skin lesions were considered mild if PASI was <10, moderate if PASI was <20, and severe if PASI was ≥20.

Target Psoriasis Therapy
The monotherapy of the selective PDE4 inhibitor apremilast (Otezla ® , Celgene International Sarl, Boudry, Switzerland) was carried out for all patients according to the Russian Federal Clinical Recommendations on PV patient treatment. The drug was prescribed in the form of tablets at a dose of 30 mg, taken orally 2 times a day (morning and evening) with an interval of about 12 h. In accordance with the manufacturer's recommendations, during the first 5 days of therapy, the dose was titrated with a stepwise increase from 10 mg to 30 mg. A clinical assessment of PDE-4 inhibitor efficacy was carried out using delta PASI indices (in percent) at week 26 after the beginning of targeted therapy. In order to search for probable predictors of response to apremilast with high effectiveness of targeted therapy, two comparison groups were formed based on the therapy outcome: (delta PASI 75% values and higher; 14 patients) and with its insufficient effectiveness (delta PASI less than 75%; 20 patients).

Sample Collection
A total of 34 blood samples were collected from the ulnar vein of PV patients living in different regions of Russia. The DNA samples were extracted from venous-blood leukocytes using a QIAmp DNA Mini Kit (Qiagen, Hilden, Germany). The quality of DNA samples was assessed using agarose gel electrophoresis, and the DNA quantity was determined using Qubit 3.0. All samples obtained were confirmed to be suitable for genotyping.

Genotyping
Samples were prepared using Illumina Infinium Global Screening Array-24 v2.0 (Illumina, San Diego, CA, USA). Array scanning on the HiScan (Illumina) platform were conducted according to the Infinium HTS Assay Guide. Genotyping was created using Illumina Auto Convert Software v2.0.1 (Illumina, San Diego, CA, USA). An Autosomal Call Rate Threshold equal to 0.01 was selected. Genotyping quality was assessed using call rate, and it appeared to be higher than 0.99 in all analyzed samples. Allele designations followed the Genome Reference Consortium Human Build 38 (GRCh38) (https://www. ncbi.nlm.nih.gov/assembly/GCF_000001405.38).

Data Analyses
The data obtained were inputwith BEAGLE 5.0 software, which had a Haplotype Reference Consortium reference panel and was filtered with DR2 > 0.7. An association test was conducted using PLINK v1.9 (classical filters were applied: geno 0.05, maximum percentage of missed values of polymorphisms; mind 0.1, maximum percentage of missed sample values; hwe 1e-5, exclusion of the polymorphisms (p-value threshold was 1e* −5 ) deviated significantly from Hardy-Weinberg equilibrium; mac 5, exclusion of SNPs by MAF (minimal number of minor alleles in a cohort is five)). The list of SNPs with statistically significant distinctions between patient groups with moderate or effec-tive therapy outcomes had a coincidence probability less than 10 −5 and a represented in Supplementary Table S1. The assessment of allele frequency differences at the pre-selected SNP sets between groups was conducted using the χ2 probability test. Therein, we found that the consuming significance probability level was 0.05 (Table 1). The polygenic risk score was 38 SNP and was calculated according to Kisiel et al. [15]. AUC values were calculated with the R-package "pROC" software. The random Russian population sample of 15,776 individuals was consulted with Genotek Ltd. and accordingly genotyped to be used as the control group [24]. Table 1. The data obtained for 78 pre-selectedsingle-nucleotide polymorphisms(SNPs)  in patient groups with advanced and lowered target apremilast therapy clinical outcomes. SNP with statistically significant differences are labeled with a star *. The population was from European origin and were healthy. The"1000 Genomes" project is shown for reference.

Results and Discussion
It should be noted that the effectiveness of therapy in this study was slightly lower than in international clinical trials and real-world studies [61][62][63][64][65]. So, 44% of patients did not reach PASI 50 by 14 weeks, and by 26 weeks-47%. This is probably due to the fact that more than half (53%) of the patients included in the study had severe psoriasis, and also had a burdened history in the form of the duration of the disease and the ineffectiveness of previous therapy.
Before the pharmacogenomic analysis, we checked whether the cohort possessed genome peculiarities that differed from the general Russian population. To evaluate the difference, GRS with 38 SNP were calculated both in the patient cohort (34 patients with moderate-to-severe and severe psoriasis) and in the population control group according to Kiesel et al. [15]. Figure 2 indicates violine plots obtained for the following groups: AUC value assessment for the differentiation was 0.6821 (95% CI: 0.5816-0.7837) and significantly differed from AUC = 0.5 for the non-informative model (Mann-Whitney-Wilcoxon p-value was 0.000175). This difference level suggested the possibility of using GRS to determine PV predisposition. In the original Kiesel et al. study [15], the AUC value for wGRS reached 0.782. Therefore, the lower AUC in this study can be explained by our low cohort sample size or by the inclusion of individuals with early stage PV. The usefulness of GRS for assessing an individual's PV predisposition does not mean that the same GRS is suitable for finding differences in patient cohorts with different clinical outcome of apremilast treatment. Calculating GRS in these patient groups revealed no distribution differences, thus uncovering nonsignificant influence of the genomic regions associated with PV predisposition on target treatment outcomes (Mann-Whitney-Wilcoxon p-value was 0.8766. AUC was 0.5179). The association test-which based all genetic data obtained with the Illumina Infinium Global Screening Array from patient cohorts with different clinical outcomes of apremilast treatment-revealed genome regions associated with therapy effectiveness. The restricted list of associations, limited by the 10 −6 p-value, was provided in Table S1. In detail, 72 SNPs, situated on four chromosomes, were associated with differences in PV apremilast therapy outcomes. There were11 SNPs at chromosome 2, covering 37 kbp region; 54 SNPs at chromosome 4 covered 44 kbp region; 6 SNPs at chromosome 5 covered6.4 kbp region and the only SNP at X-chromosome. These regions situated at the genome loci with an unrecognized function were not situated at or nearby certain genes, yet SNP rs35084576 at X-chromosome was in the ARSF (aryl sulfatase F gene). At first glance, SNP flanking regions show no remarkable peculiarities; however, SNP rs371208912 on chromosome 5 possesses minimal p-value, as it flanks the region with mutations. Therefore, extended deletions can occur (Live Ref SNPs, dbSNPb154v2, NCBI, 2020) (https://www.ncbi.nlm.nih.gov/snp/docs/RefSNP_about/). This region might be considered as a mutation hot point; however, further studies are required to resolve this issue. According to the genome-wide association analysis, there were a few genome regions highly associated with different apremilast therapy clinical outcomes. Considering the low sample size of the patient cohort, it was hard to form a haplotype analysis within the regions due to the high level of possible misinterpretation. It was also difficult to speculate about the less significant genome regions because of the huge amount of associations found. However, group differences, confirmed with a probability of 10 −6 , led us to conclude that our data can be used for a pilot model to predict apremilast therapy outcomes.
To explore the genomic regions associated with apremilast clinical outcomes, 78 SNPs strongly associated with PV or psoriatic arthritis were selected from the available The association test-which based all genetic data obtained with the Illumina Infinium Global Screening Array from patient cohorts with different clinical outcomes of apremilast treatment-revealed genome regions associated with therapy effectiveness. The restricted list of associations, limited by the 10 −6 p-value, was provided in Table S1. In detail, 72 SNPs, situated on four chromosomes, were associated with differences in PV apremilast therapy outcomes. There were11 SNPs at chromosome 2, covering 37 kbp region; 54 SNPs at chromosome 4 covered 44 kbp region; 6 SNPs at chromosome 5 covered6.4 kbp region and the only SNP at X-chromosome. These regions situated at the genome loci with an unrecognized function were not situated at or nearby certain genes, yet SNP rs35084576 at X-chromosome was in the ARSF (aryl sulfatase F gene). At first glance, SNP flanking regions show no remarkable peculiarities; however, SNP rs371208912 on chromosome 5 possesses minimal p-value, as it flanks the region with mutations. Therefore, extended deletions can occur (Live Ref SNPs, dbSNPb154v2, NCBI, 2020) (https://www.ncbi.nlm. nih.gov/snp/docs/RefSNP_about/). This region might be considered as a mutation hot point; however, further studies are required to resolve this issue. According to the genomewide association analysis, there were a few genome regions highly associated with different apremilast therapy clinical outcomes. Considering the low sample size of the patient cohort, it was hard to form a haplotype analysis within the regions due to the high level of possible misinterpretation. It was also difficult to speculate about the less significant genome regions because of the huge amount of associations found. However, group differences, confirmed with a probability of 10 −6 , led us to conclude that our data can be used for a pilot model to predict apremilast therapy outcomes.
Cytokines IL-1β and TNFα are pro-inflammatory cytokines produced by activated macrophages [66]. IL23R encodes a subunit of the receptor required for IL23A signaling. This protein associates constitutively with JAK2 and binds to transcription activator STAT3, which is one of the key signaling molecules for PV [67]. IL13 encodes a cytokine produced by activated Th2 that is involved in the maturation and differentiation of B cells. IL13 downregulates macrophage activity and inhibits proinflammatory cytokine and chemokine production. Along with the psoriasis vulgaris association, IL1β SNP rs1143633 is associated to eczema, hay fever, and asthma [68]. IL13 SNP rs20541has also been found to be associated with allergies, hay fever, asthma, and eczema [69]. IL23R SNP rs2201841 is associated with Crohn's disease [70]. TNFα(-308) SNP rs1800629 is associated with several immune diseases, i.e., asthma, Crohn's disease, system lupus erithematosus, Mediterranean spotted fever, multiple sclerosis, nasal polypes, lymphoma, leprosy, chronic obstructive pulmonary disease, and heart disease (https://www.snpedia.com/index.php/Rs1800629). Thus, the SNPs associated with apremilast therapy outcomes were also associated with immune polygenic diseases [71]. Recent studies on immunometabolism (i.e., the specific function of metabolites in immune system regulation via Krebs cycle rewiring [72]) uncovered the interconnection of intercellular metabolites (e.g., succinate, lactic acid, itakonate) levels on IL-1β and TNFα production during inflammation, as well the influence of α-ketoglutarate on the macrophage polarization process by changing IL4(IL13) dependent gene transcription [73]. The authors considered anti-inflammation effects for several drugs used in PV and RA treatment, such as dimethyl fumarate, metformin, and methotrexate. Since apremilast is an anti-inflammation drug that downregulates pro-inflammatory cytokine production after cAMP levels increase, it may also influence inflammation by changing metabolite concentration. The clinical application of apremilast results in better PA outcomes than PV outcomes. The succinate level was elevated in synovial fluid from RA patients [72] and it may be hypothesized that apremilast affects immunometabolism. SNP associations with clinical outcomes revealed indirect influences of the drug on associated cytokine levels due to cell metabolism modification. According to Zaslona and O'Neill, "We are still in the pioneering phase of gathering information about the functions of specific metabolites in immunoregulation" [73]. This offers hope for progress and future development regarding apremilast's influence on PV immunopathogenesis and therapy outcomes. An alternative way to solve the issue is genomic transcription study. cAMP signaling associated with G-protein coupled receptors (GPCR-cAMP-PKA pathway) [74], including CREBs (cAMP-response element binding protein) regulation of gene transcription under apremilast action [75].

Conclusions
Although understanding genetic backgrounds isjust a part of psoriasis pathogenesis [2], the associations found in this study will pave the way for the future creation of a model used in apremilast therapy outcome prediction. Moreover, the validation of GRS created for 38 SNPs was found to be applicable for PV susceptibility in Russians. The only obstacle to creating this model was the low patient cohort sample size used in this study. The sample size should be expanded in the future.
Supplementary Materials: The following are available online at https://www.mdpi.com/2075-4 426/11/1/20/s1. Table S1: List of top SNPs that resulted from the association test of genome-wide pharmacogenomic studies on apremilast between patients groups with distant therapy outcomes.
Author Contributions: A.A.K., supervision, study design, and manuscript editing; A.E.K., patient curation and treatment, clinical indices and therapy outcome assessment, study conceptualization; O.G.A., obtaining of informed patient consent and patient blood samples, clinical assistance, clinical indices, and calculations; D.A.V., manuscript drafting and review, DNA extraction, and data validation; D.G.D., conceptualization and methodology, and manuscript editing; V.S.S., project administration and coordination; A.K., genotyping and validation; A.R. bioinformatics analysis; A.E., formal analysis, genotyping, and bioinformatic supervision; A.C., software and visualization; All authors have read and agreed to the published version of the manuscript. Informed Consent Statement: Written informed consent was provided by all sample donors prior to participation in the study.

Data Availability Statement:
The genotypes dataset obtained is available on inquire from corresponding author.

Acknowledgments:
We are grateful for the biological samples donated by the volunteers enrolled as the subjects of this study. We thank the reviewers for valuable comments and suggestions that allowed us to improve the manuscript.

Conflicts of Interest:
The authors declare no conflict of interest.