The pituitary adenylate cyclase-activating polypeptide (PACAP) and its receptor (PAC1) are widely distributed in hypothalamic and limbic structures (1–6). This peptide and receptor help regulate responses to anxiety-provoking stimuli or stress (3, 7–10). The ADCYAP1R1 gene, which encodes the receptor PAC1, is associated with the neuroendocrine system hypothalamic–pituitary–adrenal (HPA) axis that is involved in controlling mammalian stress responses (11, 12).

There are many single-nucleotide polymorphisms (SNPs) in ADCYAP1R1. One of them, rs2267735, was reported to play a role in anxiety symptoms, fear-related cognitions, and stress-related responses (13–15). For example, rs2267735 could regulate the normal stress response by affecting the bind of estrogen receptors alpha and estrogen response element (16).

Posttraumatic stress disorder (PTSD) is a trauma- and stressor-related disorder (17) that can occur after exposure to traumatic events (18). Many studies have shown that PTSD is associated with ADCYAP1R1 and ADCYAP1R1–environment interactions [e.g., (19–22)]. In a Chinese population, a significant association between ADCYAP1R1 and the PTSD emotional numbing cluster was found, a finding that extended previous results (23).

Moreover, PTSD is correlated with other genes that are associated with the HPA axis or regulation of its activity. The FKBP5 gene regulates the HPA axis by encoding FK506 binding protein 5 (FKBP5). FKBP5 can bind to glucocorticoid receptors (24). The FKBP5 gene was detected to have a significant correlation with PTSD in a meta-analysis study (25, 26) and could reflect the risk of co-morbidity of PTSD and depression in mild trauma exposure (27). FKBP5-environment interactions predicted the risk of PTSD in adults (28–30). The CRHR1 gene and the CRHR2 gene regulate the HPA axis together by binding their encoding products, corticotropin-releasing hormone (CRH) receptors CRF1 and CRF2, to CRH. The CRHR1 gene has been suggested to be associated with PTSD (31–36). The CRHR2 gene has been shown to significantly affect PTSD (32, 37, 38) as well as through a gene–environment interaction (39).

Numerous neurogenetic studies of the ADCYAP1R1 gene have been performed. Some of them gave out a significant main effect on PTSD. For example, Ressler and her colleagues have found that ADCYAP1R1 could contribute to the diagnosis of PTSD and the severity of PTSD symptoms in females who suffered a major trauma (21). Later, Lind et al. (20) used a meta-analysis to show that the C allele at rs2267735 of ADCYAP1R1 significantly increased the risk of PTSD in women. A review supported as well that the epigenetic regulation of ADCYAP1R1 might predict PTSD risk (40). However, not all studies performed the same, and they even led to no effects. In order to replicate the finding of Ressler et al. (21), two independent samples were newly selected, and their results failed to repeat such (41). In addition, a recent investigation did not show that ADCYAP1R1 was strongly correlated with PTSD (42). The possible reasons for why these study results have not been consistent may be due to the types of trauma exposure (43, 44) and the sample characteristics (23, 41). However, the underlying factors have been most likely not singular. Therefore, we need more studies to enrich the evidence of the main or gene-environment effects of ADCYAP1R1 on PTSD.

Furthermore, the four genes ADCYAP1R1, FKBP5, CRHR1, and CRHR2 are all associated with PTSD and have a close relationship with each other by regulating the HPA axis together. Since genes functionally interact with each other, a phenomenon that is usually called epistasis (i.e., gene-gene interaction) in genetics (45), the gene-gene interactions between ADCYAP1R1 and other HPA axis genes need to be investigated. However, to our knowledge, a neurogenetic study that investigates the gene–gene interactions between ADCYAP1R1 and other HPA axis genes in PTSD has not been conducted. Thus, the correlation between the interactions of ADCYAP1R1 with other HPA axis genes (including FKBP5, CRHR1, andCRHR2) and PTSD is worth further exploring.

To further explore the role of ADCYAP1R1 in PTSD development and the physiological underpinnings of the stress response, this study examined the main, G × E (rs267735 × trauma exposure), and G × G (with other HPA axis gene polymorphisms) effects of rs2267735 on PTSD severity in a cohort of 2008 Chinese Wenchuan earthquake survivors.

In the present study, we selected 1,132 survivors of the Wenchuan earthquake in Hanwang Town, Sichuan Province, China. These participants have been described previously (38, 46), and we provide the detailed information for the cohort here again for convenience.

The participants were adults (older than 18 years old) without any mental illness or intellectual disability history. Regarding self-reported gender results, 68.19% of the participants were women, and the remaining participants were men. Moreover, most participants reported that their ethnicities were Han (99.65%), and the remaining were reported to be Qiang. Furthermore, the marital status of 13.07% of the participants was unmarried. In addition, 32.50% of the participants had an educational level involving high school or above.

The measurement occurred in a large rebuilt community of Hanwang Town from November 6–8 in 2013. The investigators were all trained clinical psychologists, psychotherapists, and psychology graduate students. Our study was approved by the institutional review board of the Institute of Psychology, Chinese Academy of Sciences, and was completely in compliance with national legislation and the Declaration of Helsinki. All the participants signed informed consent forms.

Earthquake-related trauma exposure was assessed by a self-reported questionnaire (46). The participants needed to answer 10 yes (1) or no (0) questions about whether they experienced (a) being trapped under a rubble, (b) being injured, (c) being disabled due to injuries, (d) participating in rescue efforts, (e) witnessing the death of someone, (f) being exposed to mutilated bodies, (g) the traumatic death of a family member, (h) the traumatic injury of a family member, (i) the traumatic death of a friend or neighbor, and (j) losing livelihood due to a disaster. The level of trauma exposure was defined as the sum score of the 10 items.

PTSD symptoms were assessed by a self-reported questionnaire, the PTSD Checklist of DSM-5 [PCL-5; (47, 48)]. The 20 items included on the PCL-5 used a five-point Likert scale (from 0 = never to 4 = most of the time). The Chinese version has shown good validity and reliability (49). The participants needed to complete the PCL-5 Chinese version according to their PTSD symptom occurrence frequency and severity in the past month of the earthquake. The final score indicating PTSD symptom severity was the sum of all the item scores.

Depressive symptoms were measured by a self-report questionnaire, the Center for Epidemiological Studies-Depression (CES-D) scale (50). The 20 items included on the CES-D used a four-point Likert scale (from 0 = rare or none of the time/ <1 day to 3 = most or all of the time/5–7 days). The Chinese version has been validated and widely used in Chinese populations (51). The participants were required to answer the Chinese version of the questionnaire on the basis of their personal experiences in the past week. The final score indicating depressive symptoms was reflected by the sum of all the item scores.

Four HPA axis genes with nine relevant SNPs were selected for study: ADCYAP1R1 (rs2267735), FKBP5 (rs9296158, rs3800373, rs1360780, and rs9470080), CRHR1 (rs4458044 and rs242924), and CRHR2 (rs8192496 and rs2267715). These SNPs were genotyped in a previous study (38); thus, we directly used their data.

We used PLINK (52) to perform the Hardy–Weinberg equilibrium (HWE) test (53). In addition, we also calculated the SNP call rate and minor allele frequency (MAF).

PLINK 1.07 and R 3.4.4 (https://www.r-project.org/) were used to conduct all the statistical analyses in our study. A pairwise comparison of demographic variables and the SNP rs2267735 was tested by analysis of variance (ANOVA) in R. A linkage disequilibrium (LD) analysis has been done previously (38), and we directly referenced those results.

The associations between rs2267735 (coded as minor allele counts per subject) and PTSD severity were detected by a linear regression model and adjusted for age, sex (1: women and 0: men), education level (1: high school or above and 0: other), marital status (1: single and 0: married), the environmental factor trauma exposure, and depressive symptoms. The main, gene–environment interaction (G × E), and gene–gene interaction (G × G) effects were considered. When studying the G × E effects, we further added a linear model analysis with all gene–covariate and environment–covariate interactions (54). The all-interactions model was used to control for potential confounders and compared to the single-interaction model. A subset analysis for women and men was performed to detect potential sex differences. All P-values of linear regressions in our study are two-sided.

For adjustment of possible bias and correction for multiple comparisons, the SNP-related variables were subjected to a permutation test with 100,000 cycles. The results with significant P and permutation P-values are reported. Moreover, we calculated the effect size (semipartial correlation) as well.

The summary of trauma exposure, PTSD severity (total symptom) scores, and depressive symptoms (CES-D scale total score) are shown in Table 1. The mean value of depressive symptoms in all samples is 37.00, and the standard error of the mean was 8.63. According to the DSM-5, we inferred probable PTSD in our participants. As our likely diagnoses showed, the probable PTSD prevalence was 13.78% (14.89% for women and 11.38% for men). All of the above-mentioned data were exactly the same as those of our previous study (38) since they were based on the same cohort, and we provided the data here for convenience.

The HWE test indicated that the genotype frequencies of rs2267735 agreed with the Hardy–Weinberger equilibrium (P = 0.812). The call rate was 100%, and the MAF (for allele G) was 0.485. The pairwise comparisons of the demographic variables and SNP rs2267735 are shown in Supplementary Table 1. There was no significant result between the SNP and depression.

Based on previous LD analysis results (38), during the following analyses, we used rs2267735 × rs1360780 to index the G × G effect of ADCYAP1R1–FKBP5 and rs2267735 × rs242924 to index the G × G effect of ADCYAP1R1–CRHR1.

As shown, the minor allele G of rs2267735 of ADCYAP1R1 was not significantly associated with PTSD severity. The G × E ADCYAP1R1–environment (rs2267735 × trauma exposure) effect on PTSD severity was also not significant. The trauma exposure was measured by a 10-item scale in the present study [see Supplementary Tables 2, 3 (single-interaction model) and Supplementary Table 4 (all-interaction model) for details].

The G × G ADCYAP1R1–FKBP5 (rs2267735 × rs1360780) effect on PTSD severity was significant in all samples (P = 0.049). The G × G ADCYAP1R1–CRHR1 (rs2267735 × rs242924) effect on PTSD severity was significant in men (P = 0.023; see Tables 2, 3 for details). Figure 1 shows the G × G (epistasis) effects represented by the distributions of PTSD severity across genotype combinations. We also provide all the other results of the G × G analysis (Supplementary Table 5).

In this study, two gene–gene interactions (G × G) were identified. An ADCYAP1R1–FKBP5 effect was associated with PTSD severity in all subjects, and an ADCYAP1R1–CRHR1 effect was found in men.

An FKBP5-environment interaction and ADCYAP1R1 have been found to be associated with PTSD (19–22, 28–30), but to our knowledge, few studies have shown a G × G effect of FKBP5 and ADCYAP1R1. Therefore, this is the first study showing that the ADCYAP1R1 polymorphism rs2267735 could affect PTSD development through a novel ADCYAP1R1–FKBP5 interaction. Moreover, this finding was observed under the conditions of earthquake trauma type and in a Chinese ethnic group.

In addition, CRHR1 has been reported in prospective studies to moderate childhood maltreatment effects on PTSD symptoms (28, 30, 31). The gene–gene–environment (CRHR1 × 5–HTTLPR × childhood maltreatment) interaction could predict adult depressive symptoms among black people of lower socioeconomic status (55). Here the identified G × G effect between ADCYAP1R1 and CRHR1 on PTSD in men was interesting and further revealed the influence of gene–gene interactions on PTSD.

Earlier findings of PTSD candidate gene studies have shown that rs2267735 in ADCYAP1R1 could increase amygdala reactivity and reduce functional connectivity between the amygdala and hippocampus (21, 56). In addition, four SNPs in FKBP5 (rs9296158, rs3800373, rs1360780, and rs9470080), which are related to regulating the stress response system, were shown to increase amygdala reactivity to threat stimuli and the severity of PTSD symptoms in adulthood (28, 57). Thus, the ADCYAP1R1–FKBP5 interaction may influence PTSD development by affecting amygdala reactivity together.

Both ADCYAP1R1 and CRHR1 are involved in stress response by modulating CRF function and the release of cortisol through the adrenal cortex (58). Therefore, the ADCYAP1R1–CRHR1 interaction suggests that their gene expression may influence PTSD by regulating CRF together.

The genotypes of rs2267735 were differentially distributed in females and males (Supplementary Table 2). To explore sex differences in the two G × G (ADCYAP1R1–FKBP5 and ADCYAP1R1–CRHR1) interactions, we used Cochran's Q statistic test by PLINK. To investigate heterogeneity between results from women and men (i.e., the sex–G × G interactions), we also calculated the heterogeneity index I² (range, 0–100). The ADCYAP1R1-FKBP5 interactions were not different between women and men (Cochran's Q statistic P < 0.668 and I² = 0), which indicated that the effects were sex non-specific. The ADCYAP1R1-CRHR1 interactions were different between sexes (Cochran's Q statistic P < 0.027 and I² = 79.4), which showed that sex might moderate the effects. Many studies have shown that sex could regulate the relationship between ADCYAP1R1 and PTSD. Women with the C allele at rs2267735 of ADCYAP1R1 had a higher PTSD risk than those with the G allele (20, 21). Moreover, a genome-wide association study indicated that CRHR1 is a genetic factor related to PTSD reexperiencing symptom (33). Our findings suggested that ADCYAP1R1 might also affect PTSD in men in a specific manner by interacting with CRHR1.

Our study had some limitations. One was the scales used in the current study. The PTSD was measured by a self-reported scale, and a clinical assessment could be used in the future. The trauma exposure was measured by a 10-item scale. Although it has shown great reliability and validity (38, 46, 59–61), other questionnaires with more items could be tried to be applied further. Depression was regarded as a covariate, so it was mentioned less in the study. The second was that the sample size was relatively small, and the number of women was approximately two times that of men. Large sample sizes and sex-balanced samples should be considered. The single ethnic group and trauma exposure type also need to be broadened. The third was about the definition of the G × E effect. Because this was a candidate gene study, so the “G” was only considered as specific genotype instead of genetics, and the “E” was only concentrated on earthquake, which was the most related environment factor to the PTSD symptom measured in the present study. Other possible related factors, such as rearing environment, social support, socioeconomic status, parenting styles, and education, could be taken into consideration in the future. The fourth was the sex difference of ADCYAP1R1–CRHR1 interactions which might be due to the sex hormone level which was different between men and women. This point could be detected through physiological experiments and indications.