Participant characteristics

Study samples were pulled from the Roots of Autism and ADHD Twin Study in Sweden (RATSS)23. Participant recruitment and tooth collection within RATSS have been previously described24. Briefly, we collected and analyzed teeth from 74 participants: 30 complete twin pairs, 1 triplet group, and 11 individuals from twin pairs whose sibling did not donate a tooth. For our primary analysis, we measured metals in 13 cases of ADHD, 8 ASD, 12 ADHD/ASD combined, and 41 neurotypical controls (see Table 1 for participant characteristics). This sample size accounts for ~50% of the whole RATSS cohort of tooth shedding age. The clinical aspect of the study and sample collection were approved by the Swedish National Ethical Review Board. All participants gave informed consent. Analyses were also approved by the Institutional Review Board of the Icahn School of Medicine.

Table 1 RATSS Participant Characteristics. (A) Participant numbers, gestational age and birth weight. (B) Mean IQ (standrd deviation) by diagnostic category Full size table

Clinical assessment

Participants were diagnosed using a consensus process with several experienced clinicians, and according to DSM-5 criteria, endorsed by information from the following standardized instruments: Kiddie Schedule for Affective Disorders and Schizophrenia (K-SADS25), Diagnostic Interview for ADHD in Adults (DIVA 2.026), Autism Diagnostic Observation Schedule 2nd Edition (ADOS-227), and Autism Diagnostic Interview- Revised (ADI-R28). IQ testing was performed with the following measures depending on participant’s age and verbal capacities: Wechsler Intelligence Scale for Children-IV (WISC-IV29), Wechsler Adult Intelligence Scale-IV (WAIS-IV30), and the Leiter International Performance Scale-Revised30).

Laboratory analysis

Our approach to laser ablation-inductively coupled plasma-mass spectrometry (LA-ICP-MS) tooth metals analysis and assigning developmental times has been detailed elsewhere31,32. Briefly, teeth are sectioned vertically along the labio-lingual/buccal-lingual plane and sampled parallel to the dentine-enamel junction from the dentine horn tip towards the tooth cervix. Depending on scan length, each tooth is sampled at 152 locations on average. Temporal information is assigned to sampling points using the neonatal line, a histological feature formed in enamel and dentine at the time of birth, and additional pre- and postnatal incremental markings. A New Wave Research NWR-193 (ESI, USA) laser ablation unit equipped with a 193 nm ArF excimer laser was connected to an Agilent Technologies 8800 triple-quad ICP-MS (Agilent Technologies, USA). Ablation was carried out under a helium atmosphere which is mixed with argon via Y-piece before introduction to the ICP-MS. Instrument sensitivity (maximum analyte ion counts), oxide formation (232Th16O+/232Th+, <0.3%) and fractionation (232Th+/238U+, 100 ± 5%) were monitored daily using NIST SRM 612 (trace elements in glass). Data were collected as metal to calcium ratios (e.g., 208Pb:43Ca) to control for any variations in the mineral content within a tooth and between samples.

Recurrence quantification analysis

As previously described in Curtin et al.19,33, we used non-linear methods—recurrence quantification analysis (RQA) and cross-recurrence quantification analysis (CRQA)—to characterize dynamic, cyclical properties of environmental exposures and their metabolism. Briefly, time-series data sampled from teeth were used to construct recurrence plots of single elements or cross-recurrence plots of two elements. We investigated the temporal features generated by recurrence plots, focusing on determinism, mean diagonal length, entropy, and recurrence time. These features measure diagonal line structures, or cyclical events, in recurrence plots. Determinism measures the ratio of diagonal lines (cyclical events) to vertical and/or horizontal lines (non-cyclical events) effectively defining the periodicity of an elemental time-series. Mean diagonal length measures the mean duration of diagonal lines or cycles, and recurrence time measures the interval between these cycles. Entropy captures the complexity of cyclic activity by measuring the variability in cycle lengths. These measures were similarly generated for CRQA, but captured the temporal relationship between two elemental signals. These methods are summarized in Fig. 1 and were processed with the Cross-Recurrence Toolbox v5.1634 in Matlab v2016b (Mathworks).

Fig. 1: Overview of study design. a Collected deciduous teeth were analyzed using laser ablation-inductively coupled plasma-mass spectrometry to generate temporal profiles of metal uptake during fetal and postnatal development. b Example cobalt exposure profile in a typically-developed (TD) control subject ranging from −106 to 110 days since birth. c Example cobalt exposure profile in a subject diagnosed with ADHD ranging from −125 to 187 days since birth. d Recurrence plot generated from control element trace in panel b. This graphical analytical tool presents cycIical processes as diagonal lines; in recurrence quantification analysis (RQA), the duration (mean diagonal line length, MDL), complexity (entropy), and determinism (proportion of diagonal elements) of cyclical processes are analyzed to characterize dynamical features in elemental time series. e Recurrence plot generated from ADHD case in panel c. The relative abundance of laminar states (black structures) and attenuated diagnoal length indicate reduced periodicity Full size image

Statistical methods

Prior to applying inferential statistical testing, all variables were evaluated to confirm assumptions of normality in value distributions. Values of ±2 SDs were excluded from some variables to meet assumptions of normality. For analyses of recurrence features (RQA/CRQA, described above), linear mixed models were used to test main dichotomous effects of the presence/absence of an ADHD-diagnosis on these features. Twin pairs were modeled as random variables, while also controlling for sex, gestational days, birth weight, and comorbid ASD status. Additional covariates, such as zygosity and IQ, were initially included in model construction but were ultimately excluded as these yielded no improvement in model fit (Akaike Information Criteria, AIC), had no significant effects, or caused no changes in the significance of other covariates included in modeling. False discovery rate (FDR) adjustments, stratified by metal pathway (or cross recurrence), were applied to raw p-values relating to ADHD effects; unless otherwise noted, all reported p-values reflect FDR adjustment. Adjusted p-values less than or equal to 0.05 are reported as significant.

We additionally implemented dimensionality-reduction techniques to investigate the utility of these methods in characterizing neurotypical, ADHD, ASD, and ADHD/ASD comorbid phenotypes. We first applied principal component analysis (PCA), an unsupervised dimensionality-reduction technique, to investigate relationships among features derived from RQA, and evaluate the efficacy of derived components. Measures derived from RQA of single elements and cross-recurrence (CRQA) analyses, including Determinism, Entropy, Mean Diagonal Length, and Recurrence Time, were centered and scaled for PCA. Following inspection of a scree plot, 15 components were retained for subsequent analyses, as these explained >80% of total variance and all estimated eigenvalues were greater than 1. Linear models were used to evaluate if derived components differed between subjects diagnosed as neurotypical, ASD, ADHD, or comorbid with ASD and ADHD. PCA was implemented in R v 3.5.2 with the factoextra package.

We next complemented the unsupervised analysis described above with a linear discriminant analysis (LDA), a supervised dimensionality-reduction technique, to evaluate the separation of ADHD, ASD, comorbid ASD/ADHD and neurotypical behavioral phenotypes on the basis of elemental metabolic features extracted via RQA. This analysis simultaneously evaluated all 60 features derived from eight elemental pathways and seven cross-recurrences to separate behavioral diagnoses along a lower dimensional space. Subjects with missing values for any of these features were excluded from this analysis. For each discriminant axis calculated, the correlation between raw RQA features and discriminant scores was calculated to provide a standardized measure of variable loading/importance as it related to a given axis. Analyses were implemented in R v 3.5.2 with the mass package.