Introduction

Post-traumatic stress disorder (PTSD) is a psychiatric disorder that can occur after trauma exposure. Diagnosis of PTSD requires at least six symptoms distributed across four symptom clusters and persisting for at least one month after trauma exposure1. With the lifetime prevalence of 5.6% among trauma-exposed population2, the development of PTSD is affected by both genetic factors and environment triggers3,4,5. Recent studies also indicated that inherited and acquired risk of genetic and epigenetic have significant contributions to response to trauma exposure6. To explore genetic mechanisms of PTSD, many genetic studies have been performed7, including the most recent genome-wide association study (GWAS) of Psychiatric Genomics Consortium (PGC)8. Most of the PTSD genetic studies have utilized single nucleotide polymorphisms (SNPs) as markers, and focused on main effects or gene–environment interaction effects of single gene variant; however, gene–gene interaction (epistasis) has seldom been explored. Here gene–gene interaction is defined as that the effect of one locus genotype on a phenotype depends on the genotype at other loci9. Since interactions between loci contribute to the biological and biochemical pathways that underpin disease, detecting their interactions will be useful to elucidate the genetic mechanisms underlying the complex diseases such as PTSD9,10,11. Testing genetic variants of correlated genes that are involved in the same pathway or interact with each other is a feasible strategy for gene–gene interaction analysis9, with reduced searching space and increased statistical power.

The neurotransmitters, such as dopamine (DA), oxytocin (OXT) and serotonin, are endogenous chemicals that transmit signal across a synapse between two nerve cells12. DA plays important roles in neuropsychological processes including motivation, executive functions, reward and stress13,14,15. Abnormal dopaminergic signaling have been revealed to be associated with several mental disorders, including attention deficit hyperactivity disorder (ADHD), drug and alcohol dependence, schizophrenia and PTSD16,17,18,19. OXT is needed for childbirth and is involved in regulations of both social behaviors (such as social bonding, sexual activity and motivation) and nonsocial behaviors (such as brain development, feeding and memory)20. OXT has been reported to be associated with schizophrenia, autism spectrum disorders and PTSD21,22,23. Particularly, the oxytocin receptor gene (OXTR) polymorphism may predict PTSD through interaction with attachment style24. It has been reported that OXT–DA interactions affect social behavior like motivation25. Furthermore, previous studies indicated that abnormal OXT–DA interactions may contribute to behavioral disorders such as autism, sexual dysfunction, addiction and depression26. The mechanism of their interactions has been investigated for several years and the results indicated that dopamine receptor D2 (DR2) encoded by dopaminergic gene DRD2 and oxytocin receptor (OTR) encoded by oxytocinergic gene OXTR may directly work together, by forming OTR–DR2 heteromer27,28.

There have been some candidate gene association studies on PTSD with DRD2 and OXTR respectively; the results of DRD2 were not consistent (only 11 out of total 20 studies showed significant results)7 and the results of OXTR were quite a few24,29. Furthermore, DRD2 and OXTR are co-expressed in multiple tissues of normal human brain according to the results of gene co-expression search engine SEEK30. Thus, we hypothesized that the gene–gene interaction between the two genes at genetic variant level may be involved in the PTSD development. There are two SNPs, rs2268498 in OXTR and rs1801028 (Ser311Cys) in DRD2, are quite meaningful and important. Rs2268498 (risk allele C) is localized in the promoter flanking region of the OXTR gene, which was reported to correlated to negative emotionality31 that is one of the factors of PTSD1,32,33; this SNP also affects the gene expression of OXTR34. Rs1801028 (risk allele C) is a non-synonymous SNP of DRD2, which has been frequently studied and known to be associated with schizophrenia35,36 which shares genetic risk with PTSD8. Rs1801028 has been investigated in our study before while no significant result was obtained37; rs2268498 has not been investigated in PTSD yet. In this study, to investigate our hypothesis, we analyzed the OXTR SNP rs2268498 and the DRD2 SNP rs1801028 in a traumatized Chinese cohort and examined the association between the gene–gene interaction and provisional PTSD diagnosis.

Results

Association between OXTRDRD2 and provisional PTSD diagnosis

The demography of all subjects is present in Supplementary Table 1. The genotyping call rate for both SNPs is 100%. The minor allele frequencies (MAFs) of rs2268498 are 0.2724 in cases and 0.316 in controls respectively; the MAFs of rs1801028 are 0.03846 in cases and 0.03681 in controls respectively (Supplementary Table 2). The Hardy-Weinberg equilibrium (HWE) test results are also listed in Supplementary Table 2. For provisional PTSD diagnosis inferred from PCL-5, there were no significant main effects for single gene-based analysis (P > 0.05), either without or with gene–environment interaction (SNP × earthquake-related trauma exposure) included in the logistic regression model (Supplementary Table 2). Gene–gene interaction analysis showed that rs2268498 × rs1801028 interaction was associated with provisional PTSD diagnosis (odds ratio (OR) = 9.18, 95% confidence interval (CI) = 5.25–16.05, P = 7.37e-05. Table 1 shows more detailed information). The association signals came from both females and males (P = 0.000494 and 0.0200 respectively); the effect size (OR) of the interaction was much larger in males than that in females (OR = 62.51 and 7.60 respectively) (See Supplementary Tables 3 and 4 for more details).

Table 1 Summary of logistic regression model of rs2268498 × rs1801028 for PTSD diagnosis.

Rs2268498 and rs1801028 controlled association for each other

Figure 1 shows genotype frequency variability between case and control subjects, respectively for the subsets obtained based on genotype and the full set. In Fig. 1(A–C) are for rs1801028 and (D–F) are for rs2268498. When we divided samples by rs2268498 genotypes CC/CT and TT, rs1801028 minor allele C was associated with increased risk of provisional PTSD diagnosis (OR = 3.02, 95% CI = 1.32–6.94, P = 0.009163) in rs2268498 CC/CT set and with decreased risk of provisional PTSD diagnosis (OR = 0.10, 95% CI = 0.01–0.78, P = 0.02863) in rs2268498 TT set (Table 2). Rs2268498 was associated with decreased risk of provisional PTSD diagnosis (OR = 0.71, 95% CI = 0.52–0.97, P = 0.03287) in the rs1801028 GG set (Table 3). The sample size of genotype subsets in females or in males was too small to present, although the association trends were consistent with those of overall samples (details not shown). All results of subset analysis were confirmed by permutation analysis (Tables 2 and 3). Thus, rs2268498 and rs1801028 controlled association with provisional PTSD diagnosis for each other.

Figure 1
figure 1

Genotype frequency variability between case and control. (A–C) For rs1801028 in the rs2268498 CC/CT set, the rs2268498 TT set and the full set, respectively. (D–F) For rs2268498 in the rs1801028 CC/CG set, the rs1801028 GG set and the full set, respectively.

Table 2 Logistic regression analysis of rs1801028 (with minor allele C) in the rs2268498 CC/CT set, the rs2268498 TT set and the full set, respectively.
Table 3 Logistic regression analysis of rs2268498 (with minor allele C) in the rs1801028 CC/CG set, the rs1801028 GG set and the full set, respectively.

Association between OXTRDRD2 and PTSD symptoms

For PTSD symptoms, no significant main effects or G × E effects was found (Supplementary Table 5). The G × G OXTR–DRD2 is associated with PTSD severity (rs2268498 × rs1801028: beta = 3.43, 95% CI = 0.11–6.75, P = 0.043). See Table 4 for detailed information. The association signal is not strong and both females and males contribute to the association signal (See Supplementary Tables 6 and 7 for more details). Subset analysis did not obtain any significant result (details are not shown).

Table 4 Summary of linear regression model of rs2268498 × rs1801028 for PTSD symptoms.

OXTRDRD2 in different age groups

Our results showed that for provisional PTSD diagnosis, rs2268498 × rs1801028 is more significant in young adults defined as subjects with age <50 years old (OR = 9.24 and P = 0.00470) than in old adults defined as subjects with age > = 50 years old (OR = 6.76 and P = 0.0226); for PTSD symptoms, rs2268498 × rs1801028 is significant in young adults (beta = 5.11 and P = 0.0129) but not in old adults (beta = −1.08 and P = 0.7183). Supplementary Table 8 showed the details.

Discussion

Through epistasis (gene–gene interaction) analysis, our study identified the OXTRDRD2 (rs2268498 × rs1801028) interaction to be associated with provisional PTSD diagnosis and PTSD symptoms. Interestingly, rs2268498 and rs1801028 controlled each other’s association directions for provisional PTSD diagnosis. The interaction signals came from both females and males. The two SNPs have been investigated in PGC-PTSD GWAS8 and we listed the results in Supplementary Table 9. All the association test P values (with different gender and ethnic group combinations) in PGC-PTSD GWAS are far from genome-wide significance (5e-8) or suggestive significance (e-5), indicating that the two SNPs are not associated with PTSD directly and our G × G analysis provides feasibly new information. It needs to mention that there is about six times as many controls as cases in our study (978:156), which just reflects the provisional PTSD prevalence in our epidemiologic samples and is feasible.

Our results indicated that the G × G (rs2268498 × rs1801028) effects are different between young adults (age < 50 years old) and old adults (age > = 50 years old) (Supplementary Table 8). In our data, age is significantly correlated to both provisional PTSD diagnosis (OR = 1.07 and P = 4.72e-10, by logistic regression analysis) and PTSD symptoms (beta = 0.32 and P = 3.24e-23, by linear regression analysis). Furthermore, a two-sided fisher’s exact test showed a distribution difference of provisional PTSD diagnosis cases between two age groups (P = 3.24e-8) and a two-sided t-test showed a distribution difference of PTSD symptoms between the two age groups (t = −5.46 and P = 5.97e-8). Thus, the effect of age on PTSD may contribute to the differences of G × G effects between different age groups.

The OXT–DR2 interaction was supported by evidence from animal experiments of male rats28. Our study showed that the statistical interaction exists in both women and men. In details, the association of females is more significant than that of males, which may due to the larger samples size of females than males; the interaction effect size of males is much larger than that of females, showing that males may be more sensitive to the interaction. In summary, the interaction looks like generally applies for both females and males. OXTR and DRD2 have been separately reported to be correlated with stress response through gene expression or DNA methylation which may regulate gene expression. For OXTR, it was reported that OXTR expression significantly increased in the amygdala in response to chronic social instability stress for adult female rats38; DNA methylation of OXTR in human blood is with dynamic changes after acute psychosocial stress39. For DRD2, down-regulated DRD2 gene expression in striatum is correlated to stress-induced depression-like behaviors in rats40; alterations of human dopaminergic genes including DRD2 are associated with low resilience to stress41. Our statistical genetics finding further indicated that there is a gene-gene interaction between the two genes in the context of stress response.

Rs1801028 (G > C) is corresponding to the amino acid residue polymorphism S311C of human DRD2 protein and both SIFT42 and PolyPhen-243 predicted that S311C may be with deleterious impact to DR2 (Supplementary Table 10). Together with the function of the OXTR SNP rs2268498 (influencing OXTR expression)34, it is possible to infer the SNP–SNP interaction to impact biological function level. It has been found that at protein level, DR2 and OTR interacts with each other by forming heterocomplex. Our study found that a functional genetic variant interaction of rs2268498 × rs1801028 might confer risk of provisional PTSD diagnosis. The further subset analysis showed that when rs2268498 C allele corresponding to high OXTR gene expression34 presented, rs1801028 C allele (unstable DRD2 protein) increased PTSD risk. When rs1801028 GG genotype (stable DRD2 protein) presented, rs2268498 C allele (high OXTR gene expression) decreased PTSD risk. This may be explained that high level OTR and unstable DRD2 are difficult to form feasible heterocomplex, which may lead to dynamic disequilibrium of OTR–DRD2 heteromer, influencing PTSD development. This provides evidence that dopaminergic–oxytocinergic interaction may be involved in behavior/mental disorders26, which deserves further investigation in other mental disorders. Furthermore, besides OTR–DRD2 heteromer, there are many G protein-coupled receptor (GPCR) heteroreceptor complexes involved in brain network44. Our results suggest a possibility that these heteroreceptor complexes may be replicated in PTSD, as well as in other mental disorders. Still, since genetic control of gene expression can differ between blood and brain45, and our results were based on blood samples, all these inferences need to be carefully considered.

There are several limitations in our study. First, PTSD diagnosis was based on PCL-5, a self-report measurement, other than clinician-administered instruments, which will need to be used in future studies. Second, only one trauma type- exposure to the 2008 earthquake was examined. Our results would need to be confirmed in other trauma-related PTSD cohorts. Third, the sample size of our study is modest in comparison to other genetic studies. Moreover, the male sample size is relatively smaller than that of the female’s (due to the female-male gender ratio of our epidemiologic samples). Hopefully, a much larger sample with balanced sexes could be used to replicate our findings in future. Also, there is no information about pre-earthquake trauma exposure or some conditions to which PTSD symptoms may be attributed (such as substance use, bereavement, etc.) and PTSD symptoms which may preclude the conclusion that PTSD is developed after the earthquake; still, causality can be inferred given that genes inherently serve as the first temporal variable. In future study, we will take the information into account.

Our study firstly identified the OXTRDRD2 interaction at genetic variant level to be significantly associated with a higher PTSD risk based on provisional PTSD diagnosis. Since the genetic variants correspond to polymorphisms in mRNA or protein level, the genetic interaction provides a possibility of existence of functional interaction, suggesting that abnormal interaction of dopaminergic and oxytocinergic systems may be a vulnerable factor for PTSD development. Of course, we should note these interpretations were based on provisional PTSD diagnosis. In summary, our gene–gene interaction results may provide new viewpoints in terms of the biological mechanisms and genetic architectures of PTSD, which will be helpful for early diagnosis and prevention.

Materials and Methods

Study sample

The study samples we used are epidemiologic samples which are feasible for generalization46. In detail, all the subjects were recruited from Hanwang Town, Mianzhu City, China. In the 2008 Wenchuan earthquake, the town was almost destroyed. Five and a half years after the earthquake, we conducted the survey in a large rebuilt community of the town. Subjects were predominantly Chinese adults (more than sixteen years old) who survived the earthquake. All of them were without a major psychiatry history or mental retardation such as schizophrenia and organic mental disorders47,48. We assessed PTSD symptoms by the PTSD Checklist for DSM-5 (PCL-5)49,50, a reliable and repeatable 20-item self-reported scale to capture the PTSD symptoms of DSM-5. The items of PCL-5 particularly focus on the earthquake-related trauma and the Chinese version of PCL-5 has been used for earthquake-exposure Chinese before32. We instructed participants to complete the PCL-5 referring to the “Wenchuan Earthquake”. The likely DSM-5 PTSD diagnosis was then inferred based on PCL-5, including at least one intrusion symptom, one avoidance symptom, two negative alterations in cognitions and mood symptoms, and two arousal symptoms endorsed as 2 or greater. Earthquake-related trauma exposure was considered as environment factor and was measured by a ten-question scale. The sum score was calculated by adding 10 items (0: not experienced or 1: experienced for each item) as the level of trauma exposure. The details of the ten questions have been described before51. Also, we measured depression symptoms with the Center for Epidemiological Studies-Depression (CES-D) Scale52. We then collected peripheral blood samples and extracted DNA for genotyping by using a standard phenol-chloroform protocol53. For each sample, at least 120 ng DNA was used for genotyping experiment. A total of 1134 subjects (156 PTSD cases and 978 controls) with complete clinical information, and genotyping data were included in our study.

SNP genotyping and quality control

We investigated rs2268498 in oxytocinergic gene OXTR and rs1801028 (Ser311Cys) in dopaminergic gene DRD2. We conducted genotyping by using a custom-by-design 2 × 48-Plex SNPscanTM Kit (Genesky Biotechnologies Inc., Shanghai, China), which is based on a multiplex fluorescence PCR and a double ligation. The PCR primers were showed in Supplementary Table 11. To assess genotyping accuracy, we randomly chose 10% samples for replication by the same platform. The concordance was 100%, showing that the genotyping is accurate. A Hardy-Weinberg equilibrium (HWE) test was performed by exact tests with efficient computational methods54,55. For quality control (QC), we excluded subjects with missing call rate > 0.2, and removed SNPs with call rate < 0.95 or HWE test P < 0.05 for controls or MAF (minor allele frequency) < 0.01 in controls. No sample or SNP was removed for the quality control.

Statistical analyses

All of our statistical analyses were performed by using PLINK v1.0955 and R3.4.4 (https://cran.r-project.org/). We first analyzed the association for provisional PTSD diagnosis. For a single variant-based analysis, we employed an additive model of SNP minor allele count per subject in a logistic regression with gender, age, depression symptoms (CES-D scale score, the possible confounding factor) and trauma exposure as covariates. The gene–environment interaction (SNP × trauma exposure) was further included. The results indicated that the gene–environment interaction was not significant, thus it was not included in following analyses. For gene–gene interaction analysis, we utilized an additive model of minor allele count per subject in a logistic regression with predictors SNP1, SNP2, SNP1 × SNP2, gender, age, depression symptoms and trauma exposure. To explore the possible different association properties between females and males, we also carried on the gene–gene interaction analysis for females and males, respectively. We then investigated the association for PTSD symptoms (PCL-5 total score which is with rang 0–80 and treated as a continuous variable) by using the similar methods as above, except that we employed linear regression instead of logistic regression.

To explore the detailed gene–gene interaction property, we employed subset analysis. We divided all samples into two subsets according to rs2268498 genotype: the CC/CT set (samples with C alleles) and the TT set (samples with genotype TT). Then we performed regressions (with gender, age, depression symptoms and trauma exposure as covariates) for rs1801028 respectively in the two subsets. Correspondingly, we stratified all samples by rs1801028 genotype (CC/CG and GG), examining the association for rs2268498 respectively.

In our analysis, age is significant factor affecting provisional PTSD diagnosis. The age of our study subjects is in a wide range (16–73) while the median of age is 48 and more than 50% of the subjects are with ages in the much narrower range (42, 56) (Supplementary Table 1 and Supplementary Fig. 1). To explore possible different association property between diverse age groups, we divided out samples into young adults (age < 50 years old) and old adults (age > = 50 years old). We then performed association test to detecting G × G (rs2268498 × rs1801028) effect for the two groups respectively.

All the P values of the regressions in our analyses are two-sided. Threshold of P value is considered as 0.05. We employed a permutation test with 1,000,000 cycles for SNP-related variables to adjust for possible bias due to the small sample size.

Statements of ethical approval, accordance and informed consent

The Institutional Review Board of Institute of Psychology, Chinese Academy of Sciences approved this study. All study procedures involving human participants were in accordance with the national and institutional research committee’s ethical standards. All participants signed informed consents. Among them, two were under the age of 18 years and signed informed consents with permission of their legal guardians who also provided written informed consents.