An epistasis between dopaminergic and oxytocinergic systems confers risk of post-traumatic stress disorder in a traumatized Chinese cohort

Post-traumatic stress disorder (PTSD) is a psychiatric syndrome that occurs after trauma exposure. Neurotransmitters such as dopamine and oxytocin have been reported to be involved in neuropathology of PTSD. Previous studies indicated that the dopamine–oxytocin interaction may contribute to behavioral disorders. Thus, exploring the epistasis (gene–gene interaction) between oxytocinergic and dopaminergic systems might be useful to reveal the genetic basis of PTSD. In this study, we analyzed two functional single nucleotide polymorphisms (SNPs), rs2268498 for oxytocinergic gene OXTR and rs1801028 for dopaminergic gene DRD2 based on putative oxytocin receptor–dopamine receptor D2 (OTR–DR2) heterocomplex in a Chinese cohort exposed to the 2008 Wenchuan earthquake (156 PTSD cases and 978 controls). Statistical analyses did not find any single variant or gene–environment interaction (SNP × earthquake-related trauma exposure) associated with provisional PTSD diagnosis or symptoms. An OXTR–DRD2 interaction (rs2268498 × rs1801028) was identified to confer risk of provisional PTSD diagnosis (OR = 9.18, 95% CI = 3.07–27.46 and P = 7.37e-05) and further subset analysis indicated that rs2268498 genotypes controlled the association directions of rs1801028 and rs1801028 genotypes also controlled the association directions of rs2268498. Rs2268498 × rs1801028 is also associated with PTSD symptoms (P = 0.043). Our study uncovered a genetic and putative function-based contribution of dopaminergic–oxytocinergic system interaction to PTSD.

genes that are involved in the same pathway or interact with each other is a feasible strategy for gene-gene interaction analysis 9 , 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 cells 12 . DA plays important roles in neuropsychological processes including motivation, executive functions, reward and stress [13][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 PTSD [16][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 PTSD [21][22][23] . Particularly, the oxytocin receptor gene (OXTR) polymorphism may predict PTSD through interaction with attachment style 24 . It has been reported that OXT-DA interactions affect social behavior like motivation 25 . Furthermore, previous studies indicated that abnormal OXT-DA interactions may contribute to behavioral disorders such as autism, sexual dysfunction, addiction and depression 26 . 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 heteromer 27,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 few 24,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 SEEK 30 . Thus, we hypothesized that the genegene 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 emotionality 31 that is one of the factors of PTSD 1,32,33 ; this SNP also affects the gene expression of OXTR 34 . Rs1801028 (risk allele C) is a non-synonymous SNP of DRD2, which has been frequently studied and known to be associated with schizophrenia 35,36 which shares genetic risk with PTSD 8 . Rs1801028 has been investigated in our study before while no significant result was obtained 37 ; 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 OXTR-DRD2 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 geneenvironment 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).
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 www.nature.com/scientificreports www.nature.com/scientificreports/ 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.
Association between OXTR-DRD2 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).    www.nature.com/scientificreports www.nature.com/scientificreports/ OXTR-DRD2 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 OXTR-DRD2 (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 GWAS 8 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 rats 28 . 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 rats 38 ; DNA methylation of OXTR in human blood is with dynamic changes after acute psychosocial stress 39 . For DRD2, down-regulated DRD2 gene expression in striatum is correlated to stress-induced depression-like behaviors in rats 40 ; alterations of human dopaminergic genes including DRD2 are associated with low resilience to stress 41 . 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 SIFT 42 and PolyPhen-2 43 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 expression 34 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 dopaminergicoxytocinergic interaction may be involved in behavior/mental disorders 26 , 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 network 44 . 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 brain 45 , and our results were based on blood samples, all these inferences need to be carefully considered.  www.nature.com/scientificreports www.nature.com/scientificreports/ 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 OXTR-DRD2 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 generalization 46 .
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 disorders 47,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 before 32 . 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 before 51 . Also, we measured depression symptoms with the Center for Epidemiological Studies-Depression (CES-D) Scale 52 . We then collected peripheral blood samples and extracted DNA for genotyping by using a standard phenol-chloroform protocol 53 . 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 SNPscan TM 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 methods 54,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.09 55 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.