Genetic Variation among Pharmacogenes in the Sardinian Population

Pharmacogenetics (PGx) aims to identify the genetic factors that determine inter-individual differences in response to drug treatment maximizing efficacy while decreasing the risk of adverse events. Estimating the prevalence of PGx variants involved in drug response, is a critical preparatory step for large-scale implementation of a personalized medicine program in a target population. Here, we profiled pharmacogenetic variation in fourteen clinically relevant genes in a representative sample set of 1577 unrelated sequenced Sardinians, an ancient island population that accounts for genetic variation in Europe as a whole, and, at the same time is enriched in genetic variants that are very rare elsewhere. To this end, we used PGxPOP, a PGx allele caller based on the guidelines created by the Clinical Pharmacogenetics Implementation Consortium (CPIC), to identify the main phenotypes associated with the PGx alleles most represented in Sardinians. We estimated that 99.43% of Sardinian individuals might potentially respond atypically to at least one drug, that on average each individual is expected to have an abnormal response to about 17 drugs, and that for 27 drugs the fraction of the population at risk of atypical responses to therapy is more than 40%. Finally, we identified 174 pharmacogenetic variants for which the minor allele frequency was at least 10% higher among Sardinians as compared to other European populations, a fact that may contribute to substantial interpopulation variability in drug response phenotypes. This study provides baseline information for further large-scale pharmacogenomic investigations in the Sardinian population and underlines the importance of PGx characterization of diverse European populations, such as Sardinians.


Introduction
Drug treatments are characterized by substantial difference in terms of efficacy and/or safety in different patients. Adverse drug reactions (ADR), including allergic, pseudoallergic, and exaggerated pharmacological reactions to medications, are a relatively common result of drug treatment, accounting for at least 5% of hospital admissions, with an overall fatality of 0.15% and an annual cost of >500 M$, only for the UK National Health Service [1,2]. These data highlight the social and economic costs of ADRs and the urgent need to find effective strategies to ameliorate drug efficacy and reduce ADR.
The same drug, once other parameters are fixed, can have different therapeutic effects in different people due to causal genetic variants [3]. The analysis of the genetic variability modulating the individual's drug response (pharmacogenetics, PGx) has, thus, received great attention for its capacity to provide a new way to optimize drug therapies in terms of optimal dosing to improve drug efficacy and reduce toxicity risk [4]. As a result, a patient may receive the right drug at the right dose the first time they consult their doctors such that efficacy is guaranteed, and the risk of ADR is reduced. From a pharmaceutical point of view, PGx variants can influence pharmacokinetics and pharmacodynamics drugs, thus influencing dosing, formulation sensitivity and drug-hypersensitivity reactions.
An individual's drug response can be assessed through the identification, by genotyping arrays or sequencing, of well-characterized genetic variants and specific haplotypes in key genes implicated in drug processing. For example, the gene CYP2D6 is characterized by the presence of over 100 haplotypes, which share SNPs and include gene duplications and deletions, strongly influencing the metabolism and/or bioactivation of many clinically used drugs and, thus, determining a phenotype. In this example, phenotypes are assigned to haplotypes that contains specific and relevant SNPs to differentiate CYPD6 functions [5].
The interest in ameliorating drug efficacy, while reducing ADR, promotes the development of tools to properly analyze the correlation between variability in the genome and individual's drug response. For example, The Pharmacogenomics Knowledge Base (PharmGKB: http://www.pharmgkb.org (accessed on 22 March 2022) [6,7]) covers much information about pharmacogenomics and provides a convenient approach for researchers. The Pharmacogenetics of Membrane Transporters (PMT) database is another tool focused on the effect of genetic variation in the response to drugs that interact with membrane transport proteins [8,9].
Furthermore, the increasing availability of accurate classifications of pharmacogenetic variants and haplotypes, together with guidelines for their clinical translatability, allow analysis of the potential impact of pharmacogenetics programs in many populations for which large-scale genomic resources exist [10,11]. The analysis of the prevalence of PGx-risk variants in target populations, in combination with actual data on drug usage, make it possible to predict the proportion of the population for which genetics could lead therapy decision. Overall, the following axes could support a coordinate pharmacogenetic program in the European healthcare systems: (i) the analysis of PGx variant prevalence, (ii) the results of clinical trials evaluating patient outcomes and cost-effectiveness of PGx-markers [12] and (iii) outcomes of implementation strategies [13,14].
Sardinians, a population for which large-scale genomic data are available, is particularly well suited for genetic studies. Sardinians are the contemporary human population that has retained the highest degree of inheritance from early European farmers who lived in the Neolithic period along with significant ancestry from western hunter gatherers who lived in the late Paleolithic period [15,16]. This is due to founder effects during the initial settlement of the island and the scarcity of gene flow from other populations during later periods [15,16]. As a result of its past evolutionary history the Sardinians are now a reservoir of ancient European genetic variants that are currently very rare elsewhere and may have relevant clinical consequences [17][18][19][20]. Genetic factors and the distinct genetic structure of the Sardinians thus present an excellent opportunity to also look for new pharmacogenetic information.
Here, we profiled pharmacogenetic variation in fourteen clinically relevant genes in 1577 unrelated sequenced Sardinians. We used PGxPOP [21], a PGx allele caller, based on the guidelines created by the Clinical Pharmacogenetics Implementation Consortium (CPIC), to identify the main phenotypes associated with the PGx alleles most represented in Sardinians. We estimated that 99.43% of Sardinian individuals might potentially respond atypically to at least one drug, and that, on average, each individual is expected to have an abnormal response to about 17 drugs. Furthermore, we highlighted differences in haplotype and diplotype frequencies of star alleles as compared to other populations and estimated that for 27 drugs the fraction of the population at risk of atypical responses to therapy is more than 40%. These findings represent the foundation for further large-scale and more detailed pharmacogenomic investigations in Sardinia, and, at the same time, underline the importance of the pharmacogenomic characterization of ethnically diverse European populations, as exemplified by Sardinians.

Haplotype and Phenotype Calling
To identify clinically relevant pharmacogenetic variation in 14 important genes for which Clinical Pharmacogenetics Implementation Consortium (CPIC) has created detailed gene/drug clinical practice guidelines, we processed the genomic sequence data from our SardiNIA cohort (1577 individuals) with PGxPOP [21].
PGxPOP is a PGx matching engine that is based on PharmCAT and uses its PGx allele definitions to characterize PGx haplotype, diplotype and phenotype frequencies.
It extends the capabilities of PharmCAT by generating diplotypes from population scale datasets [22]. In the analyzed Sardinia cohort, 99.43% of the 1577 participants carried at least one diplotype associated with a predicted non-typical response phenotype across the 14 pharmacogenes analyzed (Table 1). Furthermore, for each participant, we were able to predict an average of about 4 phenotypes of non-typical drug response (Mean = 3.44, Min = 1, Max = 8, Figure 1). These numbers were in line with what has been observed in the UK Biobank cohort, where participants were previously reported to carry on average 3.7 nontypical response diplotypes for the 14 pharmacogenes, with 99.5% of participants carrying at least 1 nontypical drug response diplotype [21].  Although, in general, we observed some agreement with the observations in UKBB, we noticed some differences in terms of haplotype and diplotype frequencies. Largest absolute discordance in star allele frequencies (delta MAF ≥ 10%) between Sardinians and individuals with European ancestry of the UK Biobank cohort were observed for CYP2D6*1 and *119 (results that can be affected by the fact that structural variants were not called in this analysis), CYP4F2 alleles *1, *2 (this being extremely rare in Sardinians) and *2 + *3, VKORC1 alleles −1639A and −1639G, and for SLCO1B1 alleles *1A and *14 (Supplementary Table S3).
In addition, out of 133 diplotypes called by PGxPOP for the 14 genes analyzed, only 13 diplotypes (10%) had an absolute difference in frequency ranging from 10% to 25%, as compared to UKBB European individuals (Supplementary Table S4). We detected a total of 22 different phenotypes for the 14 pharmacogenes, but for CYP2D6 and UGT1A1 our analysis did not detect any non-typical drug response diplotype in the SardiNIA cohort (Table 2). To understand whether the unexpected result for CYP2D6 was due to a limitation of the SardiNIA genetic map in the CYP2D6 region, we performed an exploratory analysis with 65 deeply sequenced samples from the same cohort (mean coverage > 30×, data not shown here), where CYP2D6 star alleles were called with the Aldy tool [23]. The results showed that common CYP2D6 star alleles had comparable frequencies to those reported for the European samples in the UKBiobank cohort [21]. Based on these analyses, we could predict that, for at least 27 drugs, about one third of the population was at risk of an atypical response (Table 3) and that, on average, each individual was expected to have an abnormal response to about 17 drugs ( Table 1).
If we focused on the frequency of non-typical drug response phenotypes in the Sardinian cohort compared to the European sub-population of UK Biobank, the more frequent in Sardinians were "Decreased warfarin dose" for the VKORC1 gene (26% vs. 14%) and "Intermediate metabolism" for the CYP3A5 gene (19% vs. 13%). For 5 other non-typical drug response phenotypes, the difference in frequency was between 1 and 5%, while for all other non-typical phenotypes the frequency in the Sardinian population was always lower than in the European UK Biobank cohort (full details are available in Supplementary Table S5). Table 3. For each gene the percentage of individuals in the Sardinia cohort who were at risk of atypical response is reported. Only Gene/Drug pairs with CPIC evidence of level A were considered here (according to https://cpicpgx.org/genes-drugs/ (accessed on 16 August 2022)). Estimated frequency of non-typical response phenotypes in Sardinians (SARD) were compared with the corresponding values in UKBiobank European populations (UKBB-EUR) reported by McInnes et al. [16] and differences among the two values are reported in column "Delta freq." (See extended data in Supplementary Table S5). The p-value refers to pairwise comparison of phenotype frequencies in Sardinian and European populations. CYP2D6 predictions were omitted because frequencies were likely to be affected by copy number variants that are not considered here.

Gene
Related Drugs

Actionable Pharmacogenomic Variants in Sardinian Genomes
We extended the above analysis by considering the carrier status of clinically actionable variants (PharmGKB level 1A/1B) in each of the 1577 Sardinian unrelated individuals.
Among the 3073 clinical annotations in the PharmGKB database (accessed on 22 March 2022), 141 single-nucleotide variants (SNVs) and 132 haplotype variants had the highest level of evidence (1A/1B or 2A/2B), while 2652 SNVs and 247 haplotype variants had a lower level (3 or 4). The SNV variants were then overlapped with the variants in the Sardinian population and their prevalence was evaluated based on their allele frequencies in both a dataset of 1577 unrelated Sardinians and in the gnomAD dataset (non-Finnish Europeans), which established extensive interpopulation differences. In this analysis, we considered only variants in non-cytochrome genes, yielding results as follows.
Among the variants associated with at least one of the highest levels of evidence, we identified, in our analysis, 13 variants for which the absolute difference in allele frequency between Sardinians and Europeans was at least 5%, of which 5 had a MAF in Sardinia that was at least 10% higher than the Europeans (Table 4 for a summary, Table 5 for clinical annotations). These 13 high-evidence pharmacogenetic variants were involvedaccording to 26 PharmGKB clinical annotations-in the response to the following 9 different drugs: the anticoagulants Acenocumarol, Penprocoumon, and Warfarin; the cholesterollowering agent Pravastatin; the anticancer agents Capecitabine, Fluorouracil; and the immunosuppressants Etanercept, Methotrexate, and Rituximab. Many of the variant-drug pairs (2 out of 3) were relevant for dosage and efficacy of the analyzed drug, while the remaining were relevant for toxicity (Table 5). Interestingly, about half of these differentiated pharmacogenetic variants with high priority level were relevant for drugs directed to pediatric populations. In more detail, 4 variants in 4 genes (VKORC1, CYP4F2, DPYD and MTHFR) were marked as important for pediatric prescription of 4 drugs (Phenprocoumon, Warfarin, Fluorouracil and Methotrexate). For 4 of these variants, we had at least one alert with Evidence level 1A or 1B (rs2108622, rs9923231 and rs9934438, whose minor allele was about 10% more frequent in Sardinians than in Europeans, whereas the rs1801265 minor allele was about 6% less frequent in Sardinians) ( Table 5). It should be noted that SNPs rs9934438 and rs9923231 are in high LD (r2 = 1) in both Sardinians and Europeans.  Among the other pharmacogenetic variants with lower levels of evidence (3 and 4), we identified 169 variants for which the minor allele frequency was at least 10% higher in Sardinians compared to Europeans, and a further 72 variants with an opposite trend (Table 6 for a summary, Supplementary Table S6for clinical annotations). The 169 low-evidence variants were involved in the response to 201 different drugs (a total of 405 variant-drug pairs). Most variant-drug pairs were relevant for drug efficacy (39% of annotations) and for toxicity (37% of annotations) (Supplementary Table S6). Overall, 41 of the variant-drug pairs (about 10%) were relevant for drugs directed to the pediatric population. For another 110 variant-drug pairs with lower levels of evidence (5 of which were relevant for the pediatric population), 72 unique variants, for which the minor allele frequency was at least 10% lower in Sardinians compared to Europeans, were involved in the response to 91 different drugs.

Discussion
We estimated the potential impact of the large-scale introduction of pharmacogenetic practices in the Sardinian population by evaluating the prevalence of clinically relevant pharmacogenetic variants in a core set of 1577 unrelated sequenced individuals, representative of the entire population (Sardinia has 1.5 M residents on the island and a similar number of individuals of Sardinian descent spread across the world). To this end, we used PGxPOP [21], a PGx matching engine that is based on PharmCAT and uses its PGx allele definitions, to characterize PGx allele and phenotype frequencies. Using this analysis, it was possible to estimate the theoretical number of Sardinian individuals exposed to adverse reactions to a range of drugs. In more detail, the frequencies of two phenotypes ("Decreased warfarin dose" and "Possibly decreased warfarin dose") involving warfarin, a widely used anticoagulant drug, were among the most interesting findings from this analysis. The two atypical phenotypes are determined by diplotypes of the VKORC1 gene [24] and affected a total of 1192 individuals in our cohort (i.e., about 3 of 4 individuals). Overall, common genetic variants in this gene, but also in CYP2C9, CYP4F2, and the CYP2C cluster (e.g., rs12777823), plus known nongenetic factors, account for 50% of warfarin dose variability [24].
Other phenotypes potentially affecting a large proportion of the population were the "Intermediate Metabolizer" and "Poor Metabolizer" phenotypes, which are determined by cytochrome CYP2C9 diplotypes and affected a total of 650 individuals in our cohort (about 41% of the cohort analyzed). These phenotypes have important effects on the ADME-Tox of Nonsteroidal Anti-Inflammatory Drugs, such as celecoxib, flurbiprofen, lornoxicam, and ibuprofen. According to CPIC guidelines [25], the diplotypes involved may result in a higher-than-normal risk of adverse events, especially in individuals with other factors affecting clearance of these drugs, such as hepatic impairment or advanced age. The same guidelines suggest a reduced dosage of these drugs and monitoring of adverse effects. The same cautions can be extended to other drugs, such as meloxicam, piroxicam and tenoxicam.
An important finding concerned two atypical phenotypes related to the SLCO1B1 gene ("Decreased Function", N = 324 individuals, and "Poor function", N = 51), which globally affect almost 1 in 4 individuals, and are important for the metabolism of important drugs, such as Atorvastatin (second among the top thirty active drugs both for consumption and expenditure in Italy) [26], and Fluvastatin, Lovastatin, Pitavastatin, Pravastatin, Rosuvastatin and Simvastatin. According to CPIC guidelines [27], these phenotypes can impact the starting dose and suggest an adjustment of doses based on disease-specific guidelines. According to suggestions in the same guidelines, prescribers should be aware of possible increased risk for myopathy.
We could then hypothesize that an important impact on the frequency of adverse effects could be caused by the high diffusion of atypical phenotypes ("Intermediate metabolizer", N = 450, "Poor metabolizer", N = 32, "Rapid metabolizer", N = 343 and "Ultrarapid Metabolizer", N = 37) attributable to diplotypes of the gene CYP2C19, involved in the metabolism of some of the most widely used antidepressants in Italy, including escitalopram and sertralin. According to the guidelines [28], among the problems caused by an incorrect dosage are increased risk for adverse cardiac and cerebrovascular events.
Special attention should be paid to the 103 individuals (approximately 6.5%) who are at high risk of severe toxicity due to antineoplastic drugs, such as azathioprine, mercaptopurine, and thioguanine, because of atypical phenotypes determined by diplotypes of the TPMT gene.
In a second phase of analysis, we aimed to identify the variants of pharmacogenetic interest that were more differentiated in Sardinia than in the general European population (taking as reference the genetic data of gnomAD version 2.1). In this analysis, we distinguished highly relevant PGx variants (levels of evidence 1A, 1B, 2A and 2B) from those of lower relevance (levels 3 and 4).
The strongest difference in terms of allele frequency was seen for the rs396991 variant located in the FCGR3A gene and which could be relevant for patients treated with Rituximab, according to a 2B level of evidence documented by PharmGKB [29]. In fact, the C allele of the rs396991 variant had a frequency 1.5 times higher in Sardinia (AF = 0.528) than in the rest of Europe (AF = 0.344) This difference may significantly affect the efficacy of Rituximab, used in the treatment of certain types of cancer and autoimmune disorders, including Rheumatoid Arthritis and Neuromyelitis Optica. Indeed, patients with a CC genotype may have an increased response to the drug compared to patients with AA and AC genotypes.
Another variant of special interest was rs8050894, for which the frequency of the G allele was 1.34 times higher in Sardinia (AF = 0.523) than in the general European population (AF = 0.389). This variant has a role, supported on a type 1B level of evidence, in influencing warfarin dosage. According to the guidelines, patients with the GG genotype may require a lower dose of warfarin as compared to patients with the CC genotype. The variant is part of a haplotype of variants in the VKORC1 gene, all of which are associated with warfarin dosing. Among them, the one with the strongest level of evidence was rs9923231, whose T allele was 1.309 times more frequent in Sardinians (AF = 0.509 versus 0.389). This last variant was also relevant for the pediatric population and had relevance not only for the dosage of Warfarin, Acenocoumarol and Phenprocoumon, but also for the resultant efficacy and toxicity of these drugs. Of note, the genotypes of VKORC1-1639G > A (rs9923231) are mentioned in the FDA Label of Warfarin.
Warfarin inhibits VKORC1 to prevent regeneration of a reduced form of vitamin K necessary for clotting factor activation [30]. The common variants, noted in our analysis, are located in the 5 UTR and introns of the VKORC1 gene and are associated with reduced gene expression and related effects on warfarin dosage. Warfarin and Acenocoumarol are common oral anticoagulant prescribed for the treatment and prevention of thromboembolic events for which genetic variants in several genes (CALU, calumenin; CYP, cytochrome P450 family members; GGCX, gamma-glutamyl carboxylase; NQO1, NAD(P)H quinone dehydrogenase 1; VKORC1, vitamin K epoxide reductase) have been associated with the need for carefully calibrated dosage to prevent bleeding episodes.
The influence of VKORC1 polymorphisms on vitamin K antagonist dose requirements provides a remarkable example of pharmacogenomic diversity worldwide. This is documented by the International Warfarin Pharmacogenomic Consortium (IWPC) datasets, comprising 5700-6200 patients recruited from four continents, and ascribed to three 'racial' groups, namely Asians, Blacks (mainly African Americans) and Whites [31].
Furthermore, considering the 19 HLA alleles associated with adverse events to the therapy with the highest level of evidence, mention should be made of HLA-B*58:01, which has been shown to have a strong effect on the development of severe cutaneous adverse reactions (SCARs), including Stevens-Johnson syndrome and toxic epidermal necrolysis after treatment with allopurinol, the common treatment for hyperuricemia and gout. However, the frequency of HLA-B*58:01 significantly differs between different ethnic groups. The frequency of HLA-B*5801 is the highest in Han-Chinese (20%), Korean (12%), and Thai (13%), but is much less frequent in Japanese (0.1%) The same allele, however, is also much more frequent in Sardinians (11%) than in other European populations (France 1.5%) [32]. We believe that this evidence is of particular relevance, given that Sardinia has the highest percentage of reports of adverse events [26] following allopurinol administration of the total number of adverse reports registered in Italy (1.9% compared to an average of 0.41% in the other Italian regions).
Among the most differentiated variants with lower levels of evidence, two independent variants rs3815087 (allele A) and rs3131003 (allele A) located in PSORS1C1 region, were of particular interest. Both variants were highly frequent in Sardinians compared to European populations (delta frequency > 29%), have been associated with epidermal necrolysis and Stevens-Johnson syndrome [33,34] after allopurinol therapy (evidence levels 3 and 4, respectively) and show coincident, strong association with psoriasis (p = 1.2 × 10 −294 , OR = 2.93; p = 1.4 × 10 −105 , OR = 1.64) [https://genetics.opentargets.org (accessed on 22 March 2022); rs3815087 and rs3131003 variants respectively]. They were very common in Sardinia (AR 50% and 74.6%), and, thus, screening for these variants before therapy could be important. It is, thus, not surprising that the variant rs2233945, localized in the same PSORS1C1 gene, modulates the response to etanercept, a TNF inhibitor used for psoriasis and other autoimmune disorders, including rheumatoid arthritis. Allele A in that locus has been associated with increased etanercept efficacy in comparison to allele C: and at the same time allele A has been associated with protection from psoriasis. This variant has been in linkage disequilibrium with one canonically described for allopurinol adverse events, rs9263726, the variant tag for HLA-B*58:01.

Dataset
In this study we focused on a subset of 1577 unrelated sequenced Sardinian samples from the Sardinian sequence-based reference panel. These were mainly the unrelated parents of a larger sample set of 3514 individuals (mainly trios) sequenced at low coverage (average coverage 4.2×), which also included their children. The sample set included 2090 individuals belonging to the SardiNIA cohort study in the subregion of Ogliastra and 1424 individuals deriving from a case-control study on autoimmunity collected across the Island [16,19]. All participants signed informed consent to study protocols approved by the Sardinian Regional Ethics Committee (protocol no. 2171/CE).

Genotyping and Imputation
All the genetic analyses were performed using a genetic map based on 6602 samples genotyped with 4 Illumina arrays (OmniExpress, ImmunoChip, Cardio-MetaboChip and ExomeChip), as previously described [19]. Imputation was performed on a genomewide scale using a Sardinian sequence-based reference panel of 3514 individuals and the software Minimac3 on pre-phased genotypes. After imputation, only markers with imputation quality (RSQR) > 0.3 for estimated minor allele frequency (MAF) ≥ 1% or > 0.6 for MAF < 1% were retained for further analyses, yielding~22 million variants (20,143,392 SNPs and 1,688,858 indels).

Variant Calling and Relatedness
Variant calling and phase assignment procedures leveraged the family structure of the extended sample of 3514 individuals. Allele frequencies were estimated from the subset of 1577 unrelated individuals and compared to the corresponding frequencies available in gnomAD project v2 (European non-Finnish subset) [35].
Relatedness was estimated by computing the genome-wide proportion of pairwise IBD (π) on a random set of 1 million SNPs with an MAF > 0.05 in 1000 Genomes populations (Phase 3 v5) [https://www.internationalgenome.org/home (accessed on 1 May 2019)]. Starting from the dataset of 3514 individuals, for each pair of individuals with π > 0.05, we preferentially removed the offspring if in a trio; otherwise, we removed the individual with the larger summed value of π across all other relationships with π > 0.05. Once all removals were completed, we had a total of 1577 samples for the allele frequency analyses.  [22] and reports exact matches to the allele definitions based on the provided phased genotype data.

Pharmacogenetics Resources
Gene phenotypes were determined to have a non-typical response if any CPIC guidance recommended an alternate dosage or drug for that phenotype, as defined by McInnes et al. [21] (Supplementary Table S1).
Currently, the CPIC guidelines (https://cpicpgx.org/genes-drugs/, accessed on 14 July 2022) indicate a set of 64 recommendations for the 14 genes profiled by PGxPOP, for a total of 53 drugs. In this work, we only considered gene/drug pairs whose CPIC evidence levels are class A (for which there is a sufficient level of evidence to suggest a recommendation) (Supplementary Table S2). Although the CPIC recommendations refer to a larger number of genes, as detailed in [https://cpicpgx.org/genes-drugs/ (accessed on 16 August 2022)], this work relies only on the 14 genes considered by PGxPOP at the time the analysis described here was carried out.
Frequencies in Sardinia of star allele haplotypes and diplotypes, and estimated prevalence of non-typical response phenotypes were compared with the corresponding values in UKBiobank populations reported by McInnes et al. [21] (Supplementary Tables S3-S5), where the burden of nontypical response phenotypes for each individual is estimated by counting the number of diplotypes with predicted nontypical response phenotypes across 14 genes with phenotypes. Gene phenotypes were determined to have a nontypical response if any CPIC guidance recommended an alternate dosage or drug for that phenotype.
Structural variants (SVs) were not called for any of the considered genes.

Statistical Analysis
All statistical analyses were performed using R Studio 2022.07.1 Build 554 (R Foundation for Statistical Computing, Vienna, Austria). Significant differences in phenotype frequencies between our cohort and UKBiobank European population were analyzed by Fisher's exact test. The False Discovery Rate (FDR) proposed by Benjamin and Hochberg [36] was used to correct the multiple analyses. Results were considered statistically significant when the p-value was <0.05.

Conclusions
In this work, we have completed what is, to date, the most extensive characterization of pharmacogenetic variability in an Italian population, specifically that of Sardinia. The analysis of the prevalence of PGx risk variants presented here may stimulate initiatives to implement large-scale pharmacogenetic strategies in Italy.
The impact of pharmacogenetic variation on health is, thus, patently obvious; in a future analysis it may be useful to consider that the impact falls disproportionally with age and it may be stronger on the elderly. The effects of aging result first from the relatively high usage of pharmaceutical drugs in the elderly (https://www.kff.org/health-reform/ issue-brief/data-note-prescription-drugs-and-older-adults/ (accessed on 22 March 2022); https://hpi.georgetown.edu/rxdrugs/ (accessed on 22 March 2022)) and second from the relatively greater sensitivity to drugs in the elderly, so that the dosages, which are determined on younger adults, are often excessive for them, increasing the risk of side effects.
Variants in pharmacogenes can then further exacerbate the problems. Concerning germ-line mutations, we have dealt here with the effect of single variants on single drugs, but we know little about the effect of genetics on combinations of multiple drugs. This is again especially important for the elderly, who, compared to the middle-aged individuals who normally participate in clinical trials, are more frequently simultaneously prescribed multiple drugs (often in a 'therapeutic cascade', in which the side effects of one drug are treated with another drug). Correspondingly, additive or synergistic incidence of side effects and intensification of genetic variant effects can be expected. Furthermore, with age each individual accumulates new somatic mutations, and the liver-in which the majority of the genes relevant to ADME-Tox are expressed-is one of the tissues most exposed to environmental mutagens. It, therefore, tends to accumulate somatic mutations that can potentially further alter the function of pharmacogenes [37,38], although the extent to which this occurs has not been quantified. Further analyses should define more precisely the relative load of genetic risk as a function of age and the measures, like lower doses and substitution of drugs by others with different genetic risk profiles, that may mitigate it.
There are some limitations to this study that could be met in the future. First, only by using high coverage sequencing data on the exome or genome will it be possible to define with certainty the existence and prevalence of rarer variants with important effects specific to the Sardinian population. Second, we did not assess structural variants, but for example, CYP2D6 is well-known to have structural variants including copy number variability and gene rearrangements between CYP2D7-CYP2D6 known as hybrid tandems. Third, our analysis was limited to 14 genes, that could be analyzed with PGxPOP at the time this work was prepared; this limitation could be overcome by future analysis, that uses new available information on drug-gene pairs (unfortunately not implemented in PGxPOP). Two other limitations currently preclude estimation of the potential economic cost of nonstratification of patients based on genetic characteristics. The absence of personal data on drug prescriptions (which would allow us to understand how many people, and which individuals, are really at risk of adverse reactions to the drugs they use) and global data on consumption in Sardinia (which would allow us to make pharmacoeconomic estimates).
Nevertheless, the findings demonstrate the value of characterizing allele frequencies in diverse populations and highlights the need for more PGx research on understudied populations, an important step in the corresponding refined implementation of modern personalized medicine. Funding: The Authors acknowledge funding from University of Sassari and Fondazione di Sardegna (grants "fondo di Ateneo per la ricerca 2019", "fondo di Ateneo per la ricerca 2020", and "Bando competitivo Fondazione di Sardegna-2016 per progetti di ricerca con revisione tra pari") and from NIH Contract n. 75N95021C00012 "Genetic and epidemiological factors for age-related traits and diseases in the Sardinian population (SardiNIA5)".

Institutional Review Board Statement:
The study was conducted in accordance with the Declaration of Helsinki, and approved by the Sardinian Regional Ethics Committee (protocol no. 2171/CE).

Informed Consent Statement:
Informed consent was obtained from all subjects involved in the study. All participants signed informed consent to study protocols approved by the Sardinian Regional Ethics Committee (protocol no. 2171/CE).