Author response: A sex-specific evolutionary interaction between ADCY9 and CETP (2021)

Article Figures and data Abstract Introduction Results Discussion Materials and methods Appendix 1 Data availability References Decision letter Author response Article and author information Metrics Abstract Pharmacogenomic studies have revealed associations between rs1967309 in the adenylyl cyclase type 9 (ADCY9) gene and clinical responses to the cholesteryl ester transfer protein (CETP) modulator dalcetrapib, however, the mechanism behind this interaction is still unknown. Here, we characterized selective signals at the locus associated with the pharmacogenomic response in human populations and we show that rs1967309 region exhibits signatures of positive selection in several human populations. Furthermore, we identified a variant in CETP, rs158477, which is in long-range linkage disequilibrium with rs1967309 in the Peruvian population. The signal is mainly seen in males, a sex-specific result that is replicated in the LIMAA cohort of over 3400 Peruvians. Analyses of RNA-seq data further suggest an epistatic interaction on CETP expression levels between the two SNPs in multiple tissues, which also differs between males and females. We also detected interaction effects of the two SNPs with sex on cardiovascular phenotypes in the UK Biobank, in line with the sex-specific genotype associations found in Peruvians at these loci. We propose that ADCY9 and CETP coevolved during recent human evolution due to sex-specific selection, which points toward a biological link between dalcetrapib's pharmacogene ADCY9 and its therapeutic target CETP. Introduction Coronary artery disease (CAD) is the leading cause of mortality worldwide. It is a complex disease caused by the accumulation of cholesterol-loaded plaques that block blood flow in the coronary arteries. The cholesteryl ester transfer protein (CETP) mediates the exchange of cholesterol esters and triglycerides between high-density lipoproteins (HDL) and lower density lipoproteins (Lagrost, 1994; Shinkai, 2012). Dalcetrapib is a CETP modulator that did not reduce cardiovascular event rates in the overall dal-OUTCOMES trial of patients with recent acute coronary syndrome (Schwartz et al., 2012). However, pharmacogenomic analyses revealed that genotypes at rs1967309 in the ADCY9 gene, coding for the ninth isoform of adenylate cyclase, modulated clinical responses to dalcetrapib (Tardif et al., 2015). Individuals who carried the AA genotype at rs1967309 in ADCY9 had less cardiovascular events, reduced atherosclerosis progression, and enhanced cholesterol efflux from macrophages when treated with dalcetrapib compared to placebo (Tardif et al., 2015; Tardif et al., 2016). In contrast, those with the GG genotype had the opposite effects from dalcetrapib. Furthermore, a protective effect against the formation of atherosclerotic lesions was seen only in the absence of both Adcy9 and CETP in mice (Rautureau et al., 2018), suggesting an interaction between the two genes. However, the underlying mechanisms linking CETP and ADCY9, located 50 Mb apart on chromosome 16, as well as the relevance of the rs1967309 non-coding genetic variant are still unclear. Identification of selection pressure on a genetic variant can help shed light on its importance. Adaptation to different environments often leads to a rise in frequency of variants, by favoring survival and/or reproduction fitness. An example is the lactase gene (LCT) (Bersaglieri et al., 2004; Enattah et al., 2007; Gamba et al., 2014; Itan et al., 2009; Poulter et al., 2003), where a positively selected intronic variant in MCM6 leads to an escape from epigenetic inactivation of LCT and facilitates lactase persistence after weaning (Labrie et al., 2016). Results of genomic studies for phenotypes such as adaptation to high-altitude hypoxia in Tibetans (Yi et al., 2010), fatty acid metabolism in Inuits (Fumagalli et al., 2015) or response to pathogens across populations (Hollenbach et al., 2001) have also been confirmed by functional studies (Li et al., 2018; Reynolds et al., 2020; Tashi et al., 2017; Meyer et al., 2018; Blais et al., 2012). Thus, population and regulatory genomics can be leveraged to unveil the effect of genetic mutations at a single non-coding locus and reveal the biological mechanisms of adaptation. When two or more loci interact during adaptation, a genomic scan will likely be underpowered to pinpoint the genetic determinants. In this study, we took a multi-step approach on the ADCY9 and CETP candidate genes to specifically study their interaction (Figure 1). We used a joint evolutionary analysis to evaluate the potential signatures of selection in these genes (Step 1), which revealed positive selection pressures acting on ADCY9. Sex-specific genetic associations between the two genes are discovered in Peruvians (Step 2), a population in which natural selection for high-altitude was previously found on genes related to cardiovascular health (Crawford et al., 2017). Furthermore, our know-down experiments and analyses of large-scale transcriptomics (Step 3) as well as available phenome-wide resources (Step 4) bring further evidence of a sex-specific epistatic interaction between ADCY9 and CETP. Figure 1 Download asset Open asset Flowchart of experimental design and main results. The four main steps of the analyses conducted in this study are reported along with the datasets used for each step and the genetic loci on which the analyses are performed. Green colored boxes represent analyses for which sex is considered. Abbreviations: KD = Knock down; OX = Overexpression. Results Signatures of selection at rs1967309 in ADCY9 in human populations The genetic variant rs1967309 is located in intron 2 of ADCY9, in a region of high linkage disequilibrium (LD), in all subpopulations in the 1000 Genomes Project (1000G), and harbors heterogeneous genotype frequencies across human populations (Figure 2a). Its intronic location makes it difficult to assess its functional relevance but exploring selective signals around intronic SNPs in human populations can shed light on their importance. In African populations (AFR), the major genotype is AA, which is the homozygous genotype for the ancestral allele, whereas in Europeans (EUR), AA is the minor genotype. The frequency of the AA genotype is slightly higher in Asia (EAS, SAS) and America (AMR) compared to that in Europe, becoming the most frequent genotype in the Peruvian population (PEL). Using the integrated haplotype score (iHS) (Voight et al., 2006) (Step 1 a, Figure 1), a statistic that enables the detection of evidence for recent strong positive selection (typically when |iHS| > 2), we observed that several SNPs in the LD block around rs1967309 exhibit selective signatures in non-African populations (|iHSSAS| = 2.66, |iHSEUR| = 2.31), whereas no signal is seen in this LD block in African populations (Figure 2b, Appendix 1—figure 1, Appendix 1). Our analyses suggest that this locus in ADCY9 has been the target of recent positive selection in several human populations, with multiple, possibly independent, selective signals detectable around rs1967309. However, recent positive selection as measured by iHS does not seem to explain the notable increase in frequency for the A allele in the PEL population (fA = 0.77), compared to the European (fA = 0.41), Asian (fA = 0.44), and other American populations (fA = 0.54 in AMR without PEL). Figure 2 Download asset Open asset Natural selection signature at rs1967309 in ADCY9. (a) Genotype frequency distribution of rs1967309 in populations from the 1000 Genomes (1000G) Project and in Native Americans (NAGD). (b) Significant iHS values (absolute values above 2) for 1000G continental populations and recombination rates from AMR-1000G population-specific genetic maps, in the ADCY9 gene. (c) PBS values in the ADCY9 gene, in CHB (outgroup, left panel), PEL (middle panel), and MXL (right panel). Horizontal lines represent the 95th percentile PBS value genome-wide for each population. Vertical dotted black lines define the LD block around rs1967309 (black circle) from 1000G population-specific genetic maps. Gene plots for ADCY9 showing location of its exons are presented in blue below each plot. Abbreviations: Altaic from Mongolia and Russia: ALT; Uralic Yukaghir from Russia: URY; Chukchi Kamchatkan from Russia: CHK; Northern American from Canada, Guatemala and Mexico: NOA; Central American from Costal Rica and Mexico: CEA; Chibchan Paezan from Argentina, Bolivia, Colombia, Costa Rica, and Mexico: CHP; Equatorial Tucanoan from Argentina, Brazil, Colombia, Gualana and Paraguay: EQT; Andean from Bolivia, Chile, Colombia and Peru: AND. For 1000G populations, abbreviations can be found here https://www.internationalgenome.org/category/population/. Figure 2—source data 1 Source file for genotype frequency distribution of rs1967309. This file contains the genotype frequency of rs1967309 for each subpopulation from 1000G and NAGD, with the number of individuals by subpopulation. https://cdn.elifesciences.org/articles/69198/elife-69198-fig2-data1-v2.txt Download elife-69198-fig2-data1-v2.txt Figure 2—source data 2 Source file for iHS plot in the ADCY9 gene. This file contains the iHS values for each position in the ADCY9 gene for each population of the 1000G dataset. https://cdn.elifesciences.org/articles/69198/elife-69198-fig2-data2-v2.txt Download elife-69198-fig2-data2-v2.txt Figure 2—source data 3 Source file for PBS plots in the ADCY9 gene. This file contains the PBS value for PEL, MXL, and CHB for each position in the ADCY9 gene. https://cdn.elifesciences.org/articles/69198/elife-69198-fig2-data3-v2.txt Download elife-69198-fig2-data3-v2.txt To test whether the difference between PEL and other AMR allele frequencies at rs1967309 is significant, we used the population branch statistic (PBS) (Step 1b, Figure 1). This statistic has been developed to locate selection signals by summarizing differentiation between populations using a three-way comparison of allele frequencies between a specific group, a closely related population, and an outgroup (Yi et al., 2010). It has been shown to increase power to detect incomplete selective sweeps on standing variation. Applying this statistic to investigate rs1967309 allele frequency in PEL, we used Mexicans (MXL) as a closely related group and a Chinese population (CHB) as the outgroup (Methods). Over the entire genome, the CHB branches are greater than PEL and MXL branches (meanCHB = 0.020, meanMXL = 0.008, meanPEL = 0.009), which reflects the expectation under genetic drift. However, the estimated PEL branch length at rs1967309 (Figure 2c), which reflects differentiation since the split from the MXL population (PBSPEL,rs1967309=0.051, empirical p-value = 0.014), surpasses the CHB branch length (PBSCHB,rs1967309=0.049, empirical p-value > 0.05), which reflects differentiation since the split between Asian and American populations, whereas no such effect is seen in MXL (PBSMXL,rs1967309=0.026, empirical p-value > 0.05), or for any other AMR populations. Furthermore, the PEL branch lengths at several SNPs in this LD block (Figure 2c) are in the top 5 % of all PEL branch lengths across the whole genome (PBSPEL,95th = 0.031), whereas these increased branch lengths are not observed outside of the LD block (Figure 2c). These results are robust to the choice of the outgroup and the closely related AMR population (Materials and methods). The increase in frequency of the A allele at rs1967309 is also seen in genotype data from Native American populations (Reich et al., 2012), with Andeans showing genotype frequencies highly similar to PEL (fA = 0.77, Figure 2a). The PEL population has a large Andean ancestry (Materials and methods, Appendix 1—figure 2a and b) and almost no African ancestry, strongly suggesting that the increase in AA genotype arose in the Andean population and not from admixture with Africans. The PEL individuals that harbor the AA genotype for rs1967309 do not exhibit a larger genome-wide Andean ancestry than non-AA individuals (p-value = 0.30, Mann-Whitney U test). Overall, these results suggest that the ancestral allele A at rs1967309, after dropping in frequency following the out-of-Africa event, has increased in frequency in the Andean population and has been preferentially retained in the Peruvian population's genetic makeup, potentially because of natural selection. Evidence for co-evolution between ADCY9 and CETP in Peru The pharmacogenetic link between ADCY9 and the CETP modulator dalcetrapib raises the question of whether there is a genetic relationship between rs1967309 in ADCY9 and CETP, both located on chromosome 16. Such a relationship can be revealed by analyzing patterns of long-range linkage disequilibrium (LRLD) (Lewontin and Kojima, 1960; Rohlfs et al., 2010), in order to detect whether specific combinations of alleles (or genotypes) at two loci are particularly overrepresented. To do so, we calculated the genotyped-based linkage disequilibrium (r2) (Step 2 a, Figure 1) between rs1967309 and each SNP in CETP with minor allele frequency (MAF) above 0.05. In the Peruvian population, there are four SNPs, (including two in perfect LD in PEL) that exhibit r2 values with rs1967309 that are in the top 1 % of r2 values (Figure 3a) computed for all 37,802 pairs of SNPs in ADCY9 and CETP genes with MAF >0.05 (Materials and methods). Despite the r2 values themselves being low (r2rs158477=0.080, r2rs158480;rs158617=0.089, r2rs12447620=0.090), these values are highly unexpected for these two genes situated 50 Mb apart (ADCY9/CETP empirical p-value < 0.006, Appendix 1—table 1) and thus correspond to a significant LRLD signal. This signal is not seen in other 1000G populations (Appendix 1—table 1). We also computed r2 between the four identified SNPs' genotypes and all ADCY9 SNPs with MAF above 0.05 (Figure 3b). The distribution of r2 values for the rs158477 CETP SNP shows a clear bell-shaped pattern around rs1967309 in ADCY9, which strongly suggests the rs1967309-rs158477 genetic association detected is not simply a statistical fluke, while the signal in the region for the other SNPs is less conclusive. The SNP rs158477 in CETP is also the only one that has a PEL branch length value higher than the 95th percentile, also higher than the CHB branch length value (PBSPEL,rs158477=0.062, Appendix 1—figure 3a), in line with the observation at rs1967309. Strikingly, this CETP SNP's genotype frequency distribution across the 1000G and Native American populations resembles that of rs1967309 in ADCY9 (Figure 3c). Figure 3 with 2 supplements see all Download asset Open asset Long-range linkage disequilibrium between rs1967309 and rs158477 in Peruvians from Lima, Peru. (a) Genotype correlation (r2) between rs1967309 and all SNPs with MAF >5% in CETP, for the PEL population. (b) Genotype correlation between the three loci identified in (a) to be in the 99th percentile and all SNPs with MAF >5% in ADCY9. The dotted line indicates the position of rs1967309. The horizontal lines in (a,b) represent the threshold for the 99th percentile of all comparisons of SNPs (MAF >5%) between ADCY9 and CETP. Figure 3—figure supplement 1 presents the same plots for Andeans and in the replication cohort (LIMAA) and Figure 3—figure supplement 2 compares the r2 values between PEL and LIMAA (c) Genotype frequency distribution of rs158477 in 1000G and Native American populations. (d) Genomic distribution of r2 values from 3,513 pairs of SNPs separated by between 50 and 60 Mb and 61 ± 10 cM away across all Peruvian chromosomes from the PEL sample, compared to the rs1967309-rs158477 r2 value (dotted gray line) (genome-wide empirical p-value = 0.01). The vertical black line shows the threshold for the 95th percentile threshold of all pairs. Gene plots showing location of exons for CETP (a) and ADCY9 (b) are presented in blue below each plot. Abbreviations: Altaic from Mongolia and Russia: ALT; Uralic Yukaghir from Russia: URY; Chukchi Kamchatkan from Russia: CHK; Northern American from Canada, Guatemala and Mexico: NOA; Central American from Costal Rica and Mexico: CEA; Chibchan Paezan from Argentina, Bolivia, Colombia, Costa Rica and Mexico: CHP; Equatorial Tucanoan from Argentina, Brazil, Colombia, Gualana and Paraguay: EQT; Andean from Bolivia, Chile, Colombia and Peru: AND. For 1000G populations, abbreviations can be found here https://www.internationalgenome.org/category/population/. Figure 3—source data 1 R2 values of all SNPs between ADCY9 and CETP genes in the PEL population from 1000G. This file contains the result from the geno-r2 command of the vcftools software for all SNPs (MAF >5%) of the PEL population between ADCY9 and CETP genes. The script to create Figure 3a and b can be found here . https://cdn.elifesciences.org/articles/69198/elife-69198-fig3-data1-v2.txt Download elife-69198-fig3-data1-v2.txt Figure 3—source data 2 Source file for genotype frequency distribution of rs158477. This file contains the genotype frequency of rs158477 for each subpopulation from 1000G and NAGD, with the number of individuals by subpopulation. https://cdn.elifesciences.org/articles/69198/elife-69198-fig3-data2-v2.txt Download elife-69198-fig3-data2-v2.txt Figure 3—source data 3 R2 values used for the null distribution in the PEL population from 1000G. 3,513 pairs of SNPs on chromosome 1–18 with a MAF between 15% and 30%, separated by between 50 and 60 Mb and 51–71 cM based on the PEL genetic map from 1000G. R2 values were obtained from the geno-r2 command of the vcftools software. https://cdn.elifesciences.org/articles/69198/elife-69198-fig3-data3-v2.txt Download elife-69198-fig3-data3-v2.txt Given that the Peruvian population is admixed (Harris et al., 2018), particular enrichment of genome segments for a specific ancestry, if present, would lead to inflated LRLD between these segments (Li and Nei, 1974; Nei and Li, 1973; Park, 2019; Slatkin, 2008), we thus performed several admixture-related analyses (Step 2b, Figure 1). No significant enrichment is seen at either locus and significant LRLD is also seen in the Andean source population (Figure 3—figure supplement 1a, b, Appendix 1). Furthermore, we see no enrichment of Andean ancestry in individuals harboring the overrepresented combination of genotypes, AA at rs1967309+ GG at rs158477, compared to other combinations (p-value = 0.18, Mann-Whitney U test). These results show that admixture patterns in PEL cannot be solely responsible for the association found between rs1967309 and rs158477. Finally, using a genome-wide null distribution which allows to capture the LRLD distribution expected under the admixture levels present in this sample (Appendix 1), we show that the r2 value between the two SNPs is higher than expected given their allele frequencies and the physical distance between them (genome-wide empirical p-value = 0.01, Figure 3d). Taken together, these findings strongly suggest that the AA/GG combination is being transmitted to the next generation more often (i.e. is likely selectively favored) which reveals a signature of co-evolution between ADCY9 and CETP at these loci. Still, such a LRLD signal can be due to a small sample size (Park, 2019). To confirm independently the association between genotypes at rs1967309 of ADCY9 and rs158477 of CETP, we used the LIMAA cohort (Asgari et al., 2020; Luo et al., 2019), a large cohort of 3509 Peruvian individuals with genotype information, to replicate our finding. The ancestry distribution, as measured by RFMix (Methods) is similar between the two cohorts (Appendix 1—figure 2a), however, the LIMAA cohort population structure shows additional subgroups compared to the 1000G PEL population sample (Appendix 1—figure 2c-e): to limit confounders, we excluded individuals coming from these subgroups (Appendix 1). In this dataset (N = 3,243), the pair of SNPs rs1967309-rs158477 is the only pairs identified in PEL who shows evidence for LRLD, with an r2 value in the top 1 % of all pairs of SNPs in ADCY9 and CETP (ADCY9/CETP empirical p-value = 0.003, Figure 3—figure supplement 1c, d, Figure 3—figure supplement 2, Appendix 1—table 1). The r2 test used above is powerful to detect allelic associations, but the net association measured will be very small if selection acts on a specific genotype combination rather than on alleles. In that scenario, and when power allows it, the genotypic association is better assessed by with a χ2 distributed test statistic (with four degrees of freedom, χ42) comparing the observed and expected genotype combination counts (Rohlfs et al., 2010). The test confirmed the association in LIMAA (χ42 = 82.0, permutation p-value < 0.001, genome-wide empirical p-value = 0.0003, Appendix 1). The association discovered between rs1967309 and rs158477 is thus generalizable to the Peruvian population and not limited to the 1000G PEL sample. Sex-specific long-range linkage disequilibrium signal Because the allele frequencies at rs1967309 were suggestively different between males and females (Figure 4—figure supplement 1), we performed PBS which that the LD block around rs1967309 is between in the Peruvians (Figure 4—figure supplement 2, Appendix 1). We further the effect of sex on the LRLD association found between rs1967309 and rs158477 and performed LRLD These analyses revealed that the correlation between rs1967309 and rs158477 is only seen in males in PEL (Figure and b, Appendix 1—figure and b, Appendix 1—table the r2 value to in males (ADCY9/CETP empirical p-value = genome-wide empirical p-value < = and in females (ADCY9/CETP empirical p-value = genome-wide empirical p-value = = In the Andean population, the association between rs1967309 and rs158477 is not significant when we by sex (Appendix 1—table 1), but we still see significant association signals with rs158477 at other SNPs in ADCY9 LD block in both (Figure 4—figure supplement Appendix 1—figure The LRLD result in PEL cannot be by of Andean ancestry between males and females (p-value = Mann-Whitney U test). A permutation analysis that the sex of that the observed difference between the is larger than we by (p-value = Appendix 1—figure Appendix 1). In the LIMAA we replicate this sex-specific result (Figure and d, Appendix 1—figure Appendix 1—table 1) where the r2 test is significant in males (ADCY9/CETP empirical p-value = 0.003, = but not in females (ADCY9/CETP empirical p-value = = The genotypic test the association between ADCY9 and CETP is present in males (χ42 = permutation p-value = 0.001, genome-wide empirical p-value = Appendix 1), an of rs158477 This is also the genotype combination the LRLD in In the test also shows a but significant effect (χ42 = permutation p-value = genome-wide empirical p-value = by an of a different genotype rs158477 AA, which however, not replicated in PEL possibly because of of power (Appendix 1). Figure with 3 supplements see all Download asset Open asset Sex-specific long-range linkage Genotype correlation between the loci identified in CETP in Figure 3a and all SNPs with MAF >5% in ADCY9 for (a,b) the PEL population and LIMAA cohort in males and in females Genotype frequencies sex are shown in Figure 4—figure supplement 1 and sex-specific PBS values in Figure 4—figure supplement The horizontal line shows the threshold for the 99th percentile of all comparisons of SNPs (MAF >5%) between ADCY9 and CETP. The vertical dotted line the position of rs1967309. represent the rs158477 SNPs and the other three SNPs identified in Figure 3a and which are in Figure 4—figure supplement 3 shows the same analysis in Andeans from Gene plots for ADCY9 showing location of its exons are presented in blue below each plot. Figure data 1 R2 values of all SNPs between ADCY9 and CETP genes in the PEL population from 1000G and LIMAA cohort in and This contains all of r2 values obtained from the geno-r2 command of the vcftools software for all SNPs (MAF >5%) of the PEL population by and and the LIMAA cohort by and between ADCY9 and CETP genes by to create those can be found Download effects on CETP gene expression LRLD between can suggest the of if are functional (Park, 2019). In order to be under selection, mutations to a or an such as gene We have shown previously (Rautureau et al., that CETP and Adcy9 interact in mice to several atherosclerotic To test whether these genes interact in we ADCY9 in (Step 3 a, Figure 1) and performed on KD biological and to evaluate the of ADCY9 expression on the We confirmed the KD was as ADCY9 expression is reduced in the KD (Figure which a in expression compared to the whole = Materials and methods). We also observed that CETP expression was increased in compared to (Figure an increase that is also significant = = This increased expression was by and also increased CETP protein (Materials and methods, Figure supplement 1a, b, Appendix 1), but its did not CETP expression (Figure supplement or CETP did not ADCY9 expression on (Figure supplement These experiments an interaction between ADCY9 and CETP at the gene expression and the that ADCY9 potentially the expression of CETP a genetic effect by rs1967309. Figure 5 with 2 supplements see all Download asset Open asset of ADCY9 on CETP (a) expression of ADCY9 or CETP genes on type and in from on biological in each were obtained from a and results in are presented in Figure supplement expression on the combination of rs1967309 and rs158477 genotypes in (b) (p-value = = = (c) in males (p-value = = = and in (d) artery in females (p-value = = = for individuals of European to reported were obtained from a interaction of a for the number of considered. Figure supplement 2 show the interaction on number of in the Figure data 1 expression of ADCY9 and CETP genes This file contains the expression of ADCY9 and CETP for the by and by in the sample from the is with the sample in the by the same number 1 to Download Figure data 2 of CETP expression by genotype. This contains all of CETP expression for of all (Materials and in the by and file by file by The number of in the is in the of the In each the values of CETP expression after for each The is the genotype of rs1967309 = AA, 1 = 2 = The is the genotype combination of the rs1967309 same coding that the and rs158477 = 1 = 2 = Download To test for potential interaction effects between rs1967309 and CETP, we used RNA-seq data from in the et al., the and et al., (Step Figure 1). When across in ADCY9 and CETP