Saliva RNA Biomarkers of Gastrointestinal Dysfunction in Children With Autism and Neurodevelopmental Disorders: Potential Implications for Precision Medicine

Gastrointestinal (GI) disorders are common in children with neurodevelopmental disorders such as autism spectrum disorder (ASD). A limited understanding of the biologic factors that predispose this population to GI disorders has prevented development of individualized therapies to address this important medical issue. The goal of the current study was to determine if elements of the salivary micro-transcriptome could provide insight into the biologic perturbations unique to children with ASD-related GI disturbance. This cohort study included 898 children (ages 18–73 months) with ASD, non-ASD developmental delay (DD), or typical development (TD). The saliva micro-transcriptome of each child was assessed with RNA-seq. Outputs were aligned to microbial and human databases. A Kruskal Wallis analysis of variance (ANOVA) was used to compare levels of 1821 micro-transcriptome features across neurodevelopmental status (ASD, DD, or TD) and GI presence or absence. An ANOVA was also used to compare micro-transcriptome levels among GI sub-groups (constipation, reflux, food intolerance, other GI condition, no GI condition), and to identify RNAs that differed among children taking three common GI medications (probiotics, reflux medication, or laxatives). Relationships between features identified in ANOVA testing were examined for associations with scores on the Autism Diagnostic Observation Schedule, 2nd Edition (ADOS-2) and the Vineland Adaptive Behavior Scales. GI disturbance rates were higher among children with ASD than peers with TD but were similar to those with DD. Five piwi-interacting RNAs and three microbial RNAs displayed an interaction between developmental status and GI disturbance. Fifty-seven salivary RNAs differed between GI sub-groups–with microRNA differences between food intolerance and reflux groups being most common. Twelve microRNAs displayed an effect of GI disturbance and showed association with GI medication uses and measures of behavior. These 12 microRNAs displayed enrichment for 13 physiologic pathways, including metabolism/digestion long-term depression, and neurobiology of addiction. This study identifies salivary micro-transcriptome features with differential expression among children with ASD-related GI disturbance. A subset of the RNAs displays relationships with treatment modality and are associated with autistic behaviors. The pathobiologic targets of the micro-transcriptome markers may serve as targets for individualized therapeutic interventions aimed at easing pain and behavioral difficulties seen in ASD-related GI disturbance.


INTRODUCTION
Previous work has demonstrated that the salivary microtranscriptome (non-coding RNA and microbial RNA) could be used to distinguish between children with autism spectrum disorder (ASD) (ages 2-6 years) and peers with typical development or developmental delay (1). These non-coding RNAs have regulatory roles in metabolism, cell differentiation and neuronal differentiation, by inhibiting gene expression (2). Elements of the micro-transcriptomes are up-or downregulated by cells in response to the external environment (3), which suggests that they may have a dynamic relationship with other factors. Several non-coding RNA in saliva have demonstrated that their levels are associated with adaptive and autistic behaviors in children with ASD (4,5), and are associated with socialization and autistic behaviors in young children with ASD (5). However, the relationship between non-coding RNA and comorbid conditions in ASD is not yet known.
Children with ASD appear to more frequently experience GI conditions than their neurotypical peers. Children with ASD have been reported to be diagnosed with a GI problem almost four times more often than children without ASD (6). The range of reported prevalence of GI symptoms is from 9 to 91% (7), likely a result of different methods of GI assessment. Constipation and diarrhea tend to be the most common GI diagnoses in ASD (6), with constipation the most common (8). Constipation can frequently be sufficiently severe to result in emergency department visits and hospital admissions among children with ASD (9).
Children with abdominal pain can also manifest difficult and distressing behaviors such as irritability, social withdrawal, stereotypy, and hyperactivity, as well as aggression and selfinjurious behaviors (7,10,11). Associated comorbid conditions can include seizures, anxiety, depressed mood, attentiondeficit/hyperactivity disorder, oppositional defiant disorder, sleep problems, as well as other problem behaviors (12)(13)(14)(15)(16). Problem behavior may, itself, sometimes be an indicator of GI distress in ASD, particularly among individuals with ASD with limited language (7). Younger individuals with ASD and GI disturbances display more externalizing behaviors such as aggression, and older individuals with ASD display more internalizing symptoms such as anxiety and depression (16). Stress reactivity as well as anxiety and autonomic arousal are also critically interrelated to severity of GI symptoms in ASD (17,18). Many ASD patients with GI disturbances, such as constipation, are less likely to respond to first line therapies, such as stool softeners (19). Therefore, a better understanding of the biologic processes driving GI disturbances in patients with ASD might provide mechanistic insights toward better treatments for GI pain and related behaviors.
Heterogeneity across the autism spectrum has led to failures in many of the early clinical trials attempting to target core features of ASD (20). While micro-transcriptomes appear to distinguish ASD patients from neurotypical controls (21), they may also have particular value in helping to distinguish specific subtypes of ASD, which might impact treatment. Because of the significant adverse effects of GI disorders in ASD (6,9), as well as their interrelationships with behavioral disturbances (12)(13)(14)(15)(16), we sought to examine differential micro-transcriptome expression in ASD patients with and without GI disturbances, as well as how this interrelates with behaviors. A greater understanding of the downstream targets of the differentially expressed microtranscriptomes may help guide future personalized medicine approaches to treatment of GI disturbances in ASD.

Ethics
The study was approved by the Western Institutional Review Board (IRB #20180172). Written consent was obtained for all participants, and written informed assent was documented for those capable of assent. DSM-5 criteria in association with standardized assessment tools (i.e., the Autism Diagnostic Observation Schedule, 2nd Edition; ADOS-2). DD participants included children referred for initial ASD assessment who did not meet diagnostic criteria, as well as children with negative ASD screening tools who required early intervention services for delays in gross motor, fine motor, language, or cognitive development (as reported by parental survey and confirmed through review of the medical record). TD participants included children recruited at the time of their annual well child visits who did not exhibit developmental delays on standard developmental surveillance tools (e.g., Survey of Wellbeing in Young Children, Parents' Evaluation of Developmental Status-Developmental Milestones, Modified Checklist for Autism in Toddlers Revised). Participants were further subdivided by presence (n = 184) or absence (n = 714) of gastrointestinal (GI) disturbance, based on: (2) parent report of (a) constipation (n = 84); (b) reflux (n = 46); (c) chronic diarrhea or abdominal pain (n = 22); (d) food intolerance (n = 45); (e) cyclic vomiting/dysphagia (n = 3); or (f) eosinophilic esophagitis (n = 7). Note that total numbers of specific conditions exceed the number of participants with a GI disturbance because 23 participants reported more than one type of GI disturbance. Exclusionary criteria for all participants included Ward of the State, g-tube dependence, active periodontal infection or acute upper respiratory illness.

Participant Characteristics
The ADOS-2 was administered by trained raters to children with ASD (n = 409), and children with DD (n = 121) in whom ASD was suspected. The Vineland Adaptive Behavioral Scales 3rd Edition (VABS-III) was used to measure adaptive behavior, communication, and social interaction for all participants. Additionally, medical and demographic information including sex, age, race, ethnicity, medical conditions, and medications, was collected through parental surveys and affirmed via review of the electronic health record where available.

Saliva Collection and Processing
Saliva was obtained from all participants in a non-fasting state via swab, targeting the base of the tongue and between the gums and buccal mucosa as locations for the collection using an Oracollect RNA swab (DNA Genotek, Ottawa, Canada). Nucleic acid extraction was performed using the Qiagen miRNeasy Microkit (Cat. No. 217084), a QIAzol based purification method. The RNA sequencing process included using an Illumina TrueSeq Small RNA Prep protocol for library construction, followed by sequencing on an Illumina flow-cell and a NextSeq 500 instrument (Illumina; San Diego, CA, United States). Sequencing outputs were a binary base call (BCL) sequence file per sample, which was then converted to a FASTQ file, a text-based format that includes detected bases and associated quality scores (i.e., confidence in correct detection). Alignment and quantification of known RNA sequences for each collected specimen was done using the Bowtie1 aligner (22) to the following reference databases: miRBase v22 (23), piRBase v1 (24), RefSeq v90 (25), and hg38. Quantification of the detected sequences yielded counts of known human micro-ribonucleic acids (miRNAs), long non-coding transcripts (small nucleolar RNAs), and piwiinteracting RNA (piRNAs). To determine microbial RNAs present in the sample, the leftover sequences that did not align to hg38 were aligned to the NCBI microbial database using k-SLAM, an efficient aligner used in metagenomic data. Aligned sequences were then assigned to microbial genes, which were quantified to a microbial identity (e.g., genus, species, strain). Prior to analysis and count normalization, low count RNAs were removed from further analysis so that only reliably expressed RNAs were interrogated. Tabulated counts of each RNA were compared to the total counts detected in that RNA category, and RNAs that did not account for at least 0.01% of the total were dropped. Following abundance filtering, the remaining RNAs were quantile normalized and mean-center scaled.

Statistics
The primary goals of the study were to: (2) identify human and microbial RNA levels in saliva that were associated with GI disturbance; (3) investigate whether these relationships were impacted by child developmental status; and (4) determine if specific RNA "biomarkers" displayed unique expression patterns in particular GI disturbances (e.g., constipation) or with particular treatments (e.g., probiotics). A two-way analysis of variance (ANOVA) was used to compare levels of 1821 RNA among the 898 participants based upon two factors: (2) neurodevelopmental status (ASD, DD, or TD); and (3) GI status (presence or absence of any GI condition). Interactions between neurodevelopmental status and GI status were reported. A one-way Kruskal Wallis rank sum test was used to identify RNAs that differed among GI sub-groups (constipation, reflux, food intolerance, other GI condition, no GI condition), and to identify RNAs that differed among those taking three common GI medications (probiotics, reflux medication, or laxatives). Finally, given the potential associations between underlying GI disturbance and child behaviors, relationships between RNAs identified in ANOVA testing, as well as the all oneway Kruskal Wallis testing were examined for associations with scores on the ADOS-2 and Vineland using Spearman Rank Testing. Benjamini Hochberg multiple testing correction was applied to all analyses. Functional analysis of candidate miRNAs (features displaying an interaction effect between neurodevelopmental status and GI disturbance, as well as relationships with treatment or behavior) was performed in DIANA miRPATH software v3.0 (26). The microT-cds algorithm (0.95 microT Threshold) was used to identify pathways overrepresented by putative messenger RNA targets by Fisher Exact Test with Benjamini Hochberg multiple testing correction. Additionally, the rates of different demographic features were tested in the population. To test for differences in age by diagnosis (ASD, DD, TD) or presence of a GI disturbance, a one-way ANOVA was used. To test for differences in rates of sex, race, ethnicity, GI disturbance, constipation, reflux, food intolerance, chronic abdominal pain, diarrhea, or eosinophilic esophagitis, a chi-squared test was used yielding the chi-squared test statistic (x) and the associated pvalue (p).  (5) 23 (4) 11 (5) 11 (5) Chronic abdominal pain or diarrhea, # (%) 22 (2) 16 (3) 6 (2) 1 (0.5) Cyclic vomiting or dysphagia, # (%) The number of participants with specific GI disturbances exceeds the total number of participants with any GI disturbance (n = 184) because 22 participants reported more than one GI disturbance. *Denotes significant difference (p < 0.05) compared with ASD group on chi-square testing.

Impact of GI Disturbance on Saliva RNAs
Among the 1821 RNA features interrogated, 28 displayed a significant difference (adj p < 0.05) between children with/without GI disturbance (Table 2A). These RNA features included four mature miRNAs and 24 small non-coding RNAs, but no microbial RNAs. There were eight RNA features that displayed a significant interaction effect between neurodevelopmental status (ASD/DD/TD) and presence/absence of GI condition (Table 2B). These RNA features included five piRNAs and three microbial RNAs (Figure 1). The piRNAs tended to display similar saliva levels across ASD/DD/TD groups without GI disturbance, but were lower among children with ASD and GI disturbance, relative to peers with TD and GI disturbance.

Differences in Saliva RNA Levels Among GI Phenotypes
There were 57 RNA features that differed between GI phenotypes ( Table 3). These RNA features included 12 microbial RNAs, three piRNAs, and 42 miRNAs. The largest differences tended to occur in miRNA levels, and were most common between children with reflux and food intolerance (Figure 2).

Effect of Medications on Saliva RNAs
Levels of 65 RNA features differed among children with GI disturbance on probiotics (n = 22) and children with GI disturbance not taking probiotics (n = 162) (Supplementary Table 1). These RNA features included 37 miRNAs, 75 piRNAs, one small non-coding RNA, and one microbial RNA. Levels of 53 RNA features differed among children with GI disturbance on laxatives (i.e., polyethylene glycol; n = 39) and children with GI disturbance not taking laxatives (n = 145). These RNA features included 15 microbial RNAs, seven small non-coding RNAs, four piRNAs, 27 miRNAs.

DISCUSSION
In this cohort study of 898 children, rates of GI disturbance were higher among children with ASD than peers with TD, as expected (6), but were similar to those with DD. There were five piRNAs and three microbial RNAs in saliva that displayed an interaction between developmental status and GI disturbance (Figure 1). These features may serve as biomarkers for the unique pathophysiology leading to elevated GI disturbance in children with ASD. There were many salivary RNAs whose levels differed between GI disturbance phenotypes-with miRNA differences between food intolerance and reflux groups being most common. Levels of 12 salivary miRNAs that displayed an effect of GI disturbance were also associated with GI medications and measures of child behavior (miR-1307-5p, miR-141-3p, miR-142-5p, miR-148a-5p, miR-186-5p, miR-200a-3p, miR-200a-5p, miR-23a-3p, miR-23b-3p, miR-28-3p, miR-532-5p, and miR-769-5p). The 12 salivary miRNAs that displayed relationships with GI disturbance, GI medications, and child behavior may serve as examples of biologic targets for personalized diagnostic and therapeutic approaches in children with ASD-related GI disturbance. Putative targets of these 12 miRNAs include transcripts that code for key regulators of both metabolism (e.g., steroid biosynthesis, porphyrin metabolism, ascorbate metabolism, calcium reabsorption, thyroid hormone signaling) and neurobiology (e.g., long-term depression). Intriguingly, exogenous steroids, porphyria, hypercalcemia, hypothyroidism, and depression are all associated with constipation and/or abdominal pain. It is possible that the 12 miRNAs contribute to sub-clinical perturbations in these physiologic pathways, in somuch-as they lead to GI pain without causing other overt clinical symptoms. For example, rodent models have demonstrated that restoration of miR-148a expression in the lower GI tract may reduce colitis (27), while elevations in miR-200a may lead to irritable bowel-like symptoms through inhibition of serotonin and cannabinoid transporters (28).
Our understanding of the nature of GI disturbances in ASD is only beginning to emerge. Numerous pathways, including autonomic arousal (17,18), serotonin dysregulation (29), and perturbations in gene expression (30) have all been implicated in this process. The micro-transcriptome features identified in this study provide a single mechanism through which each of these pathways may converge (Figure 3). For example, recent research has found that stress reactivity, anxiety, and autonomic arousal are interrelated with the severity of lower GI symptoms in ASD (17,18). One miRNA identified in this study, miR-142-5p, has been previously implicated in anxious behavior following prolonged stress (31). Whole blood serotonin levels have also been associated with lower GI symptoms in ASD (29). Here, we identify one miRNA (miR-23a-3p) implicated in ASD-related GI pathology, which has previously been shown to change in depressed patients treated with selective serotonin reuptake inhibitors (32). Specific genes, in particular polymorphisms of the Mesenchymal Epithelial Transition (MET) receptor kinase gene, are also associated with GI symptoms in ASD (30). The MET receptor has been shown to modulate miRNA expression (33), and the MET transcript is a putative target of two miRNAs in the present study (miR-23a-3p, miR-23b-3p).
Immunological factors have also been found to be associated with GI symptoms as well as behavior in ASD (34,35). The unique immunologic patterns associated with ASD may contribute to specific alterations in the microbiome profile that have been reported in children with ASD and GI symptoms (36). Evidence of this nature has even led to efforts at interventions based on the microbiome in ASD (37,38). In the present study, we identify several miRNAs that are implicated in immune development. For example, miR-28 has been shown to modulate T-cell differentiation and cytokine expression (39), while miR-200a-3p and miR-141-3p have been found to work together to modulate differentiation of interleukin-producing T-helper cells (40). We found minimal overlap between the specific microbes identified in this study, and those of previous GI microbiome studies (36). This may be because the current investigation examines microbial RNA levels in saliva, as opposed to the more traditional 16S approach, using stool samples.  As further research begins to reveal a clearer understanding of the pathways implicated in ASD patients with GI symptoms, including those specific to the ASD/GI population, we can begin to understand why some patients with ASD appear to respond less reliably to treatment than others with similar GI symptoms (7). Examining the downstream targets of miRNA differentially expressed in those with ASD and GI symptoms would likely contribute substantially to this understanding. Fortunately, with the ability to rapidly, and non-invasively measure miRNA in saliva (21), such information can readily be obtained from large populations. Development of such targeted approaches may provide opportunities for personalized treatment of gastrointestinal symptomatology, and lead to down-stream improvement in related behaviors (7,10,11), by impacting anxiety, mood, sleep, and attention (12)(13)(14)(15)(16).
This study harnesses, to our knowledge, the largest sample of the salivary micro-transcriptome in ASD. Its inclusion of children with non-ASD developmental delay as part of the control group is a relative strength. Inclusion of participants from multiple geographic sites also promotes generalizability of the findings. However, there are several limitations which should be noted. First, GI disturbances were identified through parent report and review of medical records but were not specifically assessed by physicians as part of the study. Second, we acknowledge that "GI disturbance" is a somewhat artificial distinction that groups together loosely related pathology that occurs in the GI tract. Important physiologic differences exist between conditions such as constipation and reflux, and these underlying biologic differences may have served to enhance false negative findings in the initial analyses. For this reason, we performed secondary analyses of the GI sub-phenotypes. However, this approach also has a trade-off of reducing the study's substantial sample size. Third, there are several pre-analytic factors which may potentially confound this study's findings, including batch effects and sample collection factors. We note that all samples were run on the same sequencing machine, using the same library preparation procedure, performed by the same laboratory technician. Although this analysis did not control for sample collection parameters, such as collection time, prandial status, or prior tooth-brushing, we have previously assessed the impact of many of these factors on the saliva microtranscriptome (21). We note that none of FIGURE 2 | The most significant (12) transcripts differing between GI phenotypes. The whisker box plots represent normalized abundance of miRNAs and microbial RNAs that displayed a significant difference between children without a GI condition (orange), children with constipation (green), food intolerance (blue), reflux (pink), or another GI condition (brown). Significance levels-p-val <0.001 (*), p-val <0.0001 (**). Gap junction 0.001139 10 7 Porphyrin and chlorophyll metabolism 0.005067 9 3 Endocrine and other factor-regulated calcium reabsorption 0.005067 6 7 Drug metabolism-cytochrome P450 0.007162 8 4 Glycosphingolipid biosynthesis-lacto and neolacto series 0.009386 4 5 Ascorbate and aldarate metabolism 0.016999 7 2 Lysine degradation 0.025194 5 7 Proteoglycans in cancer 0.026807 22 9 Thyroid hormone signaling pathway 0.027249 10 7 Long-term depression 0.030698 6 4 Morphine addiction 0.034341 9 7 the microbial features and very few of the miRNA features identified in this study have demonstrated relationships with pre-analytic factors.
We recognize that the salivary transcriptome serves as a proxy for the primary pathologic site of most GI disturbances, the lower GI tract. However, several studies have reported FIGURE 3 | Putative micro-transcriptome links between biologic mechanisms and physiologic pathways involved in gastrointestinal pathophysiology. The conceptual diagram connects biologic pathways that have been previously implicated in Autism Spectrum Disorder (ASD)-related gastrointestinal (GI) pathology with putative micro-transcriptome features (i.e., microRNAs and microbes). Micro-transcriptome features were included based on their associations with the biologic pathways of interest in the existing literature. Collectively, the microRNAs display target enrichment for several clinical findings (i.e., thyroid hormone signaling, long-term depression, etc.), that are known contributors to GI disturbance. significant overlap between saliva and stool micro-transcriptome features (41)(42)(43). Unlike the stool microbiome, we note that the saliva microbiome can be repeatedly sampled on demand and has shown resilience to antibiotic treatments (44). These characteristics make saliva an attractive source for sampling GI-related biology (particularly in patients with conditions such as reflux, eosinophilic esophagitis, or cyclic vomiting). Our previous work with parents of children with ASD has shown that they overwhelmingly prefer saliva as a clinical biofluid (45). Additionally, we note that association analyses between salivary transcriptome elements and ADOS scores rely solely on ASD and DD participants for whom these assessments were available. The lack of TD participants in this analysis could have impacted the findings. Finally, we must also consider the possibility that some of these markers may be caused by the downstream effects of the gastrointestinal symptoms or treatments, rather than serving a mechanistic role, but this would not diminish their potential use as biomarkers.
This is, to our knowledge, the first effort to examine the salivary RNA profile associated with GI symptoms in ASD, in a large population study. With the increased understanding of the critical importance of subtyping for meaningful precision medicine approaches in ASD (20), as well as the importance of GI symptomatology in behavioral issues in ASD (7,10,11), and the potential for mechanistic understanding through examination of the downstream targets of differentially regulated miRNAs, this is an important future direction of investigation. A subset of the micro-transcriptome features identified in this study displays relationships with treatment modality and are associated with autistic behaviors. The pathobiologic targets of these micro-transcriptome markers may serve as novel targets for individualized therapeutic interventions aimed at easing pain and behavioral difficulties seen in ASD-related GI disturbance.

DATA AVAILABILITY STATEMENT
The datasets presented in this article are not readily available due to the proprietary nature of the research. Requests to access the datasets should be directed to David Levitskiy, david.levitskiy@motionintel.com.

ETHICS STATEMENT
The studies involving human participants were reviewed and approved by the Western Institutional Review Board (IRB #20180172). Written informed consent to participate in this study was provided by the participants' legal guardian/next of kin.

AUTHOR CONTRIBUTIONS
FM, RG-K, RS, DB, KS, and SH contributed to conception and design of the study. DL performed the statistical analysis. AC and PT contributed to data collection and formatting. DB, KS, and SH wrote the first draft of the manuscript. All authors contributed to manuscript revision, read, and approved the submitted version.

FUNDING
This work was funded by an award from the National Institutes of Health to Quadrant Biosciences (R42MH111347).