Next-Generation Sequencing of PTGS Genes Reveals an Increased Frequency of Non-synonymous Variants Among Patients With NSAID-Induced Liver Injury

Purpose: The etiopathogenesis of drug-induced liver injury (DILI) is still far from being elucidated. This study aims to the study of genetic variations in DILI, related to the drug target, and specifically in the genes coding for the cyclooxygenase enzymes. Methods: By using Next-generation Sequencing we analyzed the genes coding for COX enzymes (PTGS1 and PTGS2) in 113 individuals, 13 of which were patients with DILI caused by COX-inhibitors. Results: The key findings of the study are the increased frequency, among DILI patients, of SNPs causing alterations in transcription factor binding sites and non-synonymous PTGS gene variants, as compared to control subjects. Moreover, the association with non-synonymous SNPs was exclusive of DILI patients with late-onset (50 days or more) Pc < 0.001 as compared to DILI patients with early onset, or with control subjects. Conclusions: Our findings suggest an interaction of long-term exposure to COX inhibitors combined with functional variants of the COX enzymes in the risk of developing DILI. This is a novel observation that might have been overlooked by previous genetic studies on DILI because of the limited coverage of PTGS genes in exome chips.


BACKGROUND
Although drug-induced liver injury (DILI) is a rare adverse drug event, it is often life-threatening because of the risk of developing acute liver failure. The mechanisms underlying DILI risk are not well understood and hence, the search for biomarkers of DILI risk is a major research field that aims to identify markers that could be used as both proof of the mechanisms involved and of the risk factors that can be used for DILI prediction, as has already been done with many pharmacogenomics biomarkers (Lucena et al., 2008(Lucena et al., , 2010Agúndez, 2009;Agundez et al., 2009Agundez et al., , 2011Andrade et al., 2009;Robles-Diaz et al., 2016;Nicoletti et al., 2017). There are presently several independent hypotheses to explain idiosyncratic DILI, but none of these is able to explain all the circumstances in which DILI occurs.
Some genetic biomarkers for DILI either mechanisticallybased using a case-control strategy or with a GWAs/exome sequencing approaches have been identified [for a review, see (Robles-Diaz et al., 2016)]. However, the involvement of genetic changes in DILI risk (for instance HLA risk alleles) has been documented for only a few drugs (Kaliyaperumal et al., 2018). On the other hand, case-control genotyping studies, GWAS and exome sequencing have important limitations because only some SNPs are tested, and most of the target sequence is not checked. To overcome this problem, deep sequencing comprising whole genes is necessary.
In this study, we analyzed the potential effect of mutations in cyclooxygenase genes (PTGS1 and PTGS2) on DILI risk related to NSAIDs. From a mechanistic point of view, such a risk could be related to genetic alteration in the arachidonic acid pathways, which are closely related to inflammation. On the other hand, adverse drug events for drugs acting on the COX enzymes (that is, COX inhibitors) may be more likely if COX activity is altered because of genetic variations. For this reason, we analyzed patients who developed DILI after the administration of COX-inhibitors and healthy individuals who tolerated COX-inhibitors.

CASE PRESENTATION
Thirteen patients (8 women and 5 men) who experienced DILI caused by COX inhibitors and 100 individuals who tolerated COX-inhibitors at standard doses were included in this study. The culprit drug for DILI and clinical details of patients are shown in Table 1. Gender-matched control individuals who tolerated COX-inhibitors (62 women and 38 men) individuals were recruited among staff and medical students of the Hospitals and the Universities participating in this study. Individuals which were considered as healthy after medical examination, to exclude pre-existing disorders and history of adverse events after the use of COX-inhibitors, were asked to participate and over 95% of these agreed to do so. We selected consecutive control subjects matched with patients for drug exposure: Fifty control subjects who have received ibuprofen within the previous month to sample collection, 20 who received diclofenac, 10 indomethacin, 10 naproxen, and 10 rofecoxib. These frequencies match with the frequencies for the DILI patients, except that no control subject received nimesulide since this drug was discontinued from the Spanish market due to liver safety. Both patients and controls were Caucasian Spanish individuals. Written informed consent for participation in this case report was obtained from all Abbreviations: COX, Cyclooxygenase, prostaglandin-endoperoxide synthase; NSAID, Non-steroidal anti-inflammatory drug; DILI, Drug-induced liver injury; GWAS, Genome-wide association study; SNP, Single nucleotide polymorphism; HLA, Human leukocyte antigen; PTGS1, Prostaglandin-Endoperoxide Synthase 1; PTGS2, Prostaglandin-Endoperoxide Synthase 2. participants. The protocol for this study was in accordance with the Declaration of Helsinki and its subsequent revisions and was approved by the respective Ethics Committees of the participating Hospitals.

DESCRIPTION OF LABORATORY INVESTIGATIONS AND DIAGNOSTIC TESTS
To achieve complete gene capture, we sequenced all exons, intron-exon boundaries as well as the 5 ′ and 3 ′ flanking regions for both genes. Referred to the GRCh37 assembly of the human genome, the sequences studied were the following: PTGS1: Chromosome 9:125.131.159 to 125.158.017; PTGS2: Chromosome 1:186.640.825 to 186.651.605. Partially overlapping amplicons with a size lower than 400 bp were designed. A total of 62 CS1/CS2 tagged primer pairs were synthesized and used to amplify 113 DNA samples using the Access Array platform (Fluidigm). During amplification, samples were labeled with standard MID barcodes designed for the FLX454 sequencing system. After amplification and MID-labeling, individual amplicon libraries were analyzed using a Bioanalyzer 2100 (Agilent) and bioanalyzer traces were used to estimate the amplicon concentration for each sample. Samples were then pooled, and libraries were purified by SPRI using Ampure beads to remove all possible traces of small molecules, primers, primer-dimers, or any other contaminants. The pooled library was again quantified and titrated so that a final amount of 1.95E+10 molecules with an enrichment percentage of 7% was loaded on a Pico Titer Plate (Roche) for a 200-cycle titanium-based sequencing run, made on FLX-454 equipment. Reads were processed using an amplicon processing pipeline and sff files were used for further analyses. Coverage averaged around 50x for the whole project. Coverage for the SNPs identified (shown in Supplemental Table 1) was always over 50x. Sequencing reads were de-multiplexed and aligned using the Amplicon Variant Analyzer software v2.8 (Roche) so that reads for each particular sample-target region combination were analyzed in search of variants. Details of the amplification and sequencing primers are available in Supplemental Table 1.
The putative effect on the non-synonymous variants identified in silico was assessed by using the Sorting Tolerant form Intolerant (SIFT) and Polymorphism Phenotyping (PolyPhen) scores as shown in the 1,000 genomes website for every SNP, as well as the online application MutationAssessor (http:// mutationassessor.org/r3/).

RESULTS
The sequencing results (summarized in Table 2) reveal that PTGS genes are well conserved. Although dozens of PTGS1 and PTGS2 single nucleotide polymorphisms (SNPs) have been described to occur in Caucasian populations (see Agúndez et al., 2015), our findings show that most of these SNPs were not identified, or were extremely rare, in this cohort.  Frontiers in Genetics | www.frontiersin.org    Interestingly, most of the PTGS1 and PTGS2 SNPs included in the Illumina human exome chip or human core exome chip (Urban et al., 2012) are also absent in this study group. This raises doubts about the coverage of exome chips to identify genetic associations related to PTGS1 and PTGS2 genes.
In the whole population study, we identified 31 single nucleotide polymorphisms (SNPs) for PTGS1, including four non-synonymous SNPs. For PTGS2 we identified 31 SNPs including one non-synonymous. We observed an increased frequency of PTGS1 and PTGS2 mutations among DILI patients, as compared to that observed in control individuals. Most of the SNPs identified in patients were rare among control individuals and were rare also according to the 1,000 genomes database (as shown in Table 2). All patients but one (case 1 in Table 2) had mutations at the PTGS1 gene and all patients but one (case 5 in Table 2) had mutations at the PTGS2 gene. Table 3 summarizes the comparison of relevant SNPs across patients with late-onset DILI, the rest of DILI patients and control individuals.

DISCUSSION OF THE UNDERLYING PATHOPHYSIOLOGY AND THE NOVELTY OR SIGNIFICANCE OF THE CASE
The most remarkable findings in this study are the presence among DILI patients of SNPs causing alterations in transcription factor binding sites such as the PTGS1 SNP rs10306225 (Agundez et al., 2014), and the PTGS2 SNPs rs4648253, rs689466, and rs20417, as well as non-synonymous SNPs such as PTGS1 rs1236913 (W 8 R), rs3842787 (P 17 L), rs5787 (R 108 Q), rs3842792 (K 185 T), and PTGS2 rs3218622 (R 228 H). These missense variants are extremely rare among European individuals (Agúndez et al., 2015). The putative effects of the most relevant SNPs shown in Table 3 have been revised elsewhere (Agúndez et al., 2015). In brief, besides the rs10306225 SNP, which is a promoter variant that causes a modification in a CDX1 binding site (Agundez et al., 2014), the rest of SNPs are non-synonymous. According to functional predictions and functional analyses (reviewed in Agúndez et al., 2015) the SNPs rs1236913, rs3842787 have a little functional effect, although clinical associations for these SNPs with urticaria induced by NSAIDs (Cornejo-Garcia et al., 2012) and myocardial infarction/stroke (Lee et al., 2008;Lemaitre et al., 2009;Gao et al., 2014), respectively, have been proposed. The functional effect of the rs5787 SNP is unknown, although functional prediction suggests a mild functional impact (see Table 2), rs3842792 SNP is predicted as functional ( Table 2), but in vitro findings suggest reduced functionality (Lee et al., 2007), and no functional impact for the PTGS2 SNP rs3218622 has been described.
No particular association of missense SNPs with culprit drug, age, gender, clinical presentation, type of liver injury, and severity of the disease was identified. However, as shown in Table 1, there is heterogeneity in the duration of treatment before DILI onset. This heterogeneity, rather than being a weakness, is a strong point in this study because it allowed discriminating the frequencies of PTGS gene variations in DILI patients with late and short-term onset. All the five DILI patients with the longest times to DILI onset (50 or more days; patients n • 3, 8, 9, 10, 12 in Table 1) had missense variants, and no patient with shorter time to DILI onset had such missense variants. The intergroup comparison values for carriers of any nonsynonymous PTGS variants were as follows: Patients with late DILI onset (50 or more days) vs. the rest of DILI patients (P < 0.001). Patients with late DILI onset vs. control individuals (P < 0.001). By turn, no significant differences for carriers of non-synonymous PTGS variants were observed among patients with DILI onset shorter than 50 days and control subjects (P = 0.325). Haplotype analyses (Table 4), and linkage disequilibrium (LD) analyses (Supplemental Table 2), show that the risk is due to the presence of rare haplotypes (containing missense variants) in the group of patients with late-onset DILI, but it is not due to LD variations for these variants. The strong association observed in this report, although it is based in five cases only, suggests a relationship of non-synonymous PTGS gene variations with DILI onset after long-term NSAID therapy. This is a novel observation that has not been raised by previous studies. Although the putative role of PTGS gene variations has been explored using the Illumina human exome chip or human core exome chip, it is of note that chip coverage was very limited for PTGS genes (Urban et al., 2012). By turn, this study has complete coverage thus allowing the identification of, as yet, disregarded SNPs. Another relevant difference with most DILI genetic studies is that in this report we stratified patients according to the time to onset. It cannot be ruled out heterogeneity in the etiopathogenesis of DILI, and it is conceivable that the mechanisms involved in DILI with a late onset might be different from those involved in immediate or short-latency reactions. This study, albeit with the inherent limitations of statistical power that case reports have, reinforces the view that a complete gene coverage and a detailed phenotype stratification of DILI patients could be essential to gain strength in further genetic association studies.

AUTHOR CONTRIBUTIONS
ML and EG-M participated in the design of the study, in data acquisition, and in critical revision for important intellectual content. AD, MB, and RA participated in the analysis and interpretation of the data and critical revision for important intellectual content. JA participated in the conception, design, data analysis and interpretation, the drafting of the manuscript and critical revision for important intellectual content. All authors approved the final version of the manuscript and all agree to be accountable for all aspects of the work in ensuring that questions related to the accuracy or integrity of any part of the article are appropriately investigated and resolved.