Peripheral Blood Mononuclear Cell Oxytocin and Vasopressin Receptor Expression Positively Correlates with Social and Behavioral Function in Children with Autism

The peptide hormone oxytocin is an established regulator of social function in mammals, and dysregulated oxytocin signaling is implicated in autism spectrum disorder (ASD). Several clinical trials examining the effects of intranasal oxytocin for improving social and behavioral function in ASD have had mixed or inclusive outcomes. The heterogeneity in clinical trials outcomes may reflect large inter-individual expression variations of the oxytocin and/or vasopressin receptor genes OXTR and AVPR1A, respectively. To explore this hypothesis we examined the expression of both genes in peripheral blood mononuclear cells (PBMC) from ASD children, their non-ASD siblings, and age-matched neurotypical children aged 3 to 16 years of age as well as datamined published ASD datasets. Both genes were found to have large inter-individual variations. Higher OXTR and AVPR1A expression was associated with lower Aberrant Behavior Checklist (ABC) scores. OXTR expression was associated with less severe behavior and higher adaptive behavior on additional standardized measures. Combining the sum expression levels OXTR, AVPR1A, and IGF1 yielded the strongest correlation with ABC scores. We propose that future clinical trials in ASD children with oxytocin, oxytocin mimetics and additional tentative therapeutics should assess the prognostic value of their PBMC mRNA expression of OXTR, AVPR1A, and IGF1.

The peptide hormone oxytocin is strongly implicated in social behavior and is conserved in mammals. Oxytocin has been extensively studied in the context of social function in animals and humans, in particular in children and adults with autism spectrum disorder (ASD) [1][2][3] . Intranasal oxytocin was shown to improve social function and empathy in healthy individuals [4][5][6] , and is thus a candidate therapeutic for children with ASD.
It has been suggested that genetic variations of the OXTR alleles, promoter methylation, and/or expression levels in brain areas implicated in social function may contribute to the variable response observed in clinical trials with intranasal oxytocin in ASD 7,8 . In addition, differences between and within animal species in sociality may in part be explained by substantial differences in OXTR expression levels within the social brain network, such as the amygdala, the anterior cingulate cortex, the prefrontal cortex and the nucleus accumbens 9 . In primates, OXTR is also expressed in the superior colliculus, the pulvinar, and the primary visual cortex, where oxytocin functions as a modulator of visual processing and allocation of attention 10 . Thus, variations in OXTR expression seem critical for the diversity of social behaviors across and within mammal and in particular primate species, including humans. Indeed, inter-individual variations in OXTR expression levels were shown to be associated with resilience to the effects of neonatal isolation on adult social attachment in female prairie voles 11 .

Results
The age, sex, and key ASD scores of study participants are presented in Table 1. As both oxytocin and vasopressin are implicated in ASD 1-3,24,25 , we searched for expression differences of the genes which code for their receptors, OXTR and AVPR1A, in PBMCs collected from the study participants. No statistically significant differences were found for the PBMC mRNA expression levels of OXTR or AVPR1A between ASD children, their neurotypical siblings, and/or neurotypical age-matched children (Fig. 1a,b). Data mining of GEO datasets derived from ASD and control PBMC (GSE111176, GSE25507) or whole blood (GSE18123, GSE89594) confirmed the lack of differences in OXTR or AVPR1A expression levels in PBMC or whole blood from ASD and neurotypical controls ( Fig. 1c-j).
Next, we searched for correlations between PBMC expression levels of OXTR and AVPR1A and the available social and behavioral scores of our study participants (see Methods; scores were recorded at the same time as blood sample collection for PBMC separation). We found that higher OXTR expression levels correlated with i) better development as measured by the VABS Adaptive behavioral composite, ii) less severe social impairment as measured by the SRS total t-score and iii) less severe behavior problems as measured by the lower total raw ABC scores and lower total CBCL t-scores ( Fig. 2a-d). The PBMC expression of AVPR1A correlated with the participants total raw ABC scores, but not with their VABS, SRS, or CBCL scores ( Fig. 2e; Supplementary Fig. 1e-g). Notably, no correlations with standardized scores were observed for the PBMC expression of CD38 ( Supplementary Fig. 1a-d), a gene that codes for an enzyme that synthesizes and hydrolyzes cyclic adenosine 5′-diphosphate-ribose and taking part in oxytocin secretion 31 . In addition, we also observed a correlation of the participants total raw ABC scores and their expression levels of IGF1, coding for insulin-like growth factor 1 (Fig. 2g).
Since combining the expression levels of several genes may afford more robust diagnostic or prognostic power [32][33][34][35][36] and as the ABC scores of the ASD children correlated with the OXTR, AVPR1A, and IGF1 expression levels, we determined whether combining the expression of these genes affords a better correlation with ABC scores as compared to correlation with individual gene expression. As shown in Fig. 2h the combined expression levels of all three genes yielded improved Pearson correlations (larger r values and smaller p values) compared to each gene separately or with two-gene combinations (Fig. 2f). In order to verify the robustness of the correlation between the summed OXTR, AVPR1A, and IGF1 expression levels with ABC scores, we randomly selected 30 individuals from our dataset. For each of 500 random selections we determined the Pearson's correlation coefficient and p-value of the correlation between ABC score and the summed gene expression. We found that 410 (82%) of these random subsets yielded a significant correlation (p-value < 0.05) between ABC score and summed gene expression (mean Spearman's correlation coefficient: 0.47). This is compared to 25 (5%) significant subsets expected by chance. This result demonstrates that the observed correlation is robust and is not underlined by specific outlier individuals.

Discussion
Dysfunctional signaling by the peptide hormones oxytocin and vasopressin were suggested as contributing to ASD. While the role of these hormones has been extensively studied in several mouse models of ASD 37,38 , and oxytocin receptor knockout mice display behavioral deficits resembling autism-related behaviors 39 , such mouse models typically reflect a single mutation or deletion, and do not reflect the real-world state and wide spectrum of ASD individuals. Our study aimed to explore correlations between development and behavior in children with ASD and PBMC expression levels of the oxytocin and vasopressin receptor genes OXTR or AVPR1A, respectively.
Our findings indicate that the severity of developmental and behavioral deficits in children with ASD aged 3 to 16 years is associated with their PBMC mRNA expression levels of OXTR (ABC, VABS, SRS, and CBCL scores), as well as with AVPR1A and IGF1 (ABC score). The PBMC expression of all three genes therefore seem to be associated with the severity of aberrant behavior in children with ASD. www.nature.com/scientificreports www.nature.com/scientificreports/ Our literature search (PubMed search of May 2019) identified 11 completed and published clinical trials (listed in ClinicalTrials.gov) of intranasal oxytocin prescribed for at least 5 days and including at least 10 ASD participants who received oxytocin 12,13,[15][16][17][18][19][20][21][22] (Table 2). Only six of the 11 published studies reported a significant favorable effect of intranasal oxytocin. Notably, of the five studies with a placebo-controlled crossover design, three showed a favorable effect ( Table 2). The low number of such clinical trials does not allow reaching a conclusion whether the study design (placebo-controlled crossover, placebo-controlled, or open label) affected the findings. Of note, three out of four trials with children showed favorable outcomes (Table 2) suggesting that intranasal oxytocin might be more likely to be efficacious for children. While findings from these 11 trials are difficult to interpret owing to their variable age groups, dosages, treatment durations and study designs, they clearly indicate that only some ASD individuals may benefit from oxytocin.
Of the 11 trials listed in Table 2, a single trial 22 examined the OXTR expression levels, and reported lack of correlation of its expression with the effect of intranasal oxytocin on behavioral scores of children diagnosed with ASD. However, this trial measured OXTR expression in whole blood and compared it with the participants Social Responsiveness Scale (SRS) Total Raw Score (while we studied PBMC mRNA and report correlations with ABC scores).
Among the studies listed in Table 2 two studies examined the effects of participants OXTR SNPs on the efficacy of intranasal oxytocin 12,13 ; both studies applied a placebo-controlled crossover design in adolescents and adults with ASD. Watanabe et al. 12 reported that intranasal oxytocin has a smaller effect for participants carrying OXTR rs53576 or rs2254298; while Kosaka et al. 13 reported lower efficacy for those carrying OXTR rs6791619. Of note, all three of these OXTR SNPs are intronic. We searched for the effects of these OXTR SNPs on its expression in datasets derived from studies with postmortem brain tissues (134 control individuals) deposited in the UK Brain Expression Consortium (UKBEC; http://braineac.org/) by the MRC Sudden Death Brain and Tissue Bank 40 . Our datamining showed that postmortem brain tissues from individuals carrying the minor (A) allele rs53576 allele had significantly higher OXTR expression levels in several brain regions, including frontal cortex, compared with tissues from GG homozygous (Fig. 3). In the general population, rs53576 was associated with general sociality 41 and empathy 42,43 . Taken together, our current findings on poorer ABC scores in individuals with lower PBMC OXTR expression, along with the findings of the oxytocin clinical trial by Watanabe et al. 12,23 on better efficacy (improved behavioral scores) in the OXTR rs53576 minor allele carriers, suggest that higher OXTR expression levels may correlate with improved efficacy of intranasal oxytocin in ASD individuals.
Unlike OXTR SNPs, OXTR mRNA levels seem to be more informative of ASD severity: besides DNA methylation and additional epigenetic modifiers of ASD severity 44,45 , they may inform on effects by other non-heritable effectors, such as the gut microbiome, suggested to affect ASD 46,47 . Based on our findings, we suggest that higher peripheral OXTR expression as detected in PBMCs may reflect differential brain expression and allow improved response to prescribed oxytocin and/or future oxytocin derivatives/agonists. The same considerations may apply also for the prognostic value of higher expression levels of AVPR1A and IGF1 (Fig. 2e-h); while, being supported by the clinical trial findings by Kosaka et al. 13 , the predictive value seems most robust for the PBMC expression of OXTR. Future clinical trials with intra-nasal oxytocin, or future oxytocin derivatives such as oxytocin 5-9 or synthetic agonists 48,49 , should explore the predictive value of the participants' PBMC mRNA of these genes. Findings from such clinical trials with ASD children or adults may allow the future stratification of ASD individuals for the specific therapeutic most likely to benefit them. Once validated by clinical trials, this may allow the future co-marketing of a drug/test combination, along the successful drug/test combinations applied for biologic therapeutics in oncology, such as the Herceptin/Her2 drug/test combination 50 . Considering the increasing rates and high societal cost of ASD, such research efforts seem timely and warranted.

participants, Materials and Methods
participants. The protocol is registered in clinicaltrials.gov as NCT02000284 and was approved by the Institutional Review Board at the University of Arkansas for Medical Sciences (Little Rock, AR). Parents of participants provided written informed consent. All methods were performed in accordance with the relevant guidelines and regulations. Children underwent a fasting blood draw in the morning. Control individuals did not have any neurological disorders or developmental delays. Thus, this study reports gene expression and social and behavioral measurements on a total of 63 children, including 43 children with ASD, 11 non-ASD siblings, and 9 unrelated neurotypical control children. This was a subset of the larger cohort of our study of mitochondrial www.nature.com/scientificreports www.nature.com/scientificreports/ function in children with ASD due to limited availability of PBMC for gene expression analysis. Our recent publications outline the methodology for rating participant characteristics of this participant cohort 51-53 .

Behavioral measurements.
As mentioned in our previous study 53 , our research staff was trained by a multispecialty team consisting of two licensed psychologists and a speech therapist prior to performing assessments. During the study a research psychologist supervised research staff and provided feedback and retraining if necessary. Observer rated measures included the VABS Survey Interview Form. Parents completed the ABC and the SRS. The VABS is a reliable and valid measure of the ability to perform age-appropriate everyday skills, including communication, daily living, social and motor skills, through a 20-30 minute structured interview with a caretaker 54 . Of note, functional abilities of children with autism are commonly measured with the VABS 55 as it was in this study. Although some have used IQ to distinguish high and low functioning autism, recent studies indicate that it is a poor predictor 56 . The ABC is a 58-item questionnaire 53 that measures disruptive behaviors and has convergent and divergent validity 57 . The SRS is a 65-item questionnaire that measures the severity of social skill deficits across five domains 58 which has been shown to have good correspondence to the gold-standard instrument 59 . The CBCL includes demographic information and ratings of positive behaviors, academic functioning, social competence, and behavior problems commonly applied for classification of behavior disorders 60 .
Blood collection and processing. Samples of 4 ml of venous blood were collected into an ethylenediaminetetraacetic acid-Vacutainer tube, chilled on ice and centrifuged at 1500 g for 10 minutes at 4 °C to separate plasma. Plasma was removed and stored at −80 °C for later analysis. Plasma was replaced with room temperature wash buffer containing Ca +2 /Mg +2 -free PBS with 0.1% BSA and 2 mM ethylenediaminetetraacetic acid. Diluted  (a,b). Corroborating our findings, no differences between ASD and control children were found in Gene Expression Omnibus (GEO) datasets from PBMC (c-f) or whole blood (g-j) samples analyzed by datamining for their OXTR and AVPR1A expression. GSE codes and cohort sizes are shown. Note the large distribution of expression levels for OXTR and AVPR1A in both our PBMC samples (a,b) and the GEO datasets (c-j). Between-group differences were analyzed by one-way ANOVA test (a,b) and by Student's t-test (c-j). www.nature.com/scientificreports www.nature.com/scientificreports/ blood was then layered on top of Histopaque-1077 (Sigma Aldrich, St. Louis, MO, USA) and centrifuged at 400 g for 30 minutes at room temperature. PBMCs were collected, washed twice with wash buffer and counted using a hemocytometer.
RnA isolation. Total RNA was extracted using RNeasy mini kit (Qiagen, Hilden, Germany) from 5 million PBMCs by following manufacturer provided protocol. RNA samples were shipped in 100% ethanol at room   www.nature.com/scientificreports www.nature.com/scientificreports/ temperature to Tel-Aviv University. Upon arrival to Tel-Aviv, 10% 3 M Sodium Acetate and 1ul GlycoBlue ™ (Thermo Fisher Scientific, MA, USA) were added and kept overnight at −80 °C. Samples were centrifuged at high speed (14,000 g, 4 °C) for 30 minutes. Pellets were dried and dissolved in 20ul DNAse/RNAse-free water. Geo datamining. The NCBI Gene Expression Omnibus (GEO) was queried for expression data sets derived from ASD blood or PBMC samples containing cohorts of at least 30 ASD children or adults and a similarly sized age and sex matched control cohort. Data sets from human cell lines, postmortem tissues, and mouse ASD models were excluded. The identified data sets were queried for the expression levels of OXTR and AVPR1A using the GEO2R tool on the NCBI server (https://www.ncbi.nlm.nih.gov/geo/geo2r/). OXTR rs53576 associated brain expression levels were downloaded from the Brain eQTL Almanac website (http://www.braineac.org/). The www.nature.com/scientificreports www.nature.com/scientificreports/ expression levels corresponding to each OXTR rs53576 genotype (AA, AG, and GG) were plotted for six brain regions (Fig. 3). Of note, the above website did not yield information for the OXTR expression associated with rs2254298 or rs6791619 (the two other SNPs mentioned in our Discussion).
Statistical analyses. qPCR data analysis was conducted using the GraphPad Prism v.6 (San Diego, CA, USA). Normality of data distribution was evaluated using the Shapiro-Wilk test; continuous variables between two groups were analyzed by Student's t-test; outliers were detected by Grubbs test. Continuous variables between three groups were analyzed by one-way ANOVA test. P-values ≤ 0.05 were considered as significant. Data for correlation analysis were log transformed and Spearman correlation test was performed. P-values ≤ 0.01 were considered significant. P values were adjusted for multiple testing using Benjamini-Hochberg FDR correction 61 . R statistical software was used to test the robustness of correlations between the sum expression levels of OXTR, AVPR1A, and IGF1 and the ABC scores (Fig. 2h). A total of 500 random subsets (each containing data from 30 individuals) were independently sampled from our data. For each random sample, the Pearson's correlation coefficient and p-value of the correlation between ABC score and the summed gene expression was determined. The frequency of correlations with a p-value smaller than 0.05 was then calculated.