**Figure 1.** Phenotypic (trait) profiles based on cumulative ADI-R scores for spoken language, nonverbal communication, play skills, social skills/interactions, and perseverative behaviors. In order to better compare the distribution of traits in each twin subgroup, the total number of individuals comprising each group was normalized to 100. The concordant autistic monozygotic (MZ) and dizygotic (DZ) twin subgroups are identified as cMZ and cDZ, respectively, and the discordant autistic (dA) and non-autistic (dNA) MZ and DZ twin subgroups are referred to as dA.MZ, dA.DZ, dNA.MZ, and dNA.DZ, respectively.

**Figure 2.** Boxplots showing range of trait scores of all twin groups. Student's *t*-test was used to show that there were significant differences between the trait scores for dNA.MZ and dNA.DZ twins. \*\*\* *p*-value ≤ 0.001.

**Figure 3.** *Cont*.

**Figure 3.** Comparison of average severity scores for dNA.MZ (red bars) and dNA.DZ (blue bars) twins across all ADI-R items (along horizontal axes) used to derive a cumulative trait score. \* *p*-value ≤ 0.05, \*\* *p*-value ≤ 0.01, \*\*\* *p*-value ≤ 0.001.

#### *2.2. Genotype–Phenotype Relationships Revealed by Quantitative Trait Association Analyses*

Scheme 1 illustrates the overall design and workflow for the genetic analyses conducted in this study. Based upon the severity differences in multiple traits exhibited in each of the subgroups of twins, we conducted a series of quantitative trait association analyses to discover SNPs that may associate with the phenotypic differences in trait severity. Quantitative trait association analyses were conducted using the cumulative trait scores for each individual in each twin subgroup and their respective genotyping data derived from the study by Wang et al. (2009). QTNs with unadjusted *p*-values ≤ 1 × 10−5 were selected for subsequent case-control association analyses. Results of the QTL analyses for each subgroup of twins are shown in Tables S1–S6 in Supplementary Materials. The combined sets of genes associated with the QTNs across all twin groups are shown in Table 1, with a Venn diagram showing the distribution of the overlapping genes among four of the five traits (Figure 4). Notably, there were no genes shared between those associated with perseverative behaviors and any of the other four traits. Hypergeometric distribution analyses of each set of QTN-associated genes indicated significant enrichment in genes from the SFARI database for all traits except non-verbal communication (Table 2).

**Scheme 1.** Workflow for two-phase genetic analyses using both quantitative trait and case-control association analyses of individuals stratified according to severity of traits.

**Figure 4.** Venn diagram showing SNP-associated gene overlap among language, nonverbal communication, play skills, and social skills. There was no gene overlap between the genes associated with perseverative behaviors and those associated with any of the other four traits.

Summary of genes associated with QTNs for each trait across all twin groups.

> **Table 1.**

**Table 2.** Hypergeometric distribution analyses of trait-associated genes and SFARI genes.

\* For the hypergeometric distribution analyses, the total number of genes in the population was estimated to be 20,406, and the total number of ASD-associated genes in the SFARI gene database was 1079.

#### *2.3. Case-Control Association Analyses Based on Individuals Stratified by Trait Severity for Each Subgroup of Twins*

The SNPs used for case-control association analyses were limited to the trait-associated SNPs derived from QTL analyses involving the corresponding twin groups. Results of these case-control association analyses are presented in Tables S7–S12. Examination of these results for each set of twins reveals significant SNPs (corrected for multiple testing by Bonferroni and/or Benjamini–Hochberg methods) that differentiate between cases and controls. Moreover, for all SNPs that were found to be significantly associated with both the lowest quartile (Q1) and highest quartile (Q4) cases diagnosed with ASD, the minor allele frequencies (MAF) relative to that of controls were noticeably different, resulting in different odds ratios. That is, a specific trait-associated SNP that may be considered as predisposing to ASD in the Q4 subgroup with a high odds ratio (>1.0) was found to have a lower odds ratio (often <1.0) in the Q1 subgroup, further suggesting relevance to severity of that specific trait. Although a SNP may occasionally have a higher odds ratio in the Q1 than in the Q4 case group, the differences in MAFs still suggest that these subgroups based on trait severity are genetically distinct with respect to that genetic variant. Notably, the majority of significant SNPs in the non-autistic subgroups (dNA.MZ and dNA.DZ) were found in Q4, at the highest severity end of the trait distribution.

#### *2.4. Enrichment of Known Autism Risk Genes from Case-Control Association Analyses*

The SNP-associated genes derived from the case-control association analyses described above were compared to the complete set of autism risk genes that comprise the SFARI gene database. Table 3 shows that genes mapped to the trait-associated SNPs that are significant by case-control association analyses using individuals in both Q1 and Q4 as cases for each subgroup of twins are enriched among SFARI genes in three of the six twin groups (dNA.MZ, dNA.DZ, and cDZ). Together with the results shown in Table 2, these findings suggest that the experimental strategy of combining quantitative trait association analyses with case-control association analyses using cases stratified by trait severity is capable of identifying significant autism risk genes even with relatively small numbers of samples. Furthermore, the SNPs identified can be associated with severity of specific traits, thereby associating genotype with phenotype.

**Table 3.** Hypergeometric analyses of case-control associated genes and SFARI genes.

\* For the hypergeometric distribution analyses, the total number of genes in the population was estimated to be 20,406; and the total number of ASD-associated genes in the SFARI gene database was 1079.

#### *2.5. Functional Analysis of Genes Associated with QTNs from Quantitative Trait and Case-Control Association Analyses*

Genes implicated by QTNs were analyzed using Ingenuity Pathway Analysis (IPA) software with a focus on identifying genes that impact nervous system development and function. Tables 4–8 show the neurological functions that are statistically enriched among genes implicated by QTNs associated with spoken language, nonverbal communication, play skills, social skills, and perseverative behaviors, respectively, for each subgroup of twins. Notably, many of the annotated functions are relevant to what is known about specific neuronal processes disrupted in ASD. Moreover, IPA analyses of genes associated with significant SNPs identified by case-control association analyses of both dNA.MZ and dNA.DZ subgroups in addition to the cDZ subgroup which are all enriched in SFARI genes also revealed significant over-representation of genes involved in neurological functions, behaviors, and disorders relevant to ASD, as shown in Tables 9–11, respectively.

**Table 4.** Neurological functions associated with genes implicated by language QTNs.

**Table 5.** Neurological functions associated with genes implicated by nonverbal QTNs.

**Table 6.** *Cont.*

**Table 7.** Neurological functions associated with genes implicated by social skills QTNs.

**Table 10.** Neurological functions implicated by QTN-associated genes in dNA.DZ.

**Table 11.** Neurological functions implicated by QTN-associated genes in cDZ.

**Table 11.** *Cont.*

#### *2.6. Case-Control Association Analysis Using dNA.DZ Twins as Controls for dA.DZ Co-Twins*

Given the differential diagnosis of discordant dizygotic twins with presumed genetic differences between the autistic and non-autistic co-twins, we conducted an exploratory case-control association analysis using the dNA.DZ twins (*n* = 36) as controls for the diagnosed dA.DZ co-twins (*n* = 38). Interestingly, over 140 genes were associated with a total of 390 SNPs that exhibited a nominal *p*-value ≤ 0.05 and odds ratio ≥ 3 (Table 12, complete case-control association data in Table S13). Hypergeometric distribution analyses showed that this set of genes was highly enriched for autism risk genes from the SFARI database (*q* = 1.84 × 10−7) and from an exome-sequencing study [23] (*q* = 4.64 × 10−8). Pathway analysis of this gene set using IPA revealed that netrin signaling (involving PRKG1, DCC, and PRKAG2) and axon guidance signaling (involving PRKCE, PDGFD, DCC, ROBO1, SLIT3, PRKAG2, NTNG1, PLCG2) were among the top 10 canonical pathways represented with Fisher's Exact *p*-values of 2.56 × 10−3 and 6.18 × 10−3, respectively. Several of the above-mentioned genes (specifically, PRKCE, PRKAG2, and PLCG2) are also significantly over-represented in additional pathways potentially relevant to ASD, including melatonin signaling (*p* = 1.18 × 10−2) and synaptic long-term potentiation (*p* = 4.45 × 10−2). Functional analysis showed significant over-representation of many neurological processes commonly impacted by autism, notably neuritogenesis, axonogenesis, neuronal migration, and axon guidance (Table 13).

**Table 12.** Genes associated with SNPs identified by case-control association analyses using dA.DZ twins as cases and dNA.DZ twins as controls.

\* Q indicates probability of enrichment in SFARI or exome-seq genes from hypergeometric distribution analysis. # Exome-seq genes were downloaded from study by Iossifov et al. [23].

**Table 13.** *Cont.*

#### **3. Discussion**

The goals of this study were to: (1) examine the influence of genetics on specific behaviors and traits of ASD; (2) identify SNPs (i.e., QTNs) that associate with a particular ASD trait; (3) determine if trait QTNs can discriminate between cases and controls, where the cases are individuals at the extremes of phenotypic severity across five autistic traits; and (4) to identify pathways and functions implicated by the genes associated with the QTNs. To accomplish these goals, we developed a novel approach which utilized MZ and DZ twins who are either concordant or discordant for diagnosis of autism to reveal differences in phenotypic severity among and within the subgroups of twins defined by diagnostic status and by zygosity. Then, we used the phenotypic (trait) profiles of each subgroup of twins to perform quantitative trait association analyses to identify SNPs that may be functionally associated with the respective ASD trait (i.e., QTNs). Case-control association analyses were then performed with the QTNs for five autistic behavioral traits using the six different subgroups of twins at the extreme ends (i.e., lowest and highest quartiles) of the respective quantitative trait distribution curves as cases and a large group of unrelated controls. We show here that this novel combination of experimental strategies allowed us to identify significant SNPs that could distinguish individuals with autism from controls. Genes harboring these SNPs are significantly enriched for autism risk genes annotated by SFARI Gene.

#### *3.1. Severity of Autistic Traits Are Strongly Influenced by Genetics*

Numerous studies using both monozygotic and dizygotic twins have sought to determine the role of genetics versus environment with varying results, as described earlier. The quantitative trait distribution profiles generated in this study demonstrate significant differences in severity between clinically discordant non-autistic MZ and DZ co-twins, with the profiles of the MZ twins being significantly more severe than DZ twins across all traits, while the profiles of all diagnosed twins are virtually indistinguishable from each other, regardless of zygosity (Figures 1 and 2). These differences in trait severity between the non-autistic MZ and DZ twins are seen even at the most fundamental level, that is, at the level of specific ADI-R items which comprise each trait, demonstrating the strong influence of genetics on a broad spectrum of ASD-associated behaviors and characteristics (Figure 3). In other words, a non-autistic identical twin exhibits a significantly higher severity of autistic traits in comparison to an undiagnosed non-identical twin who shares only a part of his/her autistic co-twin's genes.

#### *3.2. Quantitative Trait Association Analyses Reveal SNPs that Associate with a Specific ASD Trait*

QTL association analyses allowed us to discover SNPs that may be functionally relevant to autism traits. The genes associated with trait QTNs (Table 1) are over-represented in a number of neurological functions related to autism. The neurological functions associated with each set of trait genes are shown in Tables 4–8. SNPs associated with spoken language implicate genes involved in cell–cell adhesion of neurons, axon guidance and fasciculation, GABA-mediated receptor currents, and morphology of different brain regions, including the cerebral cortex, cerebellum, and corpus callosum, which are all impacted in ASD (Table 4). Another set of SNPs associated with nonverbal communication implicate genes that play a role in cell–cell adhesion of astrocytes, organization and morphology of mossy fibers, axon pruning, migration of pyramidal neurons, and synaptic transmission of collateral synapses (Table 5). SNPs associated with play skills highlight genes involved in dendritic spine formation, synaptic transmission, and long-term potentiation of Purkinje cells (Table 6). Purkinje cell morphology, neurotransmission, axonogenesis, and differentiation of motor neurons are functions over-represented among genes associated with social QTNs (Table 7). QTNs associated with perseverative behaviors point to genes involved in epileptic encephalopathy, loss of hippocampal neurons, and certain higher level functions such as social transmission of food preference, mating, spatial learning, and motor learning (Table 8). Interestingly, the genes associated with the QTNs for perseverative behaviors do not overlap with those associated with any of the other traits, while there is an overlap of 20 genes harboring QTNs for nonverbal, play and social skills (Figure 4), suggesting interrelated proficiencies involving the latter functions. These studies thus help to link genotype to genes and phenotypes and reveal some of the biological processes that may underlie specific ASD traits. The relevance of the above-described neurological functions to ASD further suggested that some of these QTNs may be able to distinguish cases from controls.

#### *3.3. Trait-Associated QTNs Can Discriminate between Cases and Controls*

Case-control association analyses using individuals at the extremes of phenotypic (i.e., trait) severity (Q1 and Q4) as cases enabled the identification of significant trait QTNs that were demonstrated to discriminate between cases and controls. This approach allowed reduction in the phenotypic heterogeneity of the individuals within each subgroup of twins for genetic analyses. In fact, genetic differences between the first and fourth quartiles of trait severity for all twin subgroups are quite clear. Notably, for SNPs associated with both Q1 and Q4 subgroups, the odds ratios for the two extreme groups are distinctly different and often in opposite directions (that is, greater or less than 1) relative to unrelated controls (Tables S7–S12), thereby supporting the associative nature of the specific SNP with that particular trait. Thus, the genetic heterogeneity of cases with respect to these trait-associated SNPs reflects the phenotypic heterogeneity of individuals at the extreme ends of the trait distribution profiles. In addition, hypergeometric distribution analyses of the genes derived from case-control analyses of twin subgroups show significant enrichment in SFARI autism risk genes in three of the twin subgroups (dNA.MZ, dNA.DZ, and cDZ) (Table 3). The over-representation of SFARI-annotated autism risk genes among those revealed here through case-control analyses indicates that the reduction in heterogeneity of cases at the extremes of the trait spectrum allows the detection of significant genetic variants differentiating cases and controls with relatively few samples, while also implicating novel ASD-associated genes that may be replicated by future studies. Although there is precedence for the use of individuals at the extremes of autistic traits for phenotype analyses [24–27], no other study

has utilized individuals within these phenotypically extreme populations for genome-wide genetic association analyses.

#### *3.4. Genetic Di*ff*erences between Discordant DZ Twins Strongly Associate with Autism Risk Genes*

The availability of genome-wide genotype data for both autistic and non-autistic DZ co-twins in our study prompted us to examine the feasibility of performing a case-control analysis in which the group of non-autistic DZ twins was used as a control for the group of autistic co-twins. Although none of the 3400 SNPs with raw *p*-values ≤ 0.05 remained significant after correction for multiple testing, there were 390 SNPs with nominal *p*-values ≤ 0.05 and odds ratios ≥ 3 (ranging from 3 to 11) (Table S13). Hypergeometric distribution analysis of the 146 SNP-associated protein-coding genes showed significant enrichment in ASD-risk genes annotated by SFARI Gene as well as in genes revealed by exome sequencing [23]. Table 13 shows that this set of genes is highly enriched for ASD-associated neuronal functions including development, migration and pathfinding of neurons, neuritogenesis, axonogenesis, and guidance of axons, as well as synaptic development and transmission.

While the use of discordant non-autistic DZ twins as a control group for their autistic co-twins in case-control association analyses may seem unconventional, there is precedence for family-based genetic studies in which concordant MZ twin pairs were compared against their parents by whole genome and exome sequencing, with each study revealing novel and rare variants in the probands [28,29]. The use of both concordant and discordant MZ and DZ twins to investigate the contribution of pre-and post-zygotic de novo CNVs to ASD and other neurodevelopmental disorders [30] further illustrates the value of using diagnostic concordance and discordance in dissecting the complex genetic underpinnings of ASD.

#### *3.5. Advantages and Limitations of Study Design*

An advantage of quantitative trait association analysis is that it allows the identification of SNPs that may potentially be more functionally related to the condition, thus reducing the number of SNPs to be considered for further case-control association analyses. A possible disadvantage is that this procedure may exclude some ASD-relevant but not trait-associated SNPs from case-control analyses. Another limitation is that we did not separately analyze male and female twins due to the relatively small and variable number of females in the different subgroups of twins, and thus cannot eliminate sex as a confounding factor. However, a prior analysis of affected males and females from the AGRE repository using cluster analyses of ADI-R scores as described previously [17] showed no obvious sex differences in subphenotypes [31].

Because of the extreme heterogeneity of ASD, there has been increasing attempts to stratify or subgroup individuals with ASD by specific phenotypes to obtain more homogeneous sample groups for large-scale genetics analyses. While some autism studies have been conducted using individuals at the top 5% of an autistic trait [24,25,27,32], no other study has performed genetic case-control association analyses on subgroups of twins stratified by trait severity in this manner. Although the use of both MZ and DZ twins at the extremes of trait severity reduces the genetic heterogeneity of 'cases' when compared against a large group of unrelated controls, sample sizes at these extremes are a limitation of this study. Another limitation is that there is no empirical validation of the significant SNPs identified by case-control analyses as this was purely an in silico study, and we did not have access to the twins' biological samples for follow-up genotype analyses. On the other hand, the identification of some significant SNPs (with different odds ratios) for samples in the first and fourth quartiles of trait severity in case-control analyses replicates the association of those SNPs with ASD in distinct subgroups of individuals. Despite these limitations, case-control analyses revealed a number of significant trait-associated SNPs that implicate genes that converge on certain neurological functions associated with ASD, suggesting that disruption of those functions or pathways may play a role in the manifestation of that specific trait in ASD. Notably, multiple genes are associated with each of the five traits reflecting the multi-genic nature of ASD, even at these fundamental levels. However, no information is available on the functional significance of the specific SNPs with respect to gene regulation or to protein function and expression, and many SNPs are in noncoding regions which require further exploration.

Case-control association analyses using clinically discordant DZ twins as both cases and controls, while underpowered, identified a large number of nominally significant SNPs associated with genes that are significantly enriched in autism risk genes. While this study shows that reduction of genetic heterogeneity by stratification of cases according to trait severity or by use of discordant DZ twins as cases and controls increases the ability to identify common genetic variants associated with ASD, future studies with greater numbers of individuals in each subgroup are required to replicate these findings.

#### **4. Summary**

We present here a novel experimental design using MZ and DZ twins, either concordant or discordant for ASD, which evaluates the role of genetics in contributing to the severity of autistic traits and identifies significant SNPs that can distinguish cases from unaffected controls by capitalizing on the extreme phenotypic differences exhibited especially by discordant MZ and DZ twins. Although prior studies have consistently indicated that genetics is a major contributor to the overall diagnosis of ASD, deep phenotypic profiling shows that almost all ASD traits are significantly more severe in non-autistic MZ twins than in non-autistic DZ twins, suggesting a dominant influence of genotype on ASD phenotypes probed by the majority of items on the ADI-R diagnostic instrument [16]. At the same time, the similarity of trait-severity profiles among all autistic subgroups suggests that the ADI-R, often considered the gold-standard instrument for autism diagnostics, captures the essential components of ASD across the ASD population.

Quantitative trait association analyses using ASD trait profiles reveal SNPs (QTNs) associated with genes that provide meaningful insight into neurological functions and pathways commonly associated with ASD, thus linking genotype to phenotype (i.e., specific traits) as well as to underlying molecular pathology. Moreover, reducing phenotypic heterogeneity by using individuals at the extremes of trait severity for case-control association analyses with the trait-associated QTNs resulted in the identification of SNPs that significantly discriminate both affected and non-autistic co-twins from a large group of unrelated controls. These results suggest an underlying genetic liability towards autistic traits, even in the undiagnosed individuals (i.e., non-autistic co-twins). Genes implicated by the QTNs that survive correction for multiple testing in case-control association analyses are statistically enriched in ASD-associated genes revealed by significant overlap with those in the SFARI Gene database. Finally, case-control association analysis using dA.DZ twins as cases and their discordant non-autistic co-twins (dNA.DZ) as 'controls' also reveals a large number of neurologically relevant genes that are significantly over-represented in both SFARI Gene and exome sequencing datasets.

#### **5. Materials and Methods**

#### *5.1. ASD Subphenotype Analysis and Generation of Trait Severity Distribution Profiles*

Detailed demographic information of the 284 pairs of twins included in this study is shown in Table S14. Of these sets of twins, 88 monozygotic (MZ) pairs were concordant for diagnosis of ASD, while 25 MZ pairs were discordant. Among the dizygotic (DZ) twins, 56 pairs were concordantly autistic, while 115 pairs were discordant for diagnosis. As mentioned earlier, we label the concordant autistic MZ and DZ twin subgroups as cMZ and cDZ, respectively, and the discordant autistic (dA) and non-autistic (dNA) MZ and DZ twin subgroups as dA.MZ, dA.DZ, dNA.MZ, and dNA.DZ, respectively. ADI-R diagnostic scoresheets for all sets of twins were downloaded from the AGRE repository. It should be noted that ADI-R scoresheets were available for both the autistic and non-autistic (discordant) co-twins in order to establish concordance. Raw scores for 88 items from the ADI-R scoresheet which related to five ASD traits (spoken language, nonverbal communication, play skills, social skills, and perservative behaviors) were extracted from each individual's scoresheet, and the ADI-R scores were

then adjusted according to previously described methods [17]. The specific items related to each trait are shown in Table S15. The average score for each ADI-R item, as well as the cumulative score for a specific trait, was calculated for each of the six twin subgroups. Quantitative trait distribution profiles for each twin subgroup were created by graphing each individual's cumulative score for each trait. To account for differences in the number of individuals in the different subgroups of twins, the total number of individuals in each subgroup was normalized to 100 to facilitate graphical comparisons between trait severity distribution profiles based on zygosity as well as on diagnostic concordance or discordance. The student's two-sample *t*-test in the StatPac statistical software package (https://statpac.com/index.htm) was used to determine significance of differences in severity scores of overall ASD traits and specific ADI-R items based on *p*-values ≤ 0.05.

#### *5.2. Source of Genetic Data*

Genome-wide genotype data for the twins and 2438 unaffected (non-autistic) controls were derived from a previously published genome-wide association study (GWAS) by Wang et al. [33], with data previously cleaned by Jennifer K. Lowe in the laboratory of Daniel H. Geschwind, M.D., Ph.D. at UCLA. Briefly, genotyping data were generated on Illumina HumanHap550 BeadChip arrays with over 550,000 SNP markers. Quality control procedures required DNA samples to have ≥ 0.95 call rate with minor allele frequency ≥ 0.05, and *p*-values testing significant deviation from Hardy–Weinberg equilibrium that were ≥ 0.001. Only subjects of European ancestry were genotyped and thus, only Caucasian twins could be included in the genetic association analyses, eliminating concerns about population stratification as a confounding factor. This publicly available dataset contained the genotypes for 74.3% of MZ twins (85 pairs) and ~85.4% of DZ twins (146 pairs) in this study. Controls were comprised of 2438 unrelated non-autistic individuals of European ancestry for whom genome-wide genotyping data was also available from the same GWAS study [33]. The workflow for the genetic association analyses performed on the twin data is shown in Scheme 1 for clarity.

#### *5.3. Quantitative Trait Loci (QTL) Association Analyses*

PLINK 1.07 software [34] was used to perform all genetic association analyses. Quantitative trait association analyses utilized the cumulative trait scores and the respective genotype data from the above-mentioned dataset for each individual within a specific twin subgroup as previously described [19]. Single nucleotide polymorphisms (SNPs) with unadjusted *p*-value ≤ 1.0 × 10−5 were identified as QTNs which were subsequently used in case-control association analyses as described below. The overlap of trait-associated genes implicated by the SNPs was determined using Venny 2.1.0, an online software package for creating Venn diagrams [35].

#### *5.4. Case-Control Association Analyses*

To reduce ASD heterogeneity for genetic analyses, only cases at the extremes of a defined ASD trait were used in separate case-control association analyses. These individuals exhibited cumulative trait severity scores within the first or fourth quartiles. The first quartile (Q1) was comprised of twins who exhibited the least severe phenotype for that trait and the fourth quartile (Q4) was comprised of twins who exhibited the highest severity for that trait. Cumulative trait scores at Q1 and Q4 for each subgroup of twins were determined using Tyers Box Plot software [36] (see Scheme 1 and Figure 2). Individuals with cumulative scores less than or equal to that determined by Q1 and greater than or equal to that determined by Q4 were used in case-control association analyses against 2438 unrelated controls from the study of Wang et al. [33], focusing on the SNPs discovered by QTL analysis. Discordant non-autistic twins for whom genome-wide genotyping data were available were also stratified and analyzed in the same way. Because only one individual of a pair of MZ twins was included in the original genome-wide genotyping analysis [33], the genotype data of each dA.MZ twin was used for his or her respective dNA.MZ co-twin in case-control analyses. SNPs were considered significant if the adjusted Benjamini–Hochberg FDR was ≤ 0.10.

#### *5.5. SNP Annotation and Functional Analyses of Associated Genes*

Significant SNPs were annotated using SNPper [37] and/or SNPnexus [38]. Both are web-based applications that allow downloading of information on SNPs according to user-defined criteria, such as chromosomal position, band, alleles, closest gene, role (e.g., coding exon, promoter, intron, 3 -UTR), amino acid change, and position. The SNP-associated genes were analyzed using Ingenuity Pathway Analysis (IPA) software (QIAGEN, Redwood City, CA) with a focus on genes involved in neurological functions and pathways. Significance for enrichment of genes associated with a particular pathway or function included in IPA's Knowledgebase was determined by right-tailed Fisher's exact test (*p*-value ≤ 0.05).

#### *5.6. Hypergeometric Distribution Analyses for Gene Enrichment*

SNP-associated genes were then compared to a list of known autism risk genes in the SFARI Gene database [39] and/or a list of ASD-associated genes identified by an exome sequencing study [23]. Hypergeometric distribution probabilities for over-representation of overlapped genes within the SFARI dataset or among ASD genes identified by exome sequencing were calculated using the CASIO Keisan Online Calculator (http://keisan.casio.com/exec/system/1180573201), with significance determined by an upper cumulative *q*-value ≤ 0.05.

#### **6. Conclusions**

Deep phenotyping analysis of monozygotic and dizygotic twins concordant or discordant for diagnosis of ASD indicates that the identical genotype of non-autistic MZ twins with their respective affected co-twins strongly influences the severity of autistic traits in comparison to the trait severity expressed by non-autistic DZ twins. Quantitative trait association analyses using trait distribution profiles and genotype data for the respective twin groups identify SNPs linked to traits. Case-control association analyses using these QTNs and individuals at the extremes of trait distribution profiles reveal significant SNPs associated with genes that are functionally significant to the known pathobiology of ASD. Finally, we suggest that more effective treatment strategies may be derived by associating specific phenotypes with genotype-implicated genes and related functional deficits manifested by a given subgroup of individuals.

**Supplementary Materials:** Supplementary materials can be found at http://www.mdpi.com/1422-0067/20/15/3804/ s1, Table S1: All trait QTNs for cMZ twins, Table S2: All trait QTNs for dA.MZ twins, Table S3: All trait QTNs for dNA.MZ twins, Table S4: All trait QTNs for cDZ twins, Table S5: All trait QTNs for dA.DZ twins, Table S6: All trait QTNs for dNA.DZ twins, Table S7: Case-control by quartile analysis of cMZ twins, Table S8: Case-control by quartile analysis of dA.MZ twins, Table S9: Case-control by quartile analysis of dNA.MZ twins, Table S10: Case-control by quartile analysis of cDZ twins, Table S11: Case-control by quartile analysis of dA.DZ twins, Table S12: Case-control by quartile analysis of dNA.DZ twins, Table S13: Case-control analysis using dA.DZ as cases and dNA.DZ as controls, Table S14: Demographic information for all twins from AGRE repository, Table S15: ADI-R items used for quantitative trait severity profiles and analyses.

**Author Contributions:** Conceptualization, V.W.H.; Phenotype analysis, C.A.D.; Genetic analyses, V.W.H.; SNP annotation and hypergeometric analyses, J.J.D. and V.W.H.; Original draft preparation, V.W.H. and C.A.D.; Review and editing, C.A.D., J.J.D., and V.W.H.; Supervision, V.W.H.

**Acknowledgments:** This research received no external funding. This study was an in silico reanalysis of existing genome-wide genotyping data. The authors are grateful to Zohreh Talebizadeh (Children's Mercy Hospital, Kansas City, MO, USA) and Olivia Veatch (Univ. of Pennsylvania, Philadelphia, PA, USA) for helpful discussion and comments on a draft of the manuscript.

**Conflicts of Interest:** The authors declare no conflict of interest.

#### **References**

© 2019 by the authors. Licensee MDPI, Basel, Switzerland. This article is an open access article distributed under the terms and conditions of the Creative Commons Attribution (CC BY) license (http://creativecommons.org/licenses/by/4.0/).

## **Significantly Elevated** *FMR1* **mRNA and Mosaicism for Methylated Premutation and Full Mutation Alleles in Two Brothers with Autism Features Referred for Fragile X Testing**

**Michael Field 1,**†**, Tracy Dudding-Byth 1,2,**†**, Marta Arpone 3,4, Emma K. Baker 3,4,5, Solange M. Aliaga 3,4, Carolyn Rogers 1, Chriselle Hickerton 3, David Francis 6, Dean G. Phelan 6, Elizabeth E. Palmer 1, David J. Amor 3,4, Howard Slater 3,4,6, Lesley Bretherton 3,4,7,8, Ling Ling 3 and David E. Godler 3,4,\***

Received: 2 July 2019; Accepted: 7 August 2019; Published: 11 August 2019

**Abstract:** Although fragile X syndrome (FXS) is caused by a hypermethylated full mutation (FM) expansion with ≥200 cytosine-guanine-guanine (CGG) repeats, and a decrease in *FMR1* mRNA and its protein (FMRP), incomplete silencing has been associated with more severe autism features in FXS males. This study reports on brothers (B1 and B2), aged 5 and 2 years, with autistic features and language delay, but a higher non-verbal IQ in comparison to typical FXS. CGG sizing using AmplideX PCR only identified premutation (PM: 55–199 CGGs) alleles in blood. Similarly, follow-up in B1 only revealed PM alleles in saliva and skin fibroblasts; whereas, an FM expansion was detected in both saliva and buccal DNA of B2. While Southern blot analysis of blood detected an unmethylated FM, methylation analysis with a more sensitive methodology showed that B1 had partially methylated PM alleles in blood and fibroblasts, which were completely unmethylated in buccal and saliva cells. In contrast, B2 was partially methylated in all tested tissues. Moreover, both brothers had *FMR1* mRNA ~5 fold higher values than those of controls, FXS and PM cohorts. In conclusion, the presence of unmethylated FM and/or PM in both brothers may lead to an overexpression of toxic expanded mRNA in some cells, which may contribute to neurodevelopmental problems, including elevated autism features.

**Keywords:** fragile X syndrome; autism; RNA toxicity; DNA methylation; mosaicism; pediatrics; MS-QMA; AmplideX

#### **1. Introduction**

Fragile X syndrome (FXS) is a leading single-gene cause of inherited intellectual disability (ID), with a prevalence of up to 1 in 4000 [1]. The vast majority (~90%) of males affected by FXS, and ~50% of females, show autism spectrum disorder (ASD) features, including speech perseveration, compulsions, echolalia, repetitive behaviors, poor eye contact, and deficits in social communication [2–5]. Many of these features are shared with idiopathic ASD [6]. The primary molecular cause of FXS is abnormal regulation of the fragile X mental retardation 1 (*FMR1*) gene due to the presence of ≥200 cytosine-guanine-guanine (CGG) repeats, termed full mutation (FM), within the *FMR1* promoter (reviewed in Kraan et al. [7]). These FM alleles are thought to silence *FMR1* transcription [8–12] through epigenetic changes at the *FMR1* promoter, including increased DNA methylation. Loss of *FMR1* mRNA in turn leads to depletion of the fragile X mental retardation protein (FMRP) [13–16], which is the primary cause of FXS symptomology, as FMRP has a critical function in synaptic plasticity and normal brain development [17–21].

In contrast, smaller CGG alleles (55–199 repeats), termed premutation (PM), usually have an unmethylated *FMR1* promoter, but abnormally increased levels of expanded PM mRNA [22,23]. In some PM carriers, the transcription of *FMR1* mRNA from expanded alleles has been associated with "RNA gain of function" toxicity, implicated in late onset disorders such as Fragile X-associated Tremor/Ataxia Syndrome (FXTAS) [22,23]. Other factors proposed to contribute to PM-related disorders are expanded repeat associated non-AUG translation, and increased transcription of *ASFMR1*/*FMR4* originating from the same locus as *FMR1* (reviewed in Kraan et al. [7]). In FXS, these mechanisms have not been comprehensively studied, but may explain significant variability in the type and severity of the FXS phenotype, beyond FMRP deficiency. This is consistent with our recent study suggesting that, in the majority (60%) of FXS males, FM mRNA can still be detected, and that this incomplete silencing in males <19 years of age is associated with more severe autism features, but not more severe ID [6]. This led us to propose that FXS can be stratified based on complete and incomplete silencing, with the two reciprocal mechanisms of RNA toxicity and FMRP deficiency contributing to overlapping aspects of FXS, namely the ID and autism phenotype in each individual that has incomplete silencing. This is particularly important as mosaicism for active and inactive expanded alleles may be more common than previously reported (10% to 40%) [5,11,15,24–29], based on the detection of *FMR1* mRNA in more than 60% of FM males [6]. The difference of >20% in the reported prevalence of mosaicism in FXS by Southern blot analysis as opposed to incomplete *FMR1* silencing detected by real-time PCR [6], is likely due to the limited analytical sensitivity of methylation-sensitive Southern blot analysis (the 'gold standard' FXS diagnostic test used in earlier studies) that does not detect mosaic alleles in less than 20% of cells [26]. Extreme examples of low level methylation mosaicism found in <20% of cells have previously been described in rare adult individuals with unmethylated FM (UFM) alleles that typically have mild ID (but are at risk of premutation phenotypes such as FXTAS) [30].

This study, for the first time, describes two brothers with UFM alleles detected by Southern blot analysis in a pediatric setting. Based on non-verbal intellectual functioning assessments, the brothers are 'higher functioning' compared to typical FXS, but have autism features ranging from mild to severe. The detailed clinical and molecular follow-up is described. This includes the detection of low-level mosaicism for methylated PM and FM alleles in different tissues. The CGG sizing and *FMR1* mRNA analysis of the two brothers are also compared with family members ascertained through cascade testing, typical FXS, age-matched children with PM, and typically developing (TD) controls.

#### **2. Results**

#### *2.1. Medical History*

B1 was born at term after an uncomplicated pregnancy and delivery. His parents became concerned with his development from an early age. He had extreme shyness, disliked new social circumstances, avoided some physical contact, and had limited eye contact. He was difficult to settle and he was often irritable and upset. His early motor milestones were normal, but language was significantly delayed. A clinical evaluation with the Griffiths Mental Developmental Scales, Extended Revised led to a diagnosis of mild global developmental delay, with significant difficulties with language and play skills. Fragile X testing by Southern blot analysis of DNA extracted from blood and microarray analysis were performed at the age of 3 years. Southern blot analysis detected an abnormal 1.42 kb pfxa3 fragment equivalent to approximately 170 CGG repeats expansion. No pathogenic copy number variants (ISCA 60K) were identified through a chromosomal microarray. The results of psychometric testing (Stanford-Binet Intelligence Scales 5th Edition) [31] at the age of 4.5 years were consistent with a mild ID, with relative strengths in nonverbal abilities. The Adaptive Behavioural Assessment System (ABAS-II) results indicated a mild deficit in everyday functioning. He was treated with fluoxetine, which improved his social anxiety. After B1's fragile X genetic testing, cascade testing was carried out for other family members. The mother of both brothers was found to have a 74 CGG expanded allele; the maternal half-sister, maternal aunt, and maternal grandfather were also all identified as having a PM allele (Figure 1A).

**Figure 1.** (**A**) Family pedigree, with expansion sizes from AmplideX testing. None of the individuals had manifested symptoms consistent with fragile X syndrome or fragile X-associated disorders, such as Fragile X-associated Tremor/Ataxia Syndrome (FXTAS) or Fragile X-associated Primary Ovarian Insufficiency (FXPOI). The brothers' have two half-sisters, one of which was identified as carrying a premutation (PM) allele and the other as carrying a normal fragile X mental retardation 1 (*FMR1*) allele. B1 s cousin (C1) was identified at 1 year and 11 months of age as carrying an *FMR1* PM allele of 66 CGG repeats. (**B**) Methylation-sensitive Southern blot analysis of the *NruI* restriction site within the *FMR1* CpG island in blood. The DNA samples from the two brothers in question are located in lanes 5 and 6. B2 has a 223 CGG unmethylated allele, while B1 has a 190 CGG repeat unmethylated allele. A typical female (CGG < 44) is in lane 1, while typical males (CGG < 44) are in lanes 3 and 4; and a full mutation (FM) female with a 100% methylated 580 CGG allele is in lane 2. The numbers superimposed on the blot indicate CGG sizing. The lower limit of detection for the Southern blot analysis of FM alleles on a normal allele background is 20%. (**C**) Photos of B1 and B2 at the ages of 5 and 2 years, respectively; written informed consent was obtained from the mother of the brothers to include these images in the publication (**D**) Photos of B1 and B2 at the ages of 10 and 7 years, respectively; written informed consent was obtained from the mother of the brothers to include these images in the publication. Note: Comparator DNA used of methylation-sensitive Southern blot analysis in (**B**) were also sized using standard CGG sizing PCR.

B2 is the younger brother of B1, who is a dizygous twin. The twins were conceived before the *FMR1* CGG expansion was identified in B1. B2 and his twin brother were delivered at term by elective caesarean section, in good condition after a pregnancy complicated by gestational hypertension. The twins were admitted for two days in the newborn nursery for blood sugar instability and required intravenous fluids. However, there were no ongoing neonatal concerns and the boys were discharged home on day four of life. There were no concerns about B2 s development in the first year of life; he fed and slept well and all his motor milestones were in the normal range. He had surgery for left cryptorchidism, but his medical history was otherwise unremarkable. B2 s mother became concerned about his development and behaviour around two years of age. Although B2 did not have typical craniofacial features of fragile X syndrome, he had a slightly longer face than his twin brother, which was visible at the time of assessment (Figure 1C) and also more recently (Figure 1D). He was normocephalic, had a normal tone and power, with a normal gait, and did not have joint hypermobility. At 9 months of age, fragile X PCR with AmplideX commercial kit was performed on their DNA extracted from blood, which identified a normal triplet repeat size (29 CGG repeats) in one of the twins, but a 174 PM CGG repeats expansion in B2.

#### *2.2. Extended FMR1 CGG Testing*

Extended *FMR1* testing was initiated to test the hypothesis that the neurodevelopmental phenotypes in B1 and B2 were due to the presence of FM alleles in tissues other than blood. For these extended studies, B1 and B2 were recruited into the FREE FX study at 5.6 and 2.7 years of age, respectively, but were not included in previous FREE FX study publications due to uncertainty about molecular diagnosis [6].

In blood, methylation-sensitive Southern blot analysis demonstrated the presence of unmethylated alleles of 190 CGGs in B1 and 223 CGGs in B2 (Figure 1B). These results are consistent with those previously reported for rare 'high functioning' FXS males with fully unmethylated FM alleles [32–34]. AmplideX testing of blood at the time of the FREE FX study assessments showed an allele of 168 CGGs in B1, while for B2, there were three alleles detected of 165, 185, and 199 CGGs (Figure 2B). AmplideX CGG sizing analysis in tissues other than blood showed that for B1, CGG sizing in saliva was consistent with that in blood. However, in fibroblasts, the sizing was different, with mosaicism for 140 CGG and 156 CGG alleles detected. In contrast, for B2, all alleles in the PM range were consistent between the tissues tested. However, for B2, in buccal epithelial cells (BEC) and saliva, FM alleles were also detected, which were not identified in blood. Subsequently, B2 received a molecular diagnosis of FXS that was not initially given by standard testing of blood.

The variability in CGG sizing between tissues was also investigated in both brothers using FastFraX 5 and 3 melting curve analysis (MCA) assays, using samples from C1 (cousin, as reference) [35]. Surprisingly, there was no amplification by the 3 MCA assay in B1 and B2, while amplification occurred by 3 MCA for C1, with the 5 MCA performed on the same samples from B1, B2, and C1 displaying amplification, as expected (Supplemental Figure S1). This suggested that one of the primer sites for the 3 MCA assay had a de novo sequence change of unknown significance in all tissues tested in both brothers that prevented primer binding and amplification, but not in their cousin. As a follow-up, sanger sequencing of the 3 MCA *FMR1* exon 1 binding site inclusive of the ATG translation start site, and surrounding regions, was performed. However, no sequence variants were detected in both brothers compared to controls (Supplemental Figures S3 and S4). This suggested that a de novo sequence change was present for the binding site of the primer anchored at the 3 end of the CGG repeat, rather than the primer overlapping with the ATG site within the exon 1.

#### *2.3. FMR1 Methylation and Gene Expression Analyses*

Methylation analysis using Methylation-Specific Quantitative Melt Analysis (MS-QMA) identified methylated alleles in both brothers in multiple tissues (Figure 3A,B) that were not present in 17 control males and 14 PM males co-run with these samples. However, the methylated peaks of ~5% in both brothers did overlap with those found in 41 FM only males and 18 typical PM/FM mosaics also co-run with the B1 and B2 samples (Supplemental Figure S2). The levels observed in both brothers were also consistent with those of the spiking reference sample, which contained 6% methylated FM and 94% unmethylated normal size male DNA (Figure 3C). Presence of the unmethylated expanded alleles in the majority of cells was consistent with *FMR1* mRNA analyses in blood for both brothers (Figure 4), where the levels in both brothers were ~5 fold higher than the levels observed in male and female controls. Interestingly, of the 13 PM male reference samples (most children or adolescents), only one had analogous *FMR1* mRNA levels. None of the FM only or PM/FM mosaic males had mRNA levels remotely close to those observed in the brothers, with the C1 PM (66 CGGs) cousin having mRNA levels just above the control range.

**Figure 2.** (**A**) Organization of the fragile X mental retardation 1 (*FMR1*) 5 region, including the cytosine-guanine-guanine (CGG) expansion and primer binding site locations for AmplideX CGG sizing PCR, FastFraX High-Resolution Melt Analysis, and Methylation-Specific Quantitative Melt Analysis (MS-QMA), in relation to Fragile X-Related Epigenetic Element 2 (FREE2); the *FMR1* CpG island and methylation sensitive restriction site *NruI* analysed using routine fragile X Southern blot testing; and two *HpaII* sites targeted by AmplideX methylation PCR. Note: The *FMR1* translation start site is indicated in red overlapping with FastFraX 3 melting curve analysis (MCA) primer. (**B**) AmplideX triplet-repeat primed long range PCR results for B1 and B2, for blood, buccal, saliva, and fibroblast DNA. Note: a reference sample co-run with the samples in question is included as the AmplideX sizing control in (**B**).

**Figure 3.** Derivative curve high-resolution melt profiles and mean methylation output ratios of CpG sites located within the Fragile X-Related Epigenetic Element 2 (FREE2) region assessed using Methylation-Specific Quantitative Melt Analysis (MS-QMA) between different tissues for (**A**) B1 and (**B**) B2.  These samples were mixed at different ratios for the MS-QMA methylation reference curve, with the expected % of methylation indicated on the plot for each derivative curve profile.

**Figure 4.** Comparison of fragile X mental retardation 1 (*FMR1*) mRNA levels in blood in the probands, control, premutation (PM), full mutation (FM) only, and PM/FM mosaic reference cohorts. Note: Control reference ranges included in blood of females (age 22 to 54 years; cytosine-guanine-guanine (CGG) < 44) and males (age 7.7–8.1 years) typically developing (TD) controls; males with PM alleles (3.4 to 23.6 years); and males with FM only alleles and PM/FM mosaic alleles (1.89 to 43.17 years). Control children had no family history of developmental delay and their mothers had CGG size < 44.

#### *2.4. Neurodevelopmental Outcomes at the Time of Recruitment*

In B1, the cognitive assessment as part of the FREE FX study highlighted receptive language difficulties, which greatly influenced his performance in the IQ test. He showed significant variability in his cognitive abilities. His verbal skills were a cognitive weakness, reflected in his verbal IQ (VIQ) score of 72 falling in the borderline range. His overall nonverbal problem-solving skills were in the low average range, with a performance IQ (PIQ) score of 84. His processing speed index (PSI) score was 71, which fell in the borderline range. His overall intellectual functioning, as reflected in the Full Scale IQ (FSIQ), was 77. Based on the Autism Diagnostic Observation Schedule Second Edition (ADOS-2) assessment, B1 met the cut-off for autism spectrum, with an overall ADOS-2 calibrated severity score (CSS) of 4. His Social Affect (SA) CSS was 3 (non-spectrum range), whilst his Restricted and Repetitive Behavior (RRB) CSS was 6 (autism range). B1's scores on the Child Behavior Checklist (CBCL) were in the normal range for the Total Problem and Externalizing problem scales, but fell in the borderline range for the Internalizing problems scale. Moreover, the Pervasive Developmental Problems subscale score fell in the clinical range and the Withdrawn subscale score fell in the borderline range. At the time of the research assessment, B1 was enrolled in a small mainstream primary school with additional learning support. He was reported to play alongside other children and although still shy, he did not have any major challenging behaviours. He had no seizures and was taking 2 mg/day of sertraline.

During assessment, as part of the FREE FX study, B2 predominantly used non-word vocalizations and jabbering; although on a number of occasions, he used single words and word approximations. Moreover, his formal developmental assessment with the Mullen Scales of Early Learning (MSEL) showed that his receptive and expressive language skills were respectively similar to that expected of an

infant aged 10 and 17 months old, a marked discrepancy from his chronological age (31 months). He did not obtain a valid score (T ≥ 20) on the receptive language subdomain of the MSEL; a further indication of the significant impairment in this area. B2 s ratio VIQ [36] was 44. In contrast, his performance in tasks assessing visual reception and fine motor skills was much better, obtaining age equivalents of 25 and 20 months, respectively, and a ratio non-verbal IQ (NVIQ) of 73. His overall MSEL Early Learning Composite score (using a default minimum T score of 20 for his Receptive Language domain) was 56, more than three standard deviations below the mean score (mean = 100; SD = 15) of similar aged peers, indicating a significant developmental delay. B2 s difficulties with social communication, his engagement in sensory seeking and repetitive behaviours, and restricted interests observed during the ADOS-2 assessment led to an ADOS-2 classification of autism, with an overall ADOS CSS of 7. His SA CSS score of 5 fell in the autism spectrum range, whereas his RRB CSS of 10 fell in the autism range, the highest severity level for this domain for children with ASD of his age and verbal abilities. B2 s CBCL Total, Internalizing, and Externalizing problem scales scores were all in the normal range. However, his score in the Withdrawn and Pervasive Developmental Problems subscales fell in the clinical range. His mother reported significant concerns regarding his very limited speech and social anxiety.

#### **3. Discussion**

This study describes two young brothers with expanded *FMR1* alleles, who were 'higher functioning' based on intellectual functioning assessments compared to a typical FXS cohort reported in previous studies [6]. Specifically, B1 had a PIQ score three standard deviations above the mean of the FREE FX typical FXS male cohort, while B2 had a non-verbal (NVIQ) score that was nearly two standard deviations above the same mean. Regarding VIQ, B1 s score was 1 SD greater than the FXS group mean, while B2 fell within 1 SD of the typical FXS cases on this measure. B2 s greater difficulties with expressive and receptive language skills may be closely linked to his more severe autism phenotype. In contrast, B1 presented mild autism features, though still met the ADOS-2 cut-off for autism spectrum. Based on methylation-sensitive Southern blot analysis of blood DNA, both brothers have been described as having UFM alleles approaching 200 CGGs in blood. Previous studies have described adults with similar Southern blot profiles, explaining their higher intellectual functioning through the *FMR1* promoter being completely unmethylated and expressing *FMR1* FM mRNA and FMRP [32–34]. These adult males with unmethylated alleles (FM and PM/FM mosaic) and incomplete silencing of *FMR1* mRNA from expanded FM alleles are at risk of FXTAS, based on clinical assessments and magnetic resonance imaging (MRI) features [30,37,38]. In these studies, it was hypothesized that expression of the UFM RNA is a contributor to the neurodegenerative features observed in some adults, through the same RNA toxicity mechanism as in PM-related disorders. However, the interplay between intellectual functioning and autism severity in the pediatric setting has not previously been characterized in individuals with UFM alleles.

In this study, additional molecular analyses with more sensitive techniques, such as MS-QMA and AmplideX PCR, indicate that a small proportion of cells in blood and other tissues have methylated PM and FM alleles that were not detected by Southern blot analysis. These findings are partly consistent with a previous study examining tissue heterogeneity in *FMR1* methylation and CGG size post-mortem, in a 79-year-old 'high-functioning' male with unmethylated PM and FM alleles detected in blood by Southern blot analysis [34]. In this male, a complete unmethylated PM/FM allele smear was found in blood by Southern blot analysis; however, in multiple other tissues, including the parietal lobe, methylated FM alleles were detected. Together, this suggests that the detection of methylated alleles varies with the tissue tested and the analytical sensitivity of the technique used. With the addition of more sensitive technologies to complement Southern blot analysis, the UFM classification in some of these cases is likely to be replaced with the term 'low-level methylation mosaicism'. Moreover, these males with low-level somatic mosaicism may not be as uncommon as previously thought. One reason for this may be that these unmethylated alleles are somatically unstable, with a significant

proportion of these cases that have retraction down to a normal repeat size currently not being detected by standard testing [26].

It is also important to note that based on CGG sizing and methylation analysis of multiple tissues, B1 was found to have multiple methylated alleles in the PM range in a small proportion of cells (e.g., 170 CGG repeats in blood, but 140 and 156 CGG alleles in fibroblasts). This is consistent with an earlier study [26] that also reported a male with a methylated PM allele of 110 CGG repeats, in the absence of an FM, with a FSIQ value of 46. Together, this suggests that as part of somatic retraction below an FM allele size, some alleles may remain methylated in a proportion of cells, and may not be fully functional, with this proportion and CGG sizing differing between tissues. For the two brothers, however, the RNA toxicity, rather than FMRP deficiency, is likely to be the primary factor contributing to the autism phenotype. This is consistent with FM RNA levels being ~5 fold greater than in the reference sample of typically developing controls, PM/FM mosaic and FM only males affected with FXS, and the PM cousin (C1) with 66 CGGs. This is also consistent with our earlier study showing that FM males with incomplete silencing of *FMR1* had more severe autism features compared to the FM group where *FMR1* mRNA could not be detected [6].

This study also detected a potential sequence variant in the *FMR1* promoter overlapping with the 3 MCA FastFraX primer binding site at the 3 end of the CGG repeat. While the functional and clinical significance of this potential variant is uncertain, it did not inhibit *FMR1* transcription in both brothers. Similarly, previous studies have reported a sequence variant located next to the CGG repeat in FXS individuals displaying somatic retractions [30]. The loss of primer binding sites that may be associated with somatic retraction is an important limitation of all standard and long-range PCR-based methods, including AmplideX and FastFraX, as it may lead to expanded alleles being missed or not amplified, as demonstrated by the 3 MCA assay for both brothers in this study, and for AmplideX previously [30].

In summary, the study reports two brothers with low-level methylation mosaicism not detected by standard testing in blood, who may be mistaken for males with rare UFM alleles, when based alone on Southern blot analysis of blood DNA. Moreover, in both brothers, *FMR1* mRNA levels were increased ~5-fold compared to typical developing controls, and significantly above the levels reported from PM, FM only, and PM/FM mosaic male cohorts from a previous study [6]. The abnormally elevated levels of *FMR1* mRNA in the two brothers may have led to *FMR1* mRNA-related cellular "toxicity". It is plausible that an FXS patho-mechanism exists whereby active unmethylated FM and/or PM alleles lead to the expression of toxic expanded mRNA in some cells, in conjunction with possible reduced *FMR1* mRNA and FMRP levels in other cells with *FMR1* methylation detected by MS-QMA. A combination of these two mechanisms in different cells may contribute to the brothers' neurodevelopmental problems, including the elevated ASD symptoms. This hypothesis is also in line with the findings reported by Baker et al. [6] showing that males with FM-only CGG expansions (aged < 19 years), who expressed FM *FMR1* mRNA, had significantly more severe ASD symptoms measured with the ADOS-2 compared to males with FM-only alleles who had completely silenced *FMR1*. The main limitations of this study are that the FMRP levels in blood and direct sequencing of the 3 MCA binding site, anchored at the 5 end of the CGG repeat, suspected to be modified in both brothers, have not been tested and performed. Future studies will address these limitations to shed light on the effect of the potential 3 MCA variant (to be confirmed through DNA sequencing) on FMRP structure and function, and possibly further explain the phenotype presented by these two young children.

#### **4. Materials and Methods**

#### *4.1. Ethics Approval and Consent*

All aspects of this study have received ethical approval by The Royal Children's Hospital Human Research Ethics Committee (Single Site Reference numbers: HREC 34227, HREC 33066—Multi site HREC Reference Number: HREC/13/RCHM/24, approved on 24 May 2013).

#### *4.2. Recruitment and Assessments*

Affected family members include two brothers, B1, age 5.6 years, and B2, age 2.7 years, and their cousin C1, age 3.4 years. All were recruited through their clinical geneticist into the FREE FX study [6]. As part of the FREE FX study, the brothers' neurodevelopmental outcomes were evaluated through a medical and developmental history questionnaire prepared ad hoc for this study, the ADOS-2 (Module 2 for B1 and Module 1 for B2) [39], the Wechsler Preschool and Primary Scale of Intelligence–Third Edition (WPPSI-III) (Australian) (B1) [40], MSEL (B2) [41], and CBCL [42]. Both brothers' and their cousin's buccal epithelial cells, saliva, and venous blood samples were collected at the time of their participation in the study. A skin biopsy for genetic testing on fibroblasts was performed on B1, as part of his clinical care. Molecular data for positive and negative control cohorts used as reference data in this study were analysed previously [6,43,44] and from a newly recruited control and PM participants as part of the FREE FX study.

#### *4.3. Sample Processing*

Up to four BEC samples were collected per participant using the Master Amp Buccal Swab Brush kit (Epicentre Technologies, Madison, WI, USA), as previously described [45]. A single saliva sample was collected for each participant using the Oragene® DNA Self-Collection Kit (DNA Genotek, Global) and was processed as per the manufacturer's instructions. As part of a clinical genetic follow-up, and not as part of his involvement in the FREE FX study, a skin biopsy was performed for B1 to obtain fibroblasts. The DNA extracts were evaluated using a NanoDrop 2000 spectrophotometer (Thermo Fisher Scientific, Foster City, CA, USA). Ten ml of venous blood, from B1 and B2 and from two control male participants (aged 7.6 and 8.1 years, respectively), was used for peripheral blood mononuclear cell (PBMC) isolation using Ficoll gradient separation [46]. Isolated PBMCs were used for RNA extraction (RNeasy kit; Qiagen Inc., Hilden Germany) for gene expression analyses [46]. Total RNA extraction, purification, and reverse transcription were performed as previously described [47].

#### *4.4. Methylation Specific-Quantitative Melt Analysis (MS-QMA)*

DNA samples were extracted from BEC, saliva, venous blood, and fibroblasts, and transferred into 96-well plates to be treated with sodium bisulphite. An EZ DNA Methylation-GoldTM kit (Zymo research, Irvine, CA, USA) was used to bisulphite convert each sample in two separate reactions, with each conversion analysed in duplicate reactions. FREE2 DNA methylation analysis was performed using MS-QMA, as previously described [43]. Specifically, ninety-six samples were bisulfite converted at a time (3 controls and 93 unknown samples per plate) and were serially diluted four times post-conversion. These included positive and negative control cohorts used as reference data in this study, as well as samples in question from B1, B2, and C1. The bisulfite converted DNA was then transferred into a 384 well format for real-time PCR analysis utilizing MeltDoctor™ high-resolution melt reagents in 10 μL reactions, as per the manufacturer's instructions (Life technologies, Foster City, CA, USA). A unique primer set was used for real-time PCR that targets specific CpG sites within the Fragile X-Related Epigenetic Element 2 (FREE2) region, at the *FMR1* exon1/intron 1 boundary [43]. The annealing temperature for the thermal cycling protocol was 650 ◦C for 40 cycles. The ViiA™ 7 Real-Time PCR System (Life technologies, Foster City, CA, USA) was then used to quantify the DNA concentration of the unknown samples using the relative standard curve method post bisulfite conversion, by measuring the rate of dye incorporation into double stranded DNA. To progress to the next stage of the MS-QMA analysis, the unknown samples had to be within this dynamic linear range.

The products from a methylated and unmethylated FREE2 sequence were then separated into single strands in the temperature range of 74 and 82 ◦C as part of the high-resolution melt step that followed the real-time PCR in a close tube format. The HRM Software Module for ViiA™ 7 System was then used to plot the rate of PCR product separation to single strands, with the difference in fluorescence converted to Aligned Fluorescence Units (AFU).

The AFU conversion to the methylation percentage was performed at 78 ◦C, and all of the above quality control steps, were analysed simultaneously for 384 reactions at a time using Q-MAX software (Curve Tomorrow, Melbourne, Australia), developed to automate the process.

This software utilizes a custom-designed computer algorithm to simultaneously perform multiple quality control checks to determine DNA concentrations and quality post-bisulfite conversion using raw RT-PCR and high-resolution melt data for all dilutions from each bisulfite reaction. The high-resolution melt data for those sample dilutions outside the QC ranges were discarded from the quantitative methylation analysis by Q-MAX software, and were not used for the final aggregate methylation ratio calculation. The high-resolution melt profiles discarded from quantitative assessments by Q-MAX, however, were used for visual assessments for the presence or absence of abnormal methylation compared to control melt curves (Figure 3). Males were considered to have positive MS-QMA results if they had >2% methylation derived through Q-MAX and/or one of more high-resolution melt curve profiles that had derivative high-resolution melt plots with the presence of a peak originating from methylated alleles, at the same melting temperature as in FXS male controls (Figure 3).

#### *4.5. Methylation Analysis and CGG Sizing Using Southern Blot and AmplideX PCR*

Methylation-sensitive Southern blot analysis of the *NruI* restriction site within the *FMR1* CpG island was performed for the two brothers' venous blood DNA samples, using a fully validated methylation-sensitive Southern blot procedure with appropriate normal and abnormal controls, as described previously [48]. The *FMR1* CGG repeat size was assessed using a fully validated PCR assay with a precision of +/− one triplet repeat across the normal and grey zone ranges, performed using a fragment analyser (MegaBACE, GE Healthcare, Chicago, Illinois, IL, USA), with the upper limit of detection of 170 repeats, as described previously [49]. CGG sizing using both brothers' BEC, saliva, venous blood, and fibroblast (only B2) DNA was also performed using the AmplideX®™ *FMR1* PCR Kit as per the manufacturer's instructions [50] (Asuragen, Austin, TX, USA). The *FMR1* CGG size in DNA extracted from venous blood from the two brothers was further investigated by Southern blot testing at the Victorian Clinical Genetics Services (VCGS) (without the use of methylation sensitive restriction enzymes), as described previously [48].

#### *4.6. Melting Curve Analysis (MCA) to Determine the Presence of Expanded FMR1 Alleles*

Five prime and 3 MCA high-resolution melt analyses were also undertaken in the BEC, saliva, and venous blood DNA of both brothers, in the fibroblasts DNA of B1, and in the venous blood DNA of C1, as per the manufacturer's instructions for the FastFraX commerical kit (Biofactory, Singapore) [35]. Specifically, 5 and 3 PCR reactions were performed for each sample, with each assay having one of the primers anchored at one end of the CGG repeat. Each assay contained 5 units of HotStarTaq DNA polymerase (Qiagen), 2.5× Q-Solution (Qiagen), 1× of the supplied PCR buffer (Qiagen), 0.1× SYBR Green I nucleic acid dye (Roche Applied Science, Upper Bavaria, Germany), and 50 ng genomic DNA. The 5 PCR used a deoxynucleoside triphosphate mix consisting of 0.2 mmol/L each of dATP, dTTP, and dCTP, and 0.1 mmol/L each of 7-deaza-2 -dGTP (7-deaza-dGTP) and dGTP (Roche Molecular Diagnostics, Upper Bavaria, Germany). The 3 -PCR used a deoxynucleoside triphosphate mix consisting of 0.2 mmol/L each of dATP, dTTP, dCTP, and dGTP, and primer sequences and concentrations as previously described [35], with primer locations indicated in Figure 2A. The 5 and 3 -PCR reactions were performed separately under identical thermocycling conditions in a ViiA™ 7 Real-Time PCR System (Life technologies, Foster City, CA). An initial denaturation step at 95 ◦C for 15 min was followed by 40 cycles of 99 ◦C for 2 min, 65 ◦C for 2 min, and 72 ◦C for 3 min, and then a final extension step at 72 ◦C for 10 min. PCR amplicons were then melted (after completion of the PCR program), consisting of denaturation at 95 ◦C for 1 min, a temperature-hold step at 60 ◦C for 1 min, and a temperature ramp from 60 ◦C to 95 ◦C at a rate of 0.01 ◦C/s. Reference male samples were co-run with the samples in question, including an FM (530 CGG), a PM (170 CGG), and a normal size (NS) control (30 CGG). Positive and negative calls were made based on the presence or absence of the

difference in the profile fluorescence at the melting temperature threshold (85 ◦C for 5 MCA; 90 ◦C for 3 MCA) compared to male reference samples.

#### *4.7. Sanger Sequencing*

Sanger sequencing of the 3 MCA *FMR1* exon 1 binding site inclusive of the ATG translation start site, and surrounding regions was performed at the Australian Genome Research Facility Ltd in blood of B1 and B2, and compared to 2 reference samples from typically developing controls co-run with these samples. The data was re-analysed at the Victorian Clinical Genetics Services (VCGS), using Mutation Surveyor®V4.0.9 software (SoftGenetics®, State College, PA, USA) to detect and report SNPs as previously described [51].

#### *4.8. FMR1 mRNA Analysis*

Gene expression analyses were performed using reverse transcription real-time quantitative PCR (RT-PCR) on a ViiaTM 7 System (Life Technologies, Global), with the relative standard curve method as described in Kraan et al. [52]. The mean 5'and 3 *FMR1* mRNA levels were normalized to the mean of *EIF4A2* and *SDHA* mRNA levels used as internal controls, expressed in arbitrary units (a.u) [52]. The summary measure used for mRNA expression levels for each participant was calculated by averaging the four arbitrary unit outputs originating from the two separate cDNA reactions, with each of these analysed in two separate RT-PCR reactions, performed for each RNA sample [52,53].

**Supplementary Materials:** The following are available online at http://www.mdpi.com/1422-0067/20/16/3907/s1.

**Author Contributions:** Conceptualization: M.F., T.D.-B., D.J.A., and D.E.G.; data curation: M.A., E.K.B., S.M.A., D.G.P., and L.L.; formal analysis: M.F., T.D.-B., M.A., E.K.B., E.E.P., and L.B.; funding acquisition: M.F., H.S., D.J.A., and D.E.G.; investigation: M.F., T.D.-B., M.A., S.M.A., C.R., C.H., D.F., E.E.P., L.L., D.G.P., and D.E.G.; resources: M.A., H.S., and D.E.G.; supervision: M.F., H.S., L.B., and D.E.G.; visualization: D.E.G.; writing—original draft preparation: M.F., T.D.-B., M.A., E.K.B., D.J.A., and D.E.G.; All authors have reviewed and provided corrections, and have agreed to be co-authors on the submitted version of the paper.

**Funding:** This work was funded by The Victorian Government's Operational Infrastructure Support Program, Murdoch Children's Research Institute, Royal Children's Hospital Foundation, Martin & E.H. Flack Trust, Pierce Armstrong Trust, Financial Markets Foundation for Children (Australia) (FMFC; grant number: 2017-361), and the National Health and Medical Research Council (NHMRC; project grant numbers: 1049299 and 1103389). D.E.G. was supported by the Next Generation Clinical Researchers Program—Career Development Fellowship Funded by the Medical Research Future Fund (grant number 1141334). M.A. was supported by the International Postgraduate Research Scholarship (IPRS) and the Research Training Program Fee offset scholarship funded by the Australian Government and awarded by the University of Melbourne, and in part by the Diagnosis and Development group of the Murdoch Children's Research Institute.

**Acknowledgments:** We sincerely thank the family described in this study and all participants and their parents/guardians for their contribution, support and generosity, which have been greatly appreciated. We also thank Annabel May Marsh for assisting with participants' psychological assessments.

**Conflicts of Interest:** David E. Godler is named as an inventor on patent applications (PCT/AU2010/000169 and PCT/AU2014/00004) related to the technology described in this article. The other authors declare that the research was conducted in the absence of any commercial or financial relationships that could be construed as a potential conflict of interest.

#### **References**

© 2019 by the authors. Licensee MDPI, Basel, Switzerland. This article is an open access article distributed under the terms and conditions of the Creative Commons Attribution (CC BY) license (http://creativecommons.org/licenses/by/4.0/).

*Article*

## **Associations between Monocyte and T Cell Cytokine Profiles in Autism Spectrum Disorders: E**ff**ects of Dysregulated Innate Immune Responses on Adaptive Responses to Recall Antigens in a Subset of ASD Children**

**Harumi Jyonouchi 1,2,\* and Lee Geng 1**

Received: 26 July 2019; Accepted: 23 September 2019; Published: 24 September 2019

**Abstract:** Changes in monocyte cytokine production with toll like receptor (TLR) agonists in subjects with autism spectrum disorders (ASD) were best reflected by the IL-1β/IL-10 ratios in our previous research. The IL-1β/IL-10 based subgrouping (low, normal, and high) of ASD samples revealed marked differences in microRNA expression, and mitochondrial respiration. However, it is unknown whether the IL-1β/IL-10 ratio based subgrouping is associated with changes in T cell cytokine profiles or monocyte cytokine profiles with non-TLR agonists. In ASD (*n* = 152) and non-ASD (*n* = 41) subjects, cytokine production by peripheral blood monocytes (PBMo) with TLR agonists and β-glucan, an inflammasome agonist, and T cell cytokine production by peripheral blood mononuclear cells (PBMCs) with recall antigens (Ags) (food and candida Ags) were concurrently measured. Changes in monocyte cytokine profiles were observed with β-glucan in the IL-1β/IL-10 ratio based ASD subgroups, along with changes in T cell cytokine production and ASD subgroup-specific correlations between T cell and monocyte cytokine production. Non-ASD controls revealed considerably less of such correlations. Altered innate immune responses in a subset of ASD children are not restricted to TLR pathways and correlated with changes in T cell cytokine production. Altered trained immunity may play a role in the above described changes.

**Keywords:** ASD; cytokine; monocyte; β-glucan; T cell cytokine; trained immunity; maternal immune activation

#### **1. Introduction**

Autism spectrum disorder (ASD) is a behaviorally defined syndrome and the effects of genetic and environmental factors that form its pathogenesis, likely vary in individuals with ASD diagnosis. However, recent research is increasingly supporting the role of immune mediated inflammation in its pathogenesis [1–3]. In fact, one of the most extensively studied animal models of ASD is generated by inducing sterile inflammation in pregnant rodents through intraperitoneal injection of stimulants of innate immunity, such as endotoxin [4]. This model which is called maternal immune activation (MIA) has been shown to create lasting neurodevelopmental and behavioral changes in offspring, as well as lasting effects on T cell functions [4–6]. Perinatal stresses causing immune activation such as maternal infection, during pregnancy, has also been implicated with pathogenesis of ASD and other neuropsychiatric conditions [7,8]. In our clinic, we have also observed that some ASD subjects reveal fluctuating behavioral symptoms and even repeated loss of once acquired cognitive activity following

each immune insult. In addition, our previous research indicated an association between fluctuating behavioral symptoms and dysregulated innate immune responses [9,10].

The MIA model relies on the induction of sterile immune activation through non-specific stimulants of innate immunity. However, the lack of antigen (Ag)-specific immune memory in innate immunity, as opposed to adaptive immunity, has made the lasting effects of puzzling. Although neuronal damage caused by innate immune activation in mothers could cause lasting effects in neuronal development, if this is the case, such neurodevelopmental changes will remain static, as seen in patients with cerebral palsy. However, we observed that ASD subjects reveal fluctuating behavioral symptoms following immune insults (mainly microbial infection) repeatedly affect behavioral symptoms along with a high frequency of co-morbid conditions. That is, these ASD subjects develop multi-system inflammation not restricted to the brain (inflammatory subtype) [11]. In addition, direct damage to the brain by MIA would not explain MIA's lasting effects on T cell functions [6].

The recent discovery of innate immune memory (IIM) may explain many aspects of the lasting effects of MIA. Reprogramming of innate immunity via metabolic and epigenetic changes following initial innate immune stimulus, was first realized as a part of the non-specific effects of vaccinations such as Bacillus Calmette–Guérin (BCG) [12]. Such effects of IIM were termed as trained immunity (TI). It has been shown that TI causes hyper-responsiveness of innate immunity against second stimuli, which can be totally different from the 1st stimuli, since TI is established through metabolic and subsequent epigenetic modulation of innate immune cells [12–15]. On the other hand, the 1st innate immune stimuli can also cause subsequent hypo-responsiveness of innate immunity (tolerance) as shown in endotoxin tolerance [16]. Innate immune tolerance is also thought to be induced by metabolic and epigenetic regulation [17,18]. Induction of innate immune memory appears to be associated with dose, kinds of stimuli, and other environmental and genetic factors [19].

We have previously reported that some ASD subjects reveal changes in monocyte activation status at the level of cytokine production, expression of microRNA (miRNAs), and mitochondrial respiration [20,21]. Our research has also revealed that the above-described changes are closely associated with production of monocyte cytokines as well as the ratio of IL-1β/IL-10 [20,21]. In these studies, we used a panel of agonists of toll like receptors (TLRs) for stimulating peripheral blood monocytes (PBMo), main innate immune cells in the peripheral blood. It was found that subgrouping of ASD subjects based on IL-1β/IL-10 ratios help identify ASD subjects who have dysregulated innate immune responses; high or low ratio ASD subgroup reveal marked differences in miRNA expression and mitochondrial respiration [20,21]. These findings support the presence of maladapted TI in these ASD children.

In the animal model of TI, with β-glucan derived from *Candida albicans* used as an initial stimulus, the IL-1β family was reported to have a major role in the development of TI [22]. Additionally, β-glucan tends to induce TI, but not tolerance, while endotoxin tends to induce tolerance when given at a low dose [23]. β-glucan is a representative dectin-1 agonist, which stimulates canonical inflammasomes of innate immune cells, leading to cleavage of proinflammatory cytokines, especially IL-1β and IL-18. Inflammasome mediated pathways are different from TLR mediated pathways, but often they exert converging effects on activation of innate immunity [24,25]. It is also noted that miRNAs have important regulatory roles in inflammasome priming and activation [25]. Thus, it can be hypothesized that TI suspected in ASD children will include multiple pathways of innate immunity, including both TLRs and inflammasome mediated pathways.

We also noted that in the low IL-1β/IL-10 ratio ASD subgroup, there is a higher frequency of non-IgE mediated food allergy (NFA) as compared to the normal ratio ASD subgroups [21]. NFA is thought to be associated with cellular immune reactivity to common food proteins, and we have previously reported that ASD children reveal increased T cell cytokine production (TNF-α, IFN-γ, and IL-12) against cow's milk proteins at a high frequency. Given these findings, it may be hypothesized that exaggerated innate immune responses due to maladapted TI could result in the alternation of

adaptive immune responses, resulting in unwanted T cell responses to benign environmental Ags, such as food proteins, and self-Ags.

This study was formulated to examine the above described hypothesis by assessing whether changes in the monocyte cytokine production are detected under β-glucan stimulated cultures in IL-1β/IL-10 based ASD subgroups, and whether T cell cytokine production in response to luminal antigens differ among IL-1β/IL-10 based ASD subgroups. This study determined monocyte cytokine profiles in responses to β glucan, in comparison with their responses to TLR agonists, and T cell cytokine profiles in response to food and candida proteins in ASD (*n* = 152) and control non-ASD (*n* = 41) subjects. Samples from ASD subjects were subdivided into the IL-1β/IL-10 based subgroups (high, normal, and low) as reported previously [20]. The results obtained support our hypotheses described above.

#### **2. Results**

#### *2.1. Clinical Characteristics in the IL-1*β/*IL-10 Based ASD Subgroups*

A summary of clinical characteristics of the ASD subgroups are shown in Table 1. There are no statistical differences in gender frequency and clinical characteristics (ASD severity, cognitive development, sleep, GI symptoms, history of non-IgE mediated food allergy (NFA), seizure disorders, specific antibody deficiency, AR, and asthma) by Chi-Square test (*p* > 0.05) (Table 1). Ages in the ASD subgroups did not differ (*p* > 0.05 by Kruskal–Wallis test). History of NFA tended to be higher in the low IL-1β/IL-10 ratio group than the normal ratio ASD subgroup, but this was not statistically significant (*p* = 0.1286 by Fisher's exact test).

**Table 1.** Demographics and clinical features in the IL-1β/IL-10 based autism spectrum disorder (ASD) subgroups.

1 Abbreviations used: AA, African Americans; ADHD, attention deficit hyperactivity disorder; AEDs, anti-epileptic drugs; F, female; GI, gastrointestinal; M, male; NFA, non-IgE mediated food allergy; SD, standard deviations; SSRIs, selective serotonin receptor inhibitors; W, Caucasians. 2 No significant differences in frequency of clinical characteristics in the IL-1β/IL-10 ratio based ASD subgroup by Chi-Square test (*p* > 0.05).

We also assessed the frequency of prescription medication use in the IL-1β/IL-10 ratio based ASD subgroups. Table 2 summarizes frequency of intake of prescription medications at the time of sample obtainment. For ASD subjects who had samples taken at multiple time points, the medications taken at each time-point was assessed. There was no difference in the frequency of medication intake between the ASD subgroups (*p* > 0.05 by Chi-Square test). Co-variance analysis rejected the effects of the above described clinical characteristics and prescription medication use on IL-1β/IL-10 ratios under all the culture conditions tested in this study.

**Table 2.** Frequency of medication use in the IL-1β/IL-10 based ASD subgroups.

1 Frequency of medication use at the time of sample obtainment. It should be noted that samples were obtained at 2–4 time points in a total of 25 ASD subjects. Frequency of medication use did not significantly differ in the IL-1β/IL-10 based ASD subgroups by Chi-Square test (*p* > 0.05). 2 Medications on which study subjects were on when blood samples were obtained.

#### *2.2. Changes in Monocyte Cytokine Production in Response to Candida Heat Extract (*β*-Glucan) the IL-1*β/*IL-10 Based ASD Subgroups*

Since TLR agonists were used to generate IL-1β/IL-10 ratios for subgrouping ASD samples previously, we also assessed the effects of β-glucan, a ligand to dectin 1 which activates innate immunity through inflammasome. We used candida heat extract as a source of β-glucan. One of the reasons that we used β-glucan is that this is widely used for TI induction [12]. Under β-glucan stimulated cultures, significant changes in IL-1β/IL-10 ratios and in the production of IL-1β, IL-10, TNF-α, sTNFRII, and CC-chemokine-ligand-2 (CCL2) were observed (Table 3). There is a tendency of high production of inflammatory cytokines (IL-1β and TNF-α) in the high ratio ASD subgroup, while the low ratio ASD subgroup revealed higher levels of counter-regulatory cytokines (IL-10, sTNFRII, and TGFβ) and CCL2, chemokine under the β-glucan stimulated cultures (Table 3). Therefore, differences were the most apparent between the high vs. low ratio ASD subgroups (Table 3). Such changes were not observed in the production of IL-6, IL-12, IL-23, and CCL7. The changes in monocyte cytokine profiles were more evident with stimulus of β-glucan as compared to monocyte cytokine production without stimuli (Table 3).

We were able to obtain samples at 2–4 different time points in 22 ASD subjects to test longitudinal changes in each subject. With the use of repeated measures of analysis of variance, significant variability was observed in IL-1β/IL-10 ratios under the cultures without stimulus and β-glucan stimulated cultures (F-ratio 5.0 and 3.88, with *p* < 0.01 and *p* < 0.02, respectively).

**Table 3.** Changes in IL-1β/IL-10 ratios and monocyte production in the IL-1β/IL-10 based ASD subgroups under β-glucan stimulated culture conditions.

1 Total sample numbers in the group; 2 The results are expressed as mean values ± SD. 3 Cytokine concentrations are shown as pg/mL; 4 *p* values by one-way ANOVA with α = 0.2 (Terry–Hoeffding test). Comparisons between groups by ANOVA with α = 0.05. 5 Values from non-ASD controls differed from ASD subgroups. 6 Differed from other ASD subgroups by ANOVA (α = 0.05).

#### *2.3. Changes in T Cell Cytokine Production in the IL-1*β/*IL-10 Based ASD Subgroups*

We then tested whether T cell cytokine production differs among the study groups, including non-ASD controls. Numerically significant differences were found between the study groups in the production of IFN-γ (medium only and in response to candida), TNF-α (in all the culture conditions), and IL-10 (in response to β-LG and gliadin), as shown in Table 4. When the whole ASD samples were compared with non-ASD controls, no significant differences were observed (*p* > 0.05). These changes are mainly due to higher production of TNF-α and IL-10 in the high IL-1β/IL-10 ratio ASD subgroup as compared to other study groups (Table 4). On the other hand, IFN-γ production was lower in the high ratio ASD subgroup, as compared to the normal and low ratio ASD subgroups in the cultures without stimuli or with candida protein (Table 4). In addition, the low ratio subgroup revealed lower IL-10 production in response to gliadin than other study groups (Table 4).

When we assessed the variability of the production T cell cytokines, no significant time-dependent variations were found under all culture conditions, except for IFN-γ production under soy-protein stimulated cultures (F-ratio 4.93, *p* < 0.01).

**Table 4.** T cell cytokine production by peripheral blood mononuclear cells (PBMCs) in response to β-LG, gliadin, or candida Ag in the IL-1β/IL-10 based ASD subgroups.

1 Total sample numbers in each study group. 2 The results are expressed as mean values±SD. Cytokine concentrations are shown as pg/mL; 3 *p* values by one-way ANOVA with α = 0.2 (Terry–Hoeffding test). Comparisons between groups by ANOVA with α = 0.05. 4 Values from non-ASD controls differed from ASD subgroups. 5 Differed from other ASD subgroups by ANOVA (α = 0.05).

#### *2.4. Associations between T Cell Cytokine Production by Peripheral Blood Mononuclear Cells (PBMCs) and Monocyte Cytokine Production*

Since we observed differences in the production of IFN-γ, TNF-α, and IL-10 by peripheral blood mononuclear cells (PBMCs) under cultures stimulated with β-LG, gliadin, and candida, we assessed whether production of these cytokines revealed any correlation with the cytokines produced by PBMo in response to TLR agonists and/or β-glucan.

When IFN-γ production by PBMCs were compared with monocyte cytokine profiles, we observed significant correlations between IFN-γ levels produced by PBMCs and levels of monocyte cytokines produced by PBMo in the ASD subgroups (Table 5). Interestingly, correlations markedly differed among the ASD subgroups, while the non-ASD controls revealed little correlations. Such correlations were most notable between IL-6 and CCL2 production by PBMo and IFN-γ production by PBMCs in response to candida proteins in the normal ratio ASD subgroup (Table 5).

**Table 5.** Correlations between IFN-γ production by PBMCs and monocyte cytokine production.

**Table 5.** *Cont.*

1 Total sample numbers in each study group; 2 Correlation coefficient by Spearman test with *p* value shown in ( ) when *p* value is at least *p* < 0.05. 3 Culture conditions monocyte cytokines produced under.

As for TNF-α production, production of TNF-α by PBMCs in response to β-LG and gliadin and that by PBMo revealed a significant correlation with TNF-α production by PBMo in response to TLR agonists and β-glucan (Table 6). Such correlations were most notable in the high ratio ASD subgroup. However, such correlations were also noted in the low ratio and non-ASD subgroups as well, but to a lesser extent (Table 6). On the other hand, few correlations were found in the normal ratio ASD subgroups (Table 6). It should be noted that PBMCs from the normal ratio ASD subgroup produced equivalent amounts of TNF-α, as compared to both the low ratio ASD subgroup and the non-ASD controls (Table 4), making it unlikely that the lack of correlations in the normal ASD subgroup is associated with low TNF-α production. The high and low ratio ASD subgroups revealed positive associations between TNF-α production by PBMCs in response to candida protein, and production of IL-1β, IL-6, IL-10, and sTNFRII by PBMo.

**Table 6.** Correlations between TNF-α production by PBMCs and monocyte cytokine production.

**Table 6.** *Cont.*

1 Total sample numbers in each study group; 2 Correlation coefficient by Spearman test with *p* value shown in ( ) when *p* value is at least *p* < 0.05. 3 Culture conditions monocyte cytokines produced under.

When we examined correlations between IL-10 production by PBMCs and monocyte cytokine production, we found different correlating patterns between the ASD subgroups in response to β-LG and gliadin, while the non-ASD controls revealed little correlations (Table 7). In response to candida protein, only the high ratio ASD subgroup revealed positive correlations between IL-10 production by PBMCs and monocyte cytokine production (Table 7). Non-ASD controls revealed positive correlations between IL-10 production by PBMCs and sTNFRII production by PBMo, and negative correlations between IL-10 production by PBMCs and CCL2 production by PBMo.

**Table 7.** Correlation between IL-10 production by PBMCs and monocyte cytokine production.

**Table 7.** *Cont.*

1 Total sample numbers in each study group; 2 Correlation coefficient by Spearman test with *p* value shown in ( ) when *p* value is at least *p* < 0.05. 3 Culture conditions monocyte cytokines produced under.

We also observed significant correlations between IL-12 production by PBMCs and monocyte cytokine production, although IL-12 production did not differ between IL-1β/IL-10 based ASD subgroups (Table 8). Namely, the low ratio ASD subgroup revealed significant correlations between PBMC IL-12 production and PBMo cytokine production, while the other ASD subgroups and non-ASD controls revealed little correlations between these two parameters.

**Table 8.** Correlation between IL-12 production by PBMCs and monocyte cytokine production.

**Table 8.** *Cont.*

1 Total sample numbers in each study group; 2 Correlation coefficient by Spearman test with *p* value shown in ( ) when *p* value is at least *p* < 0.05. 3 Culture conditions monocyte cytokines produced under.

#### **3. Discussion**

Changes in IL-1β/IL-10 ratios and monocyte cytokine profiles in the ASD subgroups, which we have reported before with the use of TLR agonists, were also found under the β-glucan stimulated cultures in this study. This finding indicates that the innate immune abnormalities observed in some ASD children involve both TLR and inflammasome mediated pathways. T cell cytokine profiles (IFN-γ, TNF-α, and IL-10) also differed among the IL-1β/IL-10 ratio based ASD subgroups under cultures

stimulated with food proteins (β-LG and gliadin) and candida proteins. Moreover, T cell cytokine production revealed close correlations with monocyte cytokine production, specific to each ASD subgroup. These findings support our hypothesis of sustained effects of maladapted TI on adaptive immunity in some ASD children.

The immune system operates with two components, innate and adaptive immunity. Innate immunity exerts rapid Ag-non-specific immune responses to contain infection until Ag-specific adaptive immunity sets in to clear hazards from the human body [26–28]. As opposed to adaptive immunity, innate immunity has been thought to lack lasting effects, since innate immune responses are Ag non-specific. However, for the past few years, it became clear that innate immunity can generate lasting effects or IIM by inducing metabolic and epigenetic changes in innate immune cells [29,30]. The key difference of IIM from adaptive immune memory is that the 2nd stimuli to provoke IIM can be totally unrelated to the 1st stimuli.

The presence of IIM was first convincingly shown as non-specific effects of BCG. Namely, BCG vaccination led to a greater reduction of infant mortality than predicted for decrease in mortality caused by tuberculosis in epidemiological and then randomized prospective studies [31,32]. These findings were further confirmed in rodent IIM models: BCG vaccination reduced from *Candida albicans* in the absence of lymphocyte that generates adaptive immune responses [33]. TI is the term that was given to describe these non-specific effects of IMM. Further studies revealed that TI can be generated in various innate immune cells including monocyte-macrophage lineage cells, natural killer (NK) cells [34], and bone marrow progenitor cells of innate immunity [35,36].

Metabolic and epigenetic changes have been extensively studied in 'in vitro' models of TI with the use of β-glucan derived from *Candida albicans*. Others have shown that β-glucan stimulates inflammasome through dectin-1, resulting in activation of signaling pathways involving Akt/PTEN/mTOR/HIF-1α [37]. This will shift oxidative phosphorylation (ATP synthesis) to glycolysis, resulting in decreased basal cellular respiration, and increased consumption of glucose, thereby leading to increase in lactate production [37]. These metabolic changes also affect synthesis of cholesterol and phospholipids [13,38], leading the replenishment of the Krebs cycle through formation of glutamate and α-keto-glutamate, and finally accumulation of fumarate [13,38]. A high concentration of fumarate hinders enzymatic actions of H3K4 demethylase, eventually resulting in epigenetic reprogramming [38–40].

Pre-administration of proinflammatory innate cytokines (IL-1, TNF-α, and IL-6) provided protection against a variety of microbes [41]. Among the cytokines administered, IL-1 showed superior effects over TNF-α or IL-6 [41]. In humans, TI associated protection has been mainly implicated with IL-1β and other IL-1 families [22]. Given the role of IL-1β in TI, excessive, dysregulated production of IL-1β is likely to cause maladapted TI, and resultant pathogenic consequences. Patients with gene mutations that lead to the over-production of IL-1β are known to reveal inflammatory symptoms involving multiple organs including the brain [22,42]. Moreover, chronic inflammatory conditions, including neuropsychiatric conditions have been implicated with maladapted TI [29,43].

As opposed to TI, tolerance is another form of IIM, causing lasting hypo-responsiveness against stimuli of innate immunity which is bested studied in endotoxin tolerance [44]. Upon induction of endotoxin tolerance, innate immune cells were reported to produce less inflammatory cytokines (TNF-α, IL-12, IL-6), but generate more counter-regulatory cytokines (IL-10 and TGF-β) with the subsequent endotoxin challenge [45,46]. LPS, a representative endotoxin, exerts its actions through TLR4, that results in activation of downstream signaling pathways the myeloid differentiation factor 88 (MyD88) and the TIR-domain-containing adaptor-inducing interferon-β (TRIF) [17]. Under the state of endotoxin tolerance, suppression of activation of these pathways are shown to be accomplished through multiple stages; down-regulation of TLR4, suppressed recruitment of MyD88 and TRIF to TLR4, reduced activation of IL-1 receptor-associated kinase (IRAK)1 and IRAK4, diminished activity of nuclear factor κ chain of B cell (NFκB), and up-regulated expression of counter-regulatory molecules such as SHIP1 (SH2 domain-containing inositol phosphatase 1) [47]. Induction of innate immune tolerance appears crucial in maintaining brain homeostasis; impaired innate immune tolerance has been implicated with pathogenesis of chronic neurodegenerative disorders including Alzheimer's disease [29].

We have been analyzing monocyte functions in ASD subjects for the past 4–5 years, and we have reported variable differences in monocyte cytokine profiles in ASD subjects, which have been best reflected in IL-1β/IL-10 ratios [20]. This finding led us to subgroup ASD PBMo samples on the basis of IL-1β/IL-10 ratios generated by stimulating PBMo with a panel of TLR agonists [20]. The IL-1β/IL-10 ratio based subgrouping revealed notable changes in miRNA expression and mitochondrial respiration across the IL-1β/IL-10 based ASD subgroups [20]. In addition, we also found the IL-1β/IL-10 ratio based ASD subgroup specific associations between metabolic changes (mitochondrial respiration) and monocyte cytokine profiles [21]. This was also true when we examined correlations between serum miRNA levels and monocyte cytokine profiles in the IL-1β/IL-10 ratio based ASD subgroups (manuscript submitted for publication). In addition, it is our observation that the ASD subjects whose PBMo revealed low or high IL-1β/IL-10 ratios tend to reveal fluctuating behavioral symptoms following immune insults.

The above described findings indicate that the innate immune responses in the high and low IL-1β/IL-10 ratio ASD subgroups may be altered as a result of maladapted TI, leading to clinical pictures of chronic inflammation affecting multi-organs. To address this possibility, we assessed whether responses to β-glucan also differed in the IL-1β/IL-10 ratio based ASD subgroups. Our results revealed significant differences of IL-1β/IL-10 ratios and monocyte cytokine profiles across the IL-1β/IL-10 based ASD subgroups under β-glucan stimulated cultures (Table 5). If the previously observed changes of monocyte cytokine profile were only limited to TLR pathways, such changes would be more likely to be associated with defects specific to TLR mediated signaling pathways. β-glucan stimulates a connonical inflammasome pathway and amplifies signaling through the TLR pathway [25]. Thus, the fact that we observed similar changes under β-glucan stimulated cultures supports the role of maladapted TI in ASD subjects; TI has been reported in multiple signaling pathways associated with innate immunity [30].

This study also addressed how T cell cytokine profiles in response to recall Ags are associated with changes in innate immune responses in the IL-1β/IL-10 based ASD subgroups. TI induced excessive innate immune responses are likely to induce more exuberant adaptive responses, thereby increasing the risk of unwanted immune reactivity to benign environmental antigens such as food antigens. However, if chronic activation of adaptive immunity through maladapted TI is sustained, this may eventually lead to suppression of adaptive immunity, regaining immune homeostasis. Such secondary suppression of adaptive immunity could occur by clonal deletion or through actions exerted by regulatory T cells [48,49]. However, how dysregulated innate immune abnormalities (presumably maladapted TI) affect adaptive immune responses in ASD children is poorly understood. Since ASD subjects have a high frequency of GI symptoms and delayed type reactivity to common food proteins and candida proteins [50], we opted to check T cell cytokine profiles against representative food proteins and candida antigens, when assessing associations between T cell cytokine profiles and monocyte cytokine profiles.

With this analysis, we observed significant differences in the production of TNF-α and IL-10 in response to β-LG and gliadin between the IL-1β/IL-10 based ASD subgroups. Namely, the high ratio groups tended to reveal higher levels of these two cytokines than the other subgroups (Table 5). On the other hand, the normal and low ratio ASD subgroups revealed higher IFN-γ production in response to candida protein (Table 5). As for TNF-α in response to candida Ag, both the high and low IL-1β/IL-10 ratio subgroups revealed a higher TNF-α production (Table 5). In contrast, IL-10 production against candida Ag was lower in the ASD subgroups, with the normal ratio group being the lowest. These results indicate complex associations between adaptive immune responses and altered innate immune responses in the ASD subjects.

When correlations between T cell cytokine production against β-LG and gliadin and monocyte cytokine profiles were assessed, positive correlations between IFN-γ production and production of monocyte cytokines (IL-1β, IL-10, and IL-6) were predominantly observed in the normal ratio ASD subgroup (Table 6). In contrast, correlations between production of TNF-α, IL-10, and IL-12 by PBMCs in responses to these food proteins and monocyte cytokines were predominantly found in either the high ratio ASD subgroup (for TNF-α and IL-10), or the low ratio ASD subgroup (for IL-12) (Tables 7–9). Non-ASD controls revealed little correlations except for TNF-α production by PBMCs and TNF-α production by PBMo (Table 7). These results indicate that responses to benign environmental antigens, such as food proteins, appear to be differently affected by innate immune responses in the IL-1β/IL-10 based ASD subgroups, perhaps reflecting maladapted TI.

**Table 9.** Demographics of ASD subjects and non-ASD controls.

1 Abbreviations used; AA; African American, ASD; autism spectrum disorder, F; female, GI; gastrointestinal, M; male, NFA: non-IgE mediated food allergy, SD; standard deviation, W; Caucasian. 2 Age at the time of study enrollment is shown.

When correlations between T cell cytokine production against candida proteins and monocyte cytokine profiles were assessed, we observed similar results; predominant correlations with monocyte cytokines in the normal ratio group for IFN-γ production, with more predominant correlations in the high and/or low ratio ASD subgroups for TNF-α, IL-10, and IL-12 production by PBMCs (Tables 6–9). Non-ASD controls only revealed positive correlations between IL-10 production by PBMCs and monocyte cytokine productions (Table 8). Since we were testing cellular immune reactivity to benign environmental Ags (food proteins and candida which is a normal component of the gut microbiome), a lack of close associations of inflammatory T cell cytokine profiles against these benign Ags in the non-ASD controls, and monocyte cytokine profiles may indicate the establishment of immune tolerance in adaptive immunity. Negative correlations between IL-10 production against candida protein and monocyte cytokine production may also indicate the establishment of immune homeostasis in non-ASD controls (Table 8). While the presence of close associations between monocyte cytokine profiles and T cell cytokine profiles (especially inflammatory cytokines—IFN-γ, TNF-α, and IL-12) in the ASD subgroups may indicate on-going effects of innate immunity in the absence of oral tolerance; this may also reflect maladapted TI.

When we assessed the clinical features in the IL-1β/IL-10 based ASD subgroups, we observed a high frequency of GI symptoms and a history of NFA with favorable responses to dietary interventions (usually gluten-free, dairy-free diet), consistent with our previous reports [20,21]. The IL-1β/IL-10 low ratio subgroup tended to have a higher frequency of history of NFA. β-GL is thought to be a major milk protein component associated with delayed type cellular immune reactivity in predisposed individuals. Interestingly, in the low IL-1β/IL-10 ratio ASD subgroup, IL-12 production by PBMCs in response to β-GL, a component was negatively correlated with counter-regulatory monocyte cytokines (IL-10, sTNFRII, and TGF-β) (Table 9). In contrast, IL-10 production in response to the β-GL was mostly positively correlated with monocyte cytokines production in the low ratio ASD subgroup (Table 9). This may indicate that altered innate immune responses may be regulating T cell cytokine production in the direction of suppressing excessive responses to recall antigens. On the other hand, in the high ratio group, there are positive correlations between IFN-γ, TNF-α, and IL-10 production in response to β-GL and monocyte cytokines (both inflammatory and counter-regulatory cytokines). In the high ratio

group, such negative regulations may not be properly in place. Alternatively, with immune stimuli, such negative regulations may be easily lost secondary to maladapted TI.

#### **4. Materials and Methods**

#### *4.1. Study Subjects*

This study protocol was approved as protocols 15:45 (approval on 28 April 2016) and 17:53 (approval on 26 April 2018) by the Institutional Review Board, Saint Peter's University Hospital, New Brunswick, NJ, USA. This study included both ASD and typically developing (TD), non-ASD control subjects. We obtained the signed consent forms prior to sample obtainment from parents or guardians when the study subjects were minor (<18 years of age) or were judged unable to give consent by him/herself due to intellectual disability. History of comorbidities including food allergy (FA), asthma, allergic rhinitis (AR), specific antibody deficiency (SAD), sleep disorder, and seizure disorders were assessed in ASD children by history taking and medical chart review. This study excluded subjects with chromosomal abnormalities, well defined gene mutations, or chronic diseases involving major organs, but did not exclude subjects with minor medical conditions highly prevalent in general population, such as seasonal allergy. In this study, both ASD and non-ASD, TD subjects were enrolled, and the signed consent forms were obtained prior to entering the study. Consent was obtained from parents/guardians if participant was a minor (<18 years old) or parents/guardians had custody. For ASD children, we also assessed whether they had a history of FA, asthma, AR, SAD, or seizure disorders. Subjects diagnosed with chromosomal abnormalities, other genetic diseases, or well characterized chronic medical conditions involving a major organ, were excluded from the study. Subjects with common minor medical conditions such as AR, mild to moderate asthma, eczema were not excluded from the study.

ASD subjects: ASD subjects (*n* = 152) were recruited from the Pediatric Allergy/Immunology Clinic. Diagnosis of ASD was made at various autism diagnostic centers, including ours. The ASD diagnosis was based on the Autism Diagnostic Observation Scale (ADOS) and/or Autism Diagnostic Interview-Revisited (ADI-R), and other standard measures. ASD subjects were also evaluated for their behavioral symptoms and sleep habits with the Aberrant Behavior Checklist (ABC) [51] and the Children's Sleep Habits Questionnaires (CSHQ) [52], respectively. Information regarding cognitive ability and adaptive skills were obtained from previous school evaluation records performed within 1 year of enrollment in the study; these results were based on standard measures such as the Woodcock-Johnson III test (for cognitive ability), and Vineland Adaptive Behavior Scale (VABS) (for adaptive skills) [53].

Non-ASD controls: A total of 41 non-ASD subjects served as controls. These subjects were recruited in the Pediatrics Subspecialty and General Pediatrics Clinics at our institution. These subjects were typically growing and satisfied our exclusion criteria.

Table 1 reveals demographics of study subjects. Gender difference did not reveal significant changes regarding monocyte cytokine profiles and IL-1β/IL-10 ratios in both ASD and non-ASD groups, as reported before [21].

Diagnosis of FA: IgE mediated FA was diagnosed with reactions to offending food, by affecting the skin, GI, and/or respiratory tract immediately (within 2 h) after intake with positive prick skin testing (PST) reactivity, and/or presence of food allergen-specific serum IgE. Non IgE mediated FA (NFA) was diagnosed if GI symptoms resolved, following implementation of a restricted diet (i.e., avoidance of offending food), and symptoms recurred with re-exposure to offending food [54]. NFA was also defined as being non-reactive to PFT and negative for serum IgE specific for food allergens [54].

Diagnosis of asthma and AR: AR and allergic conjunctivitis (AC) were diagnosed when subjects had corresponding clinical features along with positive PST reactivity and/or positive serum IgE specific to [55,56]. Asthma was diagnosed following the asthma guidelines from the Expert Panel Report 3 [57].

Antibody deficiency syndrome: When the subject revealed protective levels of antibodies in less than 11 of 14 serotypes of *Streptococcus pneumonia* after the booster dose of Pneumovax® or PCV13®, he/she was diagnosed with SAD [58]. Antibody (Ab) levels greater than 1.3 μg/mL were considered protective [58].

#### *4.2. Sample Collection*

Venous blood samples were obtained by physician in this study. We obtained one sample from each non-ASD control. As for ASD subjects, we obtained multiple blood samples from select ASD subjects (*n* = 22), in order to assess variability of T cell cytokine profiles. If parents or study subjects preferred, we applied a topical lidocaine/prilocaine cream (Emla cream®,IGI Laboratories, Buena, NJ, USA) to the site of venipuncture prior to blood sampling.

#### *4.3. Cell Cultures*

Ficoll–Hypaque density gradient centrifugation was used for separating PBMCs. From PBMCs, PBMo were further purified using magnetic beads labeled with anti-CD3, CD7, CD16, CD19, CD56, CD123, and glycophorin A (monocyte separation kit II—human, MILTENYI BIOTEC, Cambridge, MA, USA). Namely, this column depletes T, B, natural killer, and dendritic cells from PBMCs with combination of these antibodies.

Cytokine production by purified PBMo was induced by incubating cells overnight (2.5 × 105 cells/mL) with a panel of agonist of toll like receptors (TLRs). This assay system was designed to reflect the effects of microbial byproducts commonly encountered in real life. Lipopolysaccharide (LPS), a TLR4 agonist, represents a signaling pathway activated in response to a Gram negative (G (–)) bacteria. Zymosan, a TLR2/6 agonist, mimics an innate activation signal in response to G (+) bacteria and fungi. CL097, a TLR7/8 agonist, activates innate signaling pathways in response to ssRNA viruses that cause common respiratory infection. We adapted these stimuli in our assay system, since these innate immune stimuli have been widely used by others. In addition, candida heat extract as a source of β-glucan, dectin-1 agonist, was used as well as a representative C-lectin receptor agonist. PBMos were incubated overnight with LPS (0.1 μg/mL, GIBCO-BRL, Gaithersburg, MD, USA), zymosan (50 μg/mL, Sigma-Aldrich, St. Luis, Mo, USA), C097 (water-soluble derivative of imidazoquinoline, 20 μM, InvivoGen, San Diego, CA, USA), and candida heat extract (HCKA, heat killed candida albicans (107 cells/mL, InVivogen, San Diego, CA) as a source of β-glucan, a dectin 1 agonist, in RPMI 1640 with additives as previously described [59]. Overnight incubation (16–20 h) was adequate to induce the optimal responses in this setting. The culture supernatant was used for cytokine assays.

Production of T cell cytokines (IFN-γ, IL-5, TNF-α, IL-10, IL-12p40, IL-17, and TGF-β) was assessed by incubating PBMCs (106 cells/mL) with representative recall antigens including β-lactoglobulin (β-LG, 10 μg/mL, Sigma Aldrich), soy protein (5 μg/mL, Ross, Nutley, NJ, USA), gliadin (10 μg/mL, Sigma-Aldrich), milk protein (100 μg/mL, Ross), and candida protein (5 μg/mL, Greer, Lenoir, NC, USA) for 4 days in RPMI1540 with additives as reported previously [60]. As noted in our previous study, a 4 day incubation period resulted in the optimal production of these cytokines in this culture setting.

Levels of C-C chemokine ligand 2 (CCL2), CCL7, interferon-γ (IFN-γ), IL-1β, IL-5, IL-6, IL-10, IL-12p40, IL-17, transforming growth factor-β (TGF-β), tumor necrosis factor-α (TNF-α), and soluble TNF receptor II (sTNFRII) cytokines were measured by enzyme-linked immuno-sorbent assay (ELISA); 10–100 μL/well supernatants were used for ELISA. The OptEIA™ Reagent Sets (BD Biosciences, San Jose, CA, USA) were used for ELISA of IFN-γ, IL-1β, IL-6, IL-10, IL-12p40, and TNF-α. For CCL2, CCL7, IL-17 (IL-17A), sTNFRII, and TGF-β ELISA, reagents were obtained from BD Biosciences and R & D (Minneapolis, MN, USA). IL-23 ELISA kit was purchased from eBiosciences, San Diego, CA. Intraand inter-variations of cytokine levels were less than 5%.

#### *4.4. Categorizing ASD Samples on the Basis of IL-1*β/*IL-10 ratios*

Previously, we reported that changes in the IL-1β/IL-10 ratios best reflect altered cytokine profiles and miRNA expression by PBMo [20]. We divided ASD samples into subgroups based on the IL-1β/IL-10 ratios produced by ASD PBMo, following the criteria used in our previous study [20], as outlined below. In this study, we used IL-1β/IL-10 ratios generated under cultures with medium only, LPS, zymosan, CL097, and candida heat extract as a source of β-glucan.

High IL-1β/IL-10 ratio subgroup: In this subgroup, ASD PBMo showed 2 standard deviation (SD) higher IL-1β/IL-10 ratios under at least one culture condition and/or 1 SD greater IL-1β/IL-10 ratios under more than two culture conditions, as compared to control PBMo.

Normal IL-1β/IL-10 ratio subgroup: ASD PBMo revealed IL-1β/IL-10 ratios less than +1 SD, and greater than –1 SD under all the culture conditions, or IL-1β/IL-10 ratios greater than 1 SD but less than 2 SD under only one culture condition, as compared to control PBMo.

Low IL-1β/IL-10 ratio subgroup: ASD PBMo revealed IL-1β/IL-10 ratios less than –1 SD under at least one culture condition, as compared to control PBMo.

Among 25 ASD subjects in whom we obtained blood samples at 2–4 time points, most subjects revealed that blood samples from these subjects were categorized in the same group except for six ASD subjects. PBMo from these six ASD subjects showed high ratios at 1–2 time points and normal/low ratios at one time point. Eleven ASD subjects revealed low ratios at 1–2 time points and normal ratios at one time point. When comparing their clinical characteristics, they were categorized as high and low IL-1β/IL-10 ratio groups, respectively.

#### *4.5. Statistical Analysis*

We used a two tailed Mann–Whitney test for comparison of two sets of numerical data. Comparison of multiple data sets was assessed by one-way ANOVA and/or Kruskal—Wallis test. Normality of the numerical data were assessed by skewness and kurtosis (Omnibus) with α = 0.2. When assessing differences in frequency between two groups, we used the Fisher exact test. For assessing differences in frequency among multiple groups, we used the Chi-Square test and the Likelihood ratio. A linear association between two data sets was determined by Spearman test. A *p* value of less than 0.05 was considered nominally significant. Co-variance analysis was done with the use of general linear model for a fixed factor or for a variable factor. For assessing longitudinal changes, repeated measures of analysis of variance were used. NCSS19 (NCSS, LLC. Kaysville, UT, USA) was used for statistical analysis.

#### *4.6. Availability of Data and Material*

Clinical features of the ASD are available through NDAR data base (https://ndar.nih.gov/). The additional datasets used and/or analyzed during the current study are available from the corresponding author on reasonable request.

#### **5. Conclusions**

This study revealed similar changes in IL-1β/IL-10 ratios and monocyte cytokine profiles under β-glucan stimulated cultures as observed in the cultures stimulated with TLR agonist, indicating that changes in innate immune responses are not limited to TLR pathways. Correlations between the changes in monocyte cytokine profiles and T cell cytokine profiles in response to benign environmental Ags differed across the IL-1β/IL-10 based ASD subgroups. Such correlations were much less evident in non-ASD controls, supporting possible effects of maladapted TI in some ASD children. Further analysis of the status of innate immunity in association with TI or tolerance will be helpful in understanding the inflammatory subtype of ASD.

**Author Contributions:** H.J. was responsible for the study design, recruitment of the study subjects, collection of clinical information and blood samples, analysis of the overall data, and preparation of most of this manuscript. L.G. conducted cytokine production assays with the use of purified monocytes and conducted T cell cytokine production assay with the use of recall Ags.

**Funding:** This study was supported by funding from Autism Research Institute, San Diego, CA, Jonty Foundation, St. Paul, MN, and the Governor's Council for Medical Research and Treatment of Autism (CAUT16APL007), DHHS, Trenton, NJ.

**Acknowledgments:** We are thankful for L. Huguenin, Ph.D. for critically reviewing this manuscript.

**Conflicts of Interest:** There is no conflict of interest to declare. The funders had no role in the design of the study; in the collection, analyses, or interpretation of data; in the writing of the manuscript, or in the decision to publish the results.

#### **Abbreviations**

#### **References**

© 2019 by the authors. Licensee MDPI, Basel, Switzerland. This article is an open access article distributed under the terms and conditions of the Creative Commons Attribution (CC BY) license (http://creativecommons.org/licenses/by/4.0/).

*Article*

## **Gender Related Changes in Gene Expression Induced by Valproic Acid in A Mouse Model of Autism and the Correction by S-adenosyl Methionine. Does It Explain the Gender Di**ff**erences in Autistic Like Behavior?**

**Liza Weinstein-Fudim 1, Zivanit Ergaz 1, Gadi Turgeman 2, Joseph Yanai 1, Moshe Szyf 3 and Asher Ornoy 1,\***

Received: 27 September 2019; Accepted: 22 October 2019; Published: 24 October 2019

**Abstract:** In previous studies we produced autism like behavioral changes in mice by Valproic acid (VPA) with significant differences between genders. S-adenosine methionine (SAM) prevented the autism like behavior in both genders. The expression of 770 genes of pathways involved in neurophysiology and neuropathology was studied in the prefrontal cortex of 60 days old male and female mice using the NanoString nCounter. In females, VPA induced statistically significant changes in the expression of 146 genes; 71 genes were upregulated and 75 downregulated. In males, VPA changed the expression of only 19 genes, 16 were upregulated and 3 downregulated. Eight genes were similarly changed in both genders. When considering only the genes that were changed by at least 50%, VPA changed the expression of 15 genes in females and 3 in males. Only Nts was similarly downregulated in both genders. SAM normalized the expression of most changed genes in both genders. We presume that genes that are involved in autism like behavior in our model were similarly changed in both genders and corrected by SAM. The behavioral and other differences between genders may be related to genes that were differently affected by VPA in males and females and/or differently affected by SAM.

**Keywords:** ASD; epigenetics; mice; postnatal VPA injection; SAM; gene expression; nanostring

#### **1. Introduction**

Autism spectrum disorder (ASD) is a complex neurodevelopmental disorder characterized by impaired social communication and social interactions and restrictive stereotyped behaviors and interests. ASD core symptoms are frequently accompanied by developmental delay, anxiety and cognitive deficits [1–3]. Currently, the prevalence of ASD is 1 in 59 children, at an earlier onset of children prior to 3 years of age [4]. ASD is more prevalent in males, with the average male to female ratio 4:1 [4]. This male bias that can be partially explained by specific genetic differences between males and females had various impacts on both research and clinical practice [5]. Furthermore, autistic male and female demonstrate different phenotypes; females were usually reported to have lower social and communicative problems and fewer restricted and repetitive behaviors than males [6–8]. The age

of diagnosis is also on average later in autistic females than males, and females may often be missed by current diagnostic procedures [9].

Although ASD is a highly heritable disorder, prenatal and early postnatal exposures to environmental toxicants (i.e., pollutants, insecticides and pesticides) or maternal infections as well as epigenetic alterations are hypothesized to contribute to the development of ASD. Nevertheless, the underlying cause is still unclear [1–3,10]. Moreover, currently there are no biochemical or molecular markers for the diagnosis of ASD.

Animal models enable the use of preclinical tools to understand the role of genetic mutations and environmental factors in the etiology of ASD. Hence, there are several rodent models in which ASD-like symptoms were produced by exposure of the pregnant dams to a teratogenic agent, especially valproic acid (VPA), during different stages of gestation, and a few following early postnatal insults [2,10,11].

VPA is an antiepileptic drug that has also been used for bipolar disorder, migraine headaches and as a mood stabilizer [12,13]. VPA exposure during pregnancy, especially during the organogenesis period, is associated with various teratogenic effects including high risk of major anomalies, mainly spina bifida. VPA during pregnancy also leads to 5–10 fold increased rate of ASD in the offspring [3,14]. Similarly, rodents exposed to VPA during different stages of pregnancy exhibit behavioral deficits including delayed developmental milestones, stereotypic and self-injurious behavior and impaired social behavior similar to those observed in autistic patients [11,15–19]. Therefore, VPA rodent model became a robust, widely used, environmental preclinical model of ASD with face construct and predictive validity [20].

Interestingly, while the male to female ratio observed in human ASD is 4:1, the male to female ratio in ASD children exposed to VPA during pregnancy is nearly 1:1 [21] or 2:1 [22].

As in human, there are behavioral differences between male and female rodents exposed to VPA. Male mice exposed to VPA exhibited social impairment and reduced social interaction, manifested by lack of preference for a stranger mouse in the three-chamber test [23–25]. On the other hand, increased repetitive behavior and anxiety-like behaviors was reported in both males and females [23,26].

It is important to remember that most cases of ASD in human are genetic in origin and not related to the use of specific drugs or chemicals during pregnancy. There are also distinct neurobehavioral differences in human with ASD when compared to animal models. For example, there is a wide range of cognitive abilities in human from normal to severe retardation (1,2), while learning is often impaired in the VPA animal model.