Longrange PCR-based next-generation sequencing in pharmacokinetics and pharmacodynamics study of propofol among patients under general anaesthesia

The individual response of patients to propofol results from the influence of genetic factors. However, the state of knowledge in this matter still remains insufficient. The aim of our study was to determine genetic predictors of variable pharmacokinetics and pharmacodynamics of propofol within selected 9 genes coding for propofol biotransformation enzymes, receptors and transporters. Our studies are the first extensive pharmaocgenetics research of propofol using high throughput sequencing technology. After the design and optimization of long range PCR-based next-generation sequencing experiment, we screened promoter and coding sequences of all genes analyzed among 87 Polish patients undergoing general anaesthesia with propofol. Initially we found that two variants, c.516 G > T in the CYP2B6 gene and c.2677 T > G in the ABCB1 gene, significantly correlate with propofol’s metabolic profile, however after Bonferroni correction the P-values were not statistically significant. Our results suggest, that variants within the CYP2B6 and ABCB1 genes correlate stronger with propofol’s metabolic profile compared to other 7 genes. CYP2B6 and ABCB1 variants can play a potentially important role in response to this anaesthetic and they are promising object for further studies.


Results
Long range amplification and NGS. We optimized all twenty-seven LR-PCRs using three sets of polymerases and successfully amplified over 226 kb of the DNA for 87 patients anaesthetized with propofol ( Fig. 1). Full-length gels are presented in Figure S1.
All 87 DNA libraries with a mean fragment length between 300-500 bp were successfully sequenced on the Illumina MiSeq platform. On average 77.2% of the reads were aligned to the target regions. The mean depth of the reads was 166 (between 19-399) (Fig. 2). Four amplicons (numbered 7,9,10,21) showed a considerable lower reads depth (<30-fold) among an average of 66% of the patient group analyzed and were excluded from further analyses. Moreover, amplicon number 11 showed insufficient coverage of approximately 25% of patients. Finally, only those variants of which the locus was read minimum 30 times in all patients, were included in the further analysis.
We detected 1513 different sequence variations in total. They were annotated and filtered using SNP-nexus and Variant Studio (Illumina). From the whole amount of the variants detected we selected only those described in the Human Genome Mutation Database (HGMD ® Professional 2015.4) and nonsynonymous changes ( Table 1).
These changes were randomly verified with Sanger sequencing.
Pharmacokinetics and pharmacodynamics of propofol. The mean value of MRT was 82.5 min and the range was 8.0-504.1 min (SD = 107 min). Based on a percentile rank we determined three pharmacokinetic profiles of the propofol metabolism rate in our group of patients: rapid (MRT ≤ 30 min), intermediate (100 ≥ MRT > 30 min) and poor (MRT > 100 min) metabolizers, who constitute 27%, 48% and 22% of the group, respectively.
To specify the differences in the biotransformation pathway of propofol we analyzed the formation profile of two main metabolites: 4-hydroxypropofol (4-OHP) and propofol glucuronide (PG). Using a k-means clustering, we determined three significantly different (P = 0.00) groups of patients due to the formation of PG and 4-OHP (Fig. 3).
Profile 1 of PG formation was observed in 13 of patients, profile 2 and 3 in 23 and 44 individuals, respectively. Clearly, 4-OHP was produced to a minor extent, relative to PG, the second metabolite. Profile 1 of 4-OHP formation was reported in 14 of patients, profile 2 in 13 of the subjects and profile 3 was presented by 53 individuals (Fig. 3).
The pharmacodynamics of the anaesthetic were described by the awakening time of the patients after propofol anaesthesia. The mean awakening time was 12 min, the range was 0.5-45 min (SD = 7.8 min).
Correlation of genetic variants with PK and PD of propofol. The impact of genetic background on the pharmacokinetics of propofol was proved by correlation analysis of genetic changes found in genes coding for metabolizing and in transporting proteins of propofol with MRT and formation profiles of PG and 4-OHP. The statistical significance of distribution differences of selected variants in rapid, intermediate and poor metabolizers is shown in Table 2.
Two genetic changes initially have shown a correlation with the propofol metabolic profile, c.516 G > T (p.Q172H) in CYP2B6 gene and c.2677 T > G (p.S893A) in ABCB1 gene. Homozygotes c.516 T/T were more often (P = 0.046) presented in the group of rapid metabolizers, while heterozygotes c.2677 T/G were more frequently (P = 0.032) observed in patients presenting an intermediate and poor metabolism rate of the anaesthetic.  However, using Bonferroni correction, the P-values were not statistically significant (P = 0.14, P = 0.1 for c.516 G > T and c.2677 T > G, respectively).
Moreover, we found that mutation p.S400N (c.1199 C > T) in the ABCB1 gene has a substantial impact on the profile of PG formation in our anaesthetized patient group. Homozygotes c.1199 C/C determined the profile 2 of the PG formation rate. However, none of variants analyzed affected the rate of 4-OHP formation.
The pharmacodynamics of propofol was correlated with genetic changes located in gene coding for the adrenergic receptor. The statistical significance of the influence of ADRA1A gene variants for the recovery time after anaesthesia is presented in Table 3. These analyses revealed no correlation between the pharmacodynamics of propofol and selected variants in the ADRA1A gene.

Discussion
Personalized anaesthesia with propofol still remains a challenge for pharmacogenomics. Global studies indicate in particular the contribution of CYP2B6, CYP2C9 and UGT1A9 gene polymorphism to the variability of propofol pharmacokinetics and pharmacodynamics [6][7][8]23 . Considering the complexity of the propofol biotransformation pathway, transport and interactions with receptors, in this study we presented deep sequence analysis of 9   24 . In our study homozygous c.516 T/T was identified more often in the rapid metabolizers of propofol, but no statistically significant correlation was found, in the contrast to Kansaku et al. 7 and Mastrogianni et al. 8 , whose observations proved an association of a surprisingly higher propofol concentration in the serum with allele T 7,8 . On the other hand, Loryan et al. 23 and Khan et al. 25 did not find any correlation between polymorphism p.Q172H and the plasma concentration of propofol or its metabolites 23,25 . We also cannot share Murao et al. 's (2015) conclusions that allele T decreased the total required dose of propofol 6 . The results of these different studies are contradictory, however, they do confirm the importance of a particular gene and a defined locus.
Moreover, we observed that heterozygotes c.2677 T/G of substitution p.S893A (c.2677 T > G; CM033585) located in exon 21 of the ABCB1 gene were more frequent (but not significantly) in the group of intermediateand poor-propofol metabolizers. So far, there has been a lack of investigations that would analyse changes in the ABCB1 gene in the context of the reaction and the biotransformation rate of propofol. Nevertheless, the ABCB1 gene is considered to be one of the most important pharmacogenes and the impact of polymorphic variability is widely examined in the context of a number of therapies, especially of lymphocytic leukemia 26 . The change p.S893A resides in an important region between the membrane and the cytosolic N-terminal NBD in the part of TM10 of the ABCB1 protein 27 .
As evidence of the variable biotransformation process of propofol, in our experimental group we identified considerable interindividual differences in the formation of propofol metabolites, both PG and 4-OHP. Interestingly, as a causal genetic factor linked to interpatient variability in the formation rate of PG, we found the mutation p.S400N (c.1199 C > T) in the ABCB1 gene. Homozygous individuals c.1199 C/C revealed the profile 2 ("middle homogeneous") of the glucuronidation rate. This mutation (CM045736) is known as disruptive of the ABCB1 protein transport activity and therefore it can probably influence propofol disposition and anaesthesia efficiency 28 .
The concept for molecular analysis in our investigations consisted of applying an approach combining LR-PCR with NGS, which enabled us to conduct an in-depth analysis of over 20 Mb of the sequence in a single run. A great advantage of using amplicon libraries and such methodology, compared to commercial sets of probes, is the easy and inexpensive modification of targeted genetic fragments through changes at the stage of designing long-range PCR. So far, only a few studies exist worldwide, which used a similar solution successfully [29][30][31] . On the other hand, we are of course aware of the limitations of our study, such as the number of patients, possible interactions between propofol with other drugs and the analysis of NGS data. It was crucial for us to obtain sufficient coverage for the genomic regions analyzed, which we unfortunately failed to achieve for a few amplicons, therefore 4 of them, covering exons 3-6 of the SULT1A1 gene, promoter and exons 1-8 of the GABRA1 gene and exons 6-15 of the ALB gene, were excluded from the final analysis. We suspect that they also can contain important changes for propofol action. To improve the experiment, perhaps the LR-PCR products of these 4 fragments should be added to the library in higher quantitative ratios in relation to other fragments.
Because we obtained a large amount of NGS data as detected variants, we decided to limit our analysis in the first stage only to the genetic changes recognized as most important from the pharmacogenetic point of view (mutations described in the HGMD database and as yet unknown amino acid changes). Nevertheless, our investigations are the first such a wide genetic analysis of the pharmacokinetics and pharmacodynamics of propofol in  the world and may potentially allow genes to be identified that may play an important role in the interindividual variability of its anaesthetic action.

Materials and Methods
Patients and samples. Eighty-seven Polish patients (34 female and 53 men), Caucasians undergoing general anaesthesia with propofol (10 mg/mL propofol injectable emulsion, Diprivan; AstraZeneca, Macclesfield, UK) before minor laryngological surgery at the Regional Hospital in Poznan, Poland, were enrolled in this study. All participants provided written informed consent. The study was approved by the local Ethics Committee of the University of Medical Sciences in Poznan, Poland (resolution no. 653/09) and all experiments were performed in accordance with relevant guidelines and regulations of this Committee. Detailed information about the patients: age, weight, height, dose of propofol, infusion time, awakening time and adverse effects were collected (Table 4). All patients represented I or II class on The American Society of Anesthesiologists (ASA) scale. Perioperative monitoring includes heart rate, blood pressure and saturation. Anaesthesia was induced with propofol (2 mg/ kg) followed by continuous infusion at a rate of 8 mg/kg/h, as described previously 23 . From each patient, 5 mL of peripheral blood for molecular analysis and 2 mL of blood samples were taken at 5 time points as follows: at the end of anaesthesia, 5, 10, 20 and 30 minutes after. These samples were collected for pharmacokinetic study.
LR-PCR amplification. The genomic DNA of each patient was isolated from the peripheral blood samples using the method with guanidine isothiocyanate (GTC). Detection of genetic variants in the coding sequence,  Amplicon libraries for NGS analysis were prepared based on twenty-seven LR-PCR fragments. Pairs of primer sets were designed in reference to the human genomic sequence (GRCh37/hg19) using Primer Blast and Primer 3 Plus software. In total, over 226 kb of the DNA of each patient were amplified in fragments ranging between 1 and 18 kbp (Table S1). After optimization, LR-PCR reactions were carried out on an Applied Biosystems 2720 Thermal Cycler  (Table 5). Generally, most of the fragments were amplified using the GoTaq ® Long PCR Master Mix (Promega) polymerase in accordance with the standard 2-step or 3-step PCR program.

NGS Library Preparation.
After visualization on the agarose gel, the LR-PCR fragments obtained were purified with Exo-Sap (Affymetrix, Santa Clara, USA) before they were quantified using Qubit dsDNA BR Assay System (Invitrogen, Carlsbad, USA) on the Qubit ® 2.0 Fluorometer (Thermo Fisher, Waltham, USA).
In the next step, 27 amplicons of each patient were pooled in equimolar ratios. According to the manufacturer's protocol, 1 ng of the pooled DNA fragments was subjected to the Nextera XT procedure (Illumina) using transposome technology. Finally, using Nextera XT DNA Sample Preparation Kit (Illumina) and Nextera ® XT Index Kit (96) (Illumina), we obtained eighty-seven both-side indexed DNA libraries ready for high-throughput sequencing. Quality control of the libraries was performed on TapeStation 2200 Instrument (Agilent, Santa Clara, USA) using an HS D1000 ScreenTape (Agilent). Normalization of all libraries was carried out with magnetic beads, according to the Nextera XT procedure.
Sequencing and bioinformatic analysis. Sequencing on the Illumina MiSeq System was performed as paired-end targeted resequencing using a MiSeq Reagent Kit v2 (300 cycle) (Illumina). Obtained reads were mapped to the reference DNA sequence (GRCh37/hg19) by a Burrows-Wheeler Aligner (BWA, version v0.7.10-r789) algorithm with default settings. Then, variants located in regions involved in the manifest were inspected using IGV (Integrative Genomics Viewer, version v2.3.59). The changes detected were then filtered using a set of criteria: GQX (genotyping quality) ≥ 30, read depth ≥ 30 and heterozygous read ratio ≥ 35%. Variants with potential pharmacokinetic significance included in Table 1 were randomly confirmed using Sanger sequencing. Amplicons were purified using shrimp alkaline phosphatase and exonuclease I, following manufacturer's instructions. Direct sequencing was performed using BigDye Terminator v3.1 Cycle Sequencing Kit (Thermo Fisher Scientific) on the Applied Biosystems 3500 and Series Genetic Analyzers. Primer sequences and PCR conditions are available upon request. HPLC measurements. The propofol concentration, as described previously 24 , as well as its two main metabolites: 4-hydroxypropofol and propofol glucuronide (Fluorochem Ltd., UK) in plasma retrieved from patients at five different time points was measured using the HPLC/UV (P580A, Dionex, Germany) system coupled to a fluorescence (RF2000, Dionex) detector. Plasma samples (150 μL) were mixed with 150 μL of 2 M trichloroacetic acid (TCA) and certifugated at 10,000 g for 10 min. An aliquot of the supernatant was injected onto an analytical    Pharmacokinetic and pharmacodynamic parameters of propofol. The pharmacokinetics of propofol was characterized by the mean removal time (MRT), determined using a PKSolver tool based on HPLC measurements of the propofol concentration, dose of anaesthetic and infusion time 32 . Moreover, for both metabolites (propofol glucuronide and 4-hydroxypropofol) different profiles of the formation rate were defined by clustering analysis using k-means. Awakening time, characterized as the mean time of eye opening, first breath and orientation were used as pharmacodynamic parameter of propofol anaesthesia. 4,5,8,9,11,16,17,19,21,22,23,24,25