Space-Time Analysis of Burgeoning US Atrial Septal Defect Rates Driven by Cannabis

Article type: Research Article

Keywords: cannabis, cardiogenesis, genotoxicity, epigenotoxicity, transgenerational inheritance

Authors: Albert Stuart Reece ⚬, Gary Kenneth Hulse

Affiliations: Division of Psychiatry, University of Western Australia, Stirling Hwy, Crawley, Perth, WA 6009, Australia; School of Medical and Health Sciences, Edith Cowan University, Joondalup Dr., Joondalup, WA 6027, Australia

License: © 2026 by the authors. CC BY 4.0 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.

Article links: DOI: 10.3390/jox16020068  |  PubMed: 42041866  |  PMC: PMC13117238

Relevance: Moderate: mentioned 3+ times in text

Full text: PDF (1.9 MB)

1. Introduction

Recent quinquennial reports from the Centers for Disease Control (CDC)-affiliated National Birth Defects and Prevention Network (NBDPN) show an extraordinary variation in USA atrial septal defect (ASD) secundum rates (ASDR) from 0.4/10,000 live births in Maryland to 884, 849, 802, 769, 772, and 739 /10,000 live births in various ethnic groups in Nevada and Mississippi [ref. 1]. It is noteworthy that such extraordinarily high and increasing rates of teratogenicity, involving up to almost one in eleven births, have not been explored in the broader neonatal epidemiological literature.

1.1. Large Population Studies

The situation is rendered more poignant because several large population-based studies have linked the recent upsurge in cannabis use across Europe and the USA with ASD [ref. 1,ref. 2,ref. 3,ref. 4,ref. 5,ref. 6,ref. 7], suggesting the two trends may be related. The first large study reviewed 316,508 live births in the Hawaii state registry 1986–2002 [ref. 3]. Prenatal cannabis use was identified in 829 cases based upon self-reports and urine toxicology screens. In five cases where cannabis was the only drug used, the rate ratio of ASDR compared to the general control group was 6.12 (95%C.I. 1.98, 14.35). In a group of 12 cases where cannabis use occurred either alone or in concert with other drugs, the rate ratio was 6.69 (3.44, 11.76).

A review of 1,015,302 births in Colorado in the 2000–2014 period found that ASDR rose significantly from 299 in 2000 to 912 in 2012 (threefold; p < 0.0001) when the birth rate was only rising by 3.3% annually and across a period when cannabis exposure was rising but other exposure to other drugs was stationary or declining [ref. 4,ref. 8]. Using multiple regression in linear models revealed that cannabis exposure was significant, both as a main effect and in interaction with tobacco as a predictor of major cardiovascular anomalies. In regression models quartic in time, Δ9-tetrahydrocannabinol (THC) was significant as a main effect, and cannabis-time interactions were significant in first, second, third, and fourth powers of time for major cardiovascular anomalies.

A combined geospatial and causal inferential study of birth defects in the USA in 1991–2016 found that the mean nationwide ASDR rose three-fold from 27.4 to 82.8/10,000 live births during a period when there were 73,179,872 births [ref. 2]. The main focus of this study was the period 2003–2016, a period of 57,050,944 births and for which drug use data are available from the National Survey of Drug Use and Health (NSDUH). Drug use, ethnicity and median income by state were assessed. Cannabis was found to be significantly associated with ASDR in geospatial models as a main effect and in interaction with tobacco. In inverse probability-weighted robust additive and interactive generalized regression models, cannabis and ethnic cannabis exposure showed the strongest association with ASDR. Ten E-values were documented, mostly in the range of 2–3. In all cases, the lower 95% bound of the E-values was above the threshold of causality at 1.25 [ref. 9]. On the basis of consistency with other studies, the existence of multiple biologically plausible mechanisms, strength of association, fulfillment of the Hill criteria for causality, an unequivocal temporal sequence, inverse-probability weighting and pseudo-randomization of measured predictors and E-value constraint of unmeasured confounders, this relationship was found to be causal.

A study of 2,838,963 Canadian births from 1989 to 2009 found that cardiovascular defects were elevated in high-cannabis-using territories compared to the provinces with less cannabis use in the south of the country [ref. 7,ref. 10]. Since ASD is one of the most common cardiovascular anomalies, total cardiovascular anomalies were used as a surrogate marker for ASDR in this study. In sophisticated geospatial regression models, cannabis exposure was highly significant, with p-values < 10⁻¹⁵ both as a main effect and in interaction with tobacco and opioids and in a three-way interaction with tobacco and opioids together.

A study of twelve cardiovascular anomalies in 14 European nations in 2010–2019 with 77,410 anomaly rates found than ASDR was highly significantly related to cannabis exposure in inverse probability-weighted panel regression after adjustment for other drug exposure and median national income (p < 2.2 × 10⁻¹⁶), and this high level of significance was maintained under geospatiotemporal regression [ref. 5,ref. 11]. In inverse probability-weighted panel regression models, E-values for cannabis terms ranged up to 9.01 × 10¹⁵, with a 95% lower bound of 2.54 × 10¹³. In geospatial models, ASDR E-values ranged up to 1.54 × 10¹³⁹, with 1.07 × 10²¹ as the 95% lower bound. The relationship was said to be causal because of the close relationship with spatiotemporal regression, the strength of association, the fulfillment of Hill’s criteria for causality, the clear temporality, the presence of multiple plausible biological pathways, particularly from epigenomics, concordance with other large studies, the use of inverse probability weighting to control measured predictors by pseudo-randomization, and the very high E-values, which constrained unmeasured confounding within impossibly high limits.

An Australian study of 513,895 births in 2008–2015 found a prevalence ratio of 3.27 (2.63, 4.08) for ASDR in high- compared to low-cannabis-use regions [ref. 12,ref. 13]. The attributable fraction in the exposed group was 68.45% (61.69%, 75.47%). Elevated E-values suggested that uncontrolled confounding was unlikely to account for these effects, and the relationship was described as being causal.

1.2. Potential Implications

Cardiovascular anomalies are the largest group of congenital anomalies, among which ASD is the most common [ref. 14]. Early termination of pregnancy for anomaly (ETOPFA) is not practiced for ASD. Thus, ASD may act as a bellwether biomarker for other more serious anomalies such as chromosomal and limb anomalies for which early termination is more widely practiced and which are therefore more difficult to track longitudinally.

It should be emphasized that ASD is by one of many congenital anomalies that have recently been associated with cannabinoid exposure in both experimental and epidemiological studies [ref. 6,ref. 11,ref. 15,ref. 16,ref. 17]. Cannabinoid teratogenesis is itself a subset of the larger issue of cannabinoid genotoxicity which also includes intellectual disability, cancerogenesis including heritable cancerogenesis and cellular and organismal aging [ref. 18,ref. 19,ref. 20,ref. 21,ref. 22,ref. 23,ref. 24,ref. 25,ref. 26,ref. 27,ref. 28,ref. 29]. Therefore, the implications of the present work potentially extend far beyond fetal cardiogenesis and impact the larger consideration of cannabinoid genotoxicity, including heritable genetic and epigenomic lesions.

1.3. Study Objectives and Hypotheses

Given that the incidence of ASD seems to be rising across USA in general, a review of this issue is timely. This study explored several pre-specified hypotheses: (1) that cannabis is driving the US ASD rise in a manner similar to that described elsewhere; (2) that elevated cannabis use is responsible for the very high incidences reported in recent years; (3) that cannabis legalization is exacerbating this trend; and (4) that some of this effect may have been due to cannabis exposure through non-inhalational routes, potentially implicating environmental contamination. The study applies the power of space–time analysis to these questions augmented by statistical tools of formal causal inference and classical epidemiology.

2. Materials and Methods

2.1. Data

Data on ASDR by ethnicity was downloaded from 12 Periodical Reports of the NBDPN1989-1990 to 2016–2020, which is affiliated with the CDC, Atlanta, Georgia [ref. 1]. These reports were for quinquennial periods, so that he middle year was taken as the indicative year (e.g., for 2016–2020, 2018 was considered the indicative year). The states studied were Arkansas, Colorado, Delaware, Georgia, Illinois, Indiana, Iowa, Kentucky, Louisiana, Maryland, Massachusetts, Michigan, Minnesota, Mississippi, Nebraska, Nevada, New Jersey, New York, North Carolina, Oklahoma, Rhode Island, Tennessee, Texas, Utah, Virginia, West Virginia and Wisconsin. The quinquennia studied were 2003–2007, 2004–2008, 2005–2009, 2007–2011, 2008–2012, 2009–2012, 2010–2014, 2011–2015, 2012–2016 and 2016–2020. The method of birth defect ascertainment in each state was taken from annual NBDPN reports and was recorded as being “active”, “passive” or “mixed”.

Data on annual drug use exposure by US state was accessed from the National Survey of Drug Use and Health (NSDUH) conducted annually by the Substance Abuse and Mental Health Administration (SAMHSA) [ref. 30]. Data for last month cigarette (cigmon), last month alcohol use (alcmon), last month binge alcohol (bngalc), alcohol dependence (abodalc), last month cannabis (mrjmon), last year analgesic misuse (anlyr), and last year cocaine use (cocyr) was accessed for each year in 2005–2018, as 2003 was the first year the NSDUH collected data on cannabis use by state, so 2005 is the first indicative year for which analysis could be prepared. Drug use by ethnicity data at both the state and federal levels was also obtained from NSDUH. Median household income by state was obtained from the US Census. Cannabis legal status was derived from online sources [ref. 31,ref. 32]. Legal status was treated as categorial rather than ordinal data. Cannabinoid data (relating to Δ9THC, cannabigerol, cannabidiol, cannabichromene) was taken from published reports [ref. 33,ref. 34,ref. 35,ref. 36,ref. 37]. In each regression, the dependent variable was the log(ASDR). Predictor covariates were drug exposures and median state income, as indicated in each table.

2.2. Derived Data

Within each state, the numbers of people of each ethnicity using drugs were compared to the overall prevalence of that ethnic group in the population to derive an ethnic rate of drug use (e.g., mrjRel). This was then multiplied by the level of that drug use (e.g., mrjmon) to derive an estimate of the ethnic cannabis exposure (mrjRelmrj). This was further multiplied by the THC (tetrahydrocannabinol) content for that year at the national level to derive an estimate of average Δ9THC exposure for that ethnicity (mrjRelmrj9THC) [ref. 2,ref. 38,ref. 39,ref. 40]. Cannabis legal status was broken into four categories: illegal, medical, decriminalized and legal. The two groups, medical and decriminalized, were conflated due to small numbers in each group. Cannabis legal status was dichotomized as legal v. the rest. A semi-quantitative scale was used to plot cannabis legal status, where “Illegal” status = 0.1, “Medical” status = 0.3, “Decriminalized” = 0.6 and “Legal” = 0.98. Legal status was treated as a categorical covariate for analysis rather than an ordinal covariate. State-level cannabinoid exposure was estimated by multiplying Federal level cannabinoid concentrations by the level of cannabis used in that state.

2.3. Statistics

Data was processed in RStudio 2025.05.0 as GUI for R 4.5.0. Data were log transformed in the interest of normality assumptions. Such covariates included the ASD rate (ASDRt), median income, and cannabis, analgesic, and cocaine use. Data was centered for all continuous regressions in the interest of minimizing collinearity. T-tests were used to compare normally distributed data; otherwise, Wilcoxson nonparametric tests were used. Data was manipulated in the tidyverse (version 2.0.0) with “dplyr” [ref. 41]. Maps were drawn using “sf” (version 1.0-21) [ref. 42], and spatial data were manipulated with “spdep” (version 1.3-11) and “spatialreg” (version 1.3-6) [ref. 43,ref. 44] and analyzed with “splm” (version 1.6-5) [ref. 45,ref. 46]. Model reduction was achieved by the classical technique of removing the least significant term [ref. 47]. Custom color palettes were prepared for univariate plots. One map used the “viridis” color palette (version 0.6.5) [ref. 48]. Bivariate maps were prepared using the package “colorplaner” run in an R4.1.1 workspace [ref. 49]. E-Values were calculated with the package “EValue” (version 4.1.3) [ref. 50,ref. 51]. epiR (version 2.0.83) was used for the analysis of 2 × 2 tables [ref. 52].

2.4. Geospatiotemporal Analysis

Spatial analyses do not tolerate missing data. Up to two missing data points were tolerated for each state; 19 missing data points (251extant; total dataset N = 270) were imputed into the final dataset using within-state linear interpolation/extrapolation across time using the package “zoo” (version 1.8-14) [ref. 53] for a kriging rate of 7.6%. The appropriate error structure was determined by retaining significant error terms. Spatial weight matrices were computed using “spatialreg” (version 2.0.0) and “spdep” (version 2.0.0) [ref. 43,ref. 44]. Psi (ψ) is the parameter in spatial analysis for serial autocorrelation and was retained in all spatial models. In each case, the other major spatial parameters such as rho, lambda and phi were not significant and were therefore omitted from the final models. The major models of interest were full models containing all substance and income terms in additive and interactive models. This was true both for models, including cannabis and cannabinoids. Other models were used to establish or assess a bivariate relationship or otherwise assess the behavior of smaller groups of covariates and may be considered exploratory in nature.

2.5. Non-Overlapping Data

As noted above, some of the quadrennia of the data overlapped. To address the potential confusion caused by this, three non-overlapping data sets were chosen for a non-overlapping spatial analysis, namely 2003–2007, 2008–2012 and 2016–2020. These results are presented separately.

2.6. Survey Regression

The dependent variable was averaged over 5 years but the drug use data was available on an annual basis; this disparity was addressed by applying survey regression in the R package survey (version 4.4-2) [ref. 54]. Ten relevant periods of analysis were defined. Drug use in each state in each period was averaged. The state was incorporated into the design specification of the survey regression. All survey regressions were inverse probability-weighted using R package ipw (version 1.2.1) [ref. 55].

2.7. Temporal Lagging

In order to address the delayed onset of exposure, temporal lagging was also conducted. Temporal lags were introduced into the independent covariates using spatialreg and spatially regressed against the ASDR.

2.8. E-Values

E-values were derived from survey regression models without reliance on likelihood-based inference. Point estimates and confidence intervals were obtained from survey-weighted models using Taylor linearization for variance estimation. E-values were calculated using exponentiated coefficients and their confidence limits, as recommended for ratio-scale effect measures. Because the E-value framework does not require residual variance or likelihood-based quantities, the absence of such measures in design-based models does not preclude valid E-value computation. For estimates below the null, the corresponding inverse rate ratio and inverted confidence interval were used to compute E-values on the ≥1 scale.

2.9. Multicollinearity Diagnostics

Conventional variance inflation factors (VIFs) are defined for ordinary least squares regression and rely on assumptions of independent and identically distributed errors that do not hold in spatial regression models with spatially correlated disturbances. In spatial error and spatial maximum-likelihood models, the error covariance matrix is non-diagonal, which invalidates the auxiliary regressions on which standard VIF calculations are based [ref. 56,ref. 57]. Consequently, VIFs are not well-defined for assessing multicollinearity in spatial models. To evaluate collinearity among covariates independently of the spatial error structure, we employed Belsley–Kuh–Welsch collinearity diagnostics based on the eigenstructure of the standardized design matrix. This approach examines condition indices and variance-decomposition proportions to identify near-linear dependencies among regressors without relying on assumptions about the error process. p < 0.05 was considered significant.

2.10. Ethics

This study was given ethical clearance from the University of Western Australia Human Research Ethics committee, number RA/4/20/4724, on 24 September 2021.

2.11. Data Availability

Data including computational code is publicly available through the Mendeley data repository at doi 10.17632/4xp62s957c.5.

3. Results

Datasets for both ASD and drug exposure were available for the period 2003–2020 (indicative years 2005–2018), so this was the period for analysis. However, ASDR data was notably incomplete, particularly amongst ethnic minorities. For this reason, consideration was limited to the Non-Hispanic white ethnicity. The initial ASDR dataset amongst Americans of European background comprised data from 42 states over 11 different time periods. As geospatial techniques do not tolerate missing data, a decision had to be made between data retention and missing data completion. Over ten time periods, 15 states had complete data, five states had one missing datapoint, and seven states had two missing datapoints. It was decided to complete this dataset by kriging using linear regression within each state to complete missing data. In this way data for twelve states was completed with 19 results for a kriging rate of 7.6%. The size of the final dataset was thus 27 states over ten data periods in 2005–2018, for a final sample size of 270. The quinquennial span of the ten periods, including their indicative years, is shown in Supplementary Table S1. The data referred to 110,107 ASD cases amongst 17,751,437 live births identified in this period. Similar techniques have been published [ref. 2,ref. 38,ref. 39,ref. 40].

Supplementary Table S2 provides the baseline data for this analysis, with the groups stratified by legal status dichotomized as legal v. the rest. Many significant differences are noted between the two groups of states.

3.1. Sequential Map–Graphical Analysis

Figure 1 presents the data for ASDR. High rates in 2018 are noted in Nevada, New Mexico, Kentucky, Tennessee, Michigan, and Missouri followed by New York, Oregon and Vermont.

Supplementary Figures S1–S6 present the spatiotemporal data for cigarettes, alcohol consumption, binge alcohol use, alcoholism, analgesic misuse, and cocaine abuse. Overall, these plots show either a decline with time or stationarity. The cocaine plot shows a recent uptick in some states in the 2016–2020 quinquennium.

The space–time plot for last cannabis use is shown in Figure 2. High rates are noted in Oregon, Nevada, Colorado, California, New Mexico, Alaska, Michigan, Main, Vermont, Massachusetts and New Hampshire.

Supplementary Figures S7–S10 present the bivariate spatiotemporal data for the ASDR against cigarettes, alcohol use disorder, analgesic misuse and cocaine abuse. Reference to the colorplane key reveals that the colors pink and purple indicate that both covariates are elevated together. Green shows where both covariates are low. Other statuses are as indicated in the key. These four plots reflect the largely disjointed time course of both covariates and the declining rate of use of the substance.

Figure 3 reveals the emergence of positive purple or pink signals in Colorado in 2014 and in Nevada, Oregon, Alaska, Colorado, New Mexico and Vermont in 2018. However this graph is generally suppressed by the very elevated signals in Nevada in 2014 and 2018.

Figure 4 re-plots this data with the Nevada data excised and the remaining data re-scaled appropriately. This graph reveals a bipartite picture with positive pink or purple signals coming from Oregon, Colorado, Alaska, New Mexico, Michigan and Vermont but high rates of ASD only in the Midwest states of Kentucky, Tennessee, Mississippi and Missouri.

Supplementary Figure S11 plots the cannabis legal status. Distinct similarities to Figure 2 are noted.

Figure 5 is a semi-quantitative bivariate plot of cannabis legal status plotted against the ASD rate. This plot has obvious similarity to Figure 4, given the well-described elevated rate of cannabis use under more liberal legal paradigms [ref. 58,ref. 59].

3.2. Spatial Regression

Interstate geospatial links were computed to generate a spatial weight matrix for spatiotemporal regression (Supplementary Figure S12).

Both last month cannabis and ethnic cannabis exposure were significantly related to ASDR when considered alone in spatial regression (Table 1, first two models). When ethnic exposure to all the drugs and median income were considered, only cannabis remained in the final models (models 3 and 4). When cannabinoids were considered, both Δ9THC and cannabidiol had positive regression coefficients (model 5).

Introduction of interactions in these models produced the following results. When a tobacco: cannabis interaction was used, only the term for ethnic cannabis exposure remained positive (model 6). Then, when an interaction between Δ9THC: cannabidiol: cannabigerol was considered, terms for cannabidiol, cannabigerol and the cannabidiol: cannabigerol interaction were positive (model 7). When a tobacco interaction was added to this three-way cannabinoid interaction, no change was induced in the final model (model 8).

3.3. Legal Status

The levels of ASDR by legal status for the groups Illegal, Decriminalized, Medical and Legal are, respectively, 59.42 (mean, 95%C.I.; 51.71, 67.12), 55.81 (36.53, 75.10), 74.53 (54.50, 94.56), 207.00 (11.59, 402.41). These values are significantly different (ChiSqu. = 774.71, df = 711.0, p = 0.048).

When legal status is dichotomized as legal v. the remainder, the ASDR values for not legal and legal are 61.70 (54.78, 68.62) and 207.00 (11.59, 402.41), respectively (t = 1.45, df = 6.01, p = 0.19). When these data are considered as a two-by-two table, the R.R. is 2.73 (2.66, 2.80), attributable fraction in the exposed (AFE) is 3.57% (3.43%, 3.71%), Attributable Fraction in the Population (AFP) 3.55% (3.41%, 3.69%) and p = 0.0000. The E-value estimate is 4.90, and 4.76 is the lower bound of the 95% confidence interval.

When cannabis legal status is considered as a factor in a space–time regression, the results shown in Table 2 are derived, and “Legal” cannabis status is shown to be highly significant (model 1; β-est. = 0.47 (0.18, 0.75) p = 0.0013). The applicable E-values are 10.45 with 3.31 as the 95% lower bound

When dichotomized cannabis legal status is considered “Legal” cannabis status is again shown to be highly significant (β-est. = 0.34 (0.10, 0.58) p = 0.0059; model 2, Table 2). The applicable E-values are 6.47 with 2.23 as the 95% lower bound (with RR 3.51, (1.44, 8.54)).

When last month cannabis use was considered in an interactive spatial model along with cannabis legal status, the effect of cannabis legal status persisted in interactions (Table 2).

3.4. Multicollinearity Diagnostics

Belsley condition indices below conventional thresholds (10) and the absence of multiple predictors loading heavily on the same high-index dimension indicated that multicollinearity was not problematic in the geospatial models analyzed. Supplementary Table S4 presents the detailed results by covariate for each of the eight final models presented in Table 1. The Belsley collinearity diagnostics indicate no serious multicollinearity among the predictors as all condition indices are below 10, which is well under commonly accepted thresholds for concern. Variance decomposition proportions show that although the largest condition index (9.78) is associated with a relatively high proportion for CBG, no other variables share large variance decomposition proportions in that dimension, meaning no group of predictors is jointly contributing to instability. Therefore, the cannabinoid variables appear sufficiently independent for reliable regression estimation.

3.5. Robustness Analysis

The impact of the highest scoring state, Nevada, on these analyses was considered in a robustness analysis by analyzing the spatiotemporal dataset without Nevada. The results of space–time regression are shown in Supplementary Table S5, and the applicable E-values and relative risks are shown in Supplementary Table S6. Results of this analysis essentially confirms the general analysis shown in Table 1.

The ASD Rates for the 19 state–year combinations was also imputed by the last-observation-carried-forward technique. This technique confirmed the analysis and conclusions presented above.

3.6. Survey Regression

This dataset had an unusual structure since the ASDR’s were provided in quinquennial segments but the drug exposure data was available on an annual basis (Supplementary Table S1). This was addressed by the use of survey regression, with drug use data averaged within each quinquennium and the periods of analysis incorporated into the regression structure of the survey models. As shown in Table 3, these additive and interactive analyses were strongly confirmatory of the geospatial analysis but at much higher levels of statistical significance. The relative risks and E-values for all covariates from these regression summaries are shown in Supplementary Table S3.

3.7. Non-Overlapping Dataset

One issue with the above analysis is that some of the data is for overlapping time periods. Therefore, a new analysis was performed in which used the data from the 2003–2007, 2008–2012 and 2016–2020 periods were used. As shown in Supplementary Table S7 for additive and interactive models of both cannabis and cannabinoids, cannabis and cannabinoids were strongly associated with ASDR in all models.

3.8. Temporal Lagging

Another potential issue is the use of data that was not temporally lagged. Since human gestation normally takes about nine months, it may be argued that a more realistic manner in which to proceed would be to consider lagging the predictor covariates by one year to better model the biology of the exposure–outcome paradigm. The results of this analysis are presented in Supplementary Table S8 for additive and interactive models of both cannabis and cannabinoids. In each case, cannabis and cannabinoids are strongly associated with ASDR.

3.9. Case Ascertainment

It may also be hypothesized that the method of case ascertainment within each registry might contribute to the ASDR reported in that state. For this reason, case ascertainment styles were extracted from the annual NBDPN reports and tabulated against the ASDR. The results of this analysis are shown in Supplementary Figure S13 as both groups over aggregated time and jittered scatterplots across time. There were 93, 103 and 74 registries in the passive, mixed and active groups, respectively. The median (IQR) rates in the passive, mixed and active case ascertainment registries were 73.9 (26.7, 113.0); 36.2 (24.5, 68.2) and 34.9 (26.9, 56.5), and these were not found to be statistically different (ChiSqu. = 503, df = 464, p = 0.10). Thus, contrary to what might otherwise be expected, case ascertainment did not play an obvious role in systematically contributing to ASDR.

3.10. Sensitivity Analysis

E-values were calculated to assess the minimum strength of unmeasured confounding and their relationships with both the predictor and outcome covariates to explain away the observed associations, conditional on the ecological and design-based nature of the analysis. E-Values from the above regression and bivariate comparisons are summarized (Table 4). In total, 35/44 (79.5%) E-value estimates are seen to exceed 9, and all surpass the threshold of causality at 1.25; 12/44 (27.3%) lower bounds exceed 9, and all exceed 1.25. The median (IQR) of the E-values is 13.78 (9.25, 123.62), and the median (IQR) of the 95% lower confidence interval of the E-values is 3.31 (1.48, 9.53). These data demonstrate that results are robust to unmeasured confounding.

4. Discussion

4.1. Main Results

Study findings indicate an association of the ASDR and cannabis use. Indeed, in multivariate spatial models, cannabis use was the only covariate to survive model reduction. Plotting cannabis legal status across space and time phenocopied last month cannabis use, likely through the well-described effect of cannabis legalization increasing cannabis use [ref. 58,ref. 59]. This finding is concordant with earlier studies in locations including Hawaii, Colorado, Canada, Australia, the USA and Europe [ref. 2,ref. 3,ref. 4,ref. 5,ref. 6,ref. 7]. Of note, the strength of the association reported in this study was higher than in earlier reports, as revealed by the median E-values. Data were further supported and reinforced by the demonstration of closely concordant results using non-overlapping and temporally lagged spatial analyses and by analyses which omitted Nevada, the US State with the highest ASDRs. Taken together, the convergence of results across IPW-weighted survey models and spatial regression approaches, along with large E-values and associated analyses, provides convergent evidence consistent with a causal interpretation, although residual confounding cannot be completely excluded.

One important exception to this overall pattern occurred in the Midwestern states of Kentucky, Tennessee, Mississippi, and Missouri, where elevated rates of ASD were not accompanied by high rates of cannabis use in NSUDH. These States are known to grow large cannabis crops [ref. 60,ref. 61]. One speculative possibility is that cannabinoid exposure may occur through environmental contamination as an inevitable consequence of large-scale cannabis cultivation where plant refuse and debris contaminate rivers, streams and lakes and thereby the local water supply and food chain [ref. 62]. This has been documented around the Great Lakes and west coast waterways [ref. 63,ref. 64,ref. 65,ref. 66]. Demonstration of such a possibility would require further study. An alternative possibility is the existence of busy local cannabis dispensaries [ref. 67].

One of remarkable study feature is the very elevated ASDRs reported. Teratology rates of almost one in eleven are otherwise unusual in the neonatal literature. This is, however, consistent with the well-described exponential dose–response relationship seen in numerous laboratory studies of cannabinoid genotoxicity [ref. 68,ref. 69,ref. 70]. The implications of this exponential relationship for public health are far-reaching since it implies that relatively modest increases in teratogenic exposure can induce disproportionate rises in teratological outcomes. This issue of the supra-exponential rise in the ASDR has been considered formally [ref. 71]. That work also found a four-fold rise in the ASDR across the period 2003–2020.

The geospatial techniques used in this study do not allow temporal lags to be introduced. However since cardiogenesis occurs within the first three months of gestation, this technique should suffice to cover this biogenic period.

The method of state congenital anomaly registry case ascertainment did not appear to contribute systematically to ASDRs. Omitting the highest scoring state, Nevada, from the analysis did not materially impact the main conclusions.

This study reported higher ASDRs in states where cannabis was legal than elsewhere. This theme is considered in more detail in a companion paper [ref. 71].

4.2. Mechanisms

Central to any consideration of the link between defined exposures and specific outcomes is the existence of biologically plausible mechanisms that are critical to any potential causal mechanistic pathway. For this reason, a brief consideration of cardiac developmental embryology and cardiogenic epigenomics is relevant to this discussion. It is noted that this is external evidence sourced from independent investigators and does not arise from the present investigation.

During embryological development, the heart and great vessels are known to form from a complex coalescence of cells from several cardiogenic fields including the primary, secondary, and lateral heart fields, the proepicardium, the nuchal crest, and from the somites of the pharyngeal arches [ref. 72]. Atrial septal formation occurs when an initial atrial septum primum forms, degenerates later, and is superseded by the atrial septum secundum, involving complex cellular dynamics [ref. 72].

Cardiac development occurs through the sequential activation of cardiogenic gene cassettes, which are controlled by gradients of several tissue morphogens governing the development of the central cardiovasculature, along with epigenomic activation of relevant genes, both of which are widely disrupted by cannabinoids [ref. 72,ref. 73,ref. 74,ref. 75,ref. 76].

Key genes involved in cardiogenesis including sonic hedgehog (shh), VEGF, notch, Eph and Ephrin [ref. 72] have all been shown to be disrupted epigenomically by cannabis [ref. 74]. The delicate cellular choreography of cardiogenesis is controlled by gene expression and the concentration of local morphogen gradients which guide and direct cell migration and development. Sonic hedgehog (shh) is a key morphogen which is known to play a critical role in the development of the heart and other body tissues. Shh is known to be blocked directly by both Δ9THC and cannabidiol, amongst other cannabinoids [ref. 77], and epigenetically by cannabis [ref. 74]. Similarly, many other key cardiogenic morphogens [ref. 78] are inhibited both directly and epigenetically by cannabis, including fibroblast growth factor (FGF), bone morphogenetic proteins (BMP), retinoids, notch, Wnt and hippo signaling [ref. 5,ref. 17,ref. 74,ref. 75,ref. 79,ref. 80,ref. 81,ref. 82,ref. 83,ref. 84].

For example BMPs are involved in heart tube formation, formation of the cardiac outflow tract, and endocardial cushion and heart valve formation [ref. 78,ref. 85,ref. 86]. A double feedforward mechanism involving BMP4/Smad1/5, Wnt3/Tcf3 and Nodal/Smad2/3 has been described in which they cooperate to open mesoderm enhancers, affecting the expression of cardiac lineage genes [ref. 78]. Indeed, morphogens such as BMP, notch and TGF-β operate collaboratively to orchestrate cardiac outflow tract and valvular development [ref. 85].

FGF8, FGF9, FGF10, and FGF16 are involved in embryonic heart development and act as paracrine signals, influencing processes such as cardiac progenitor cell proliferation, differentiation, and patterning [ref. 87]. FGF8 is essential for the migration and differentiation of cardiac progenitor cells within the anterior heart field, which contributes to the formation of the right ventricle and outflow tract [ref. 88]. FGF10 expression is regulated by key transcription factors like ISL1, NKX2-5, and TBX1 and is crucial for the deployment of cardiogenic progenitors from the secondary heart field [ref. 89].

4.3. Causality

The Hill criteria are frequently used to bridge association to causation. Data in the present study and others establish that cannabis fulfills all of these criteria of strength of association, consistency amongst other series, temporality, external coherence, biological pathways, dose–response effects, other locations and experimentation. Key amongst these are cellular mechanisms that might give rise to the observed effects. As discussed briefly above and elsewhere, there are many mechanisms linking cannabis exposure to defective cardiogenesis [ref. 5,ref. 17,ref. 20,ref. 74,ref. 77,ref. 90].

Beyond the high E-values and low p-values reported herein, several technical features of this analysis point to a strong association between cannabis and the ASDR. Th study’s findings that cannabis was the most robust predictive covariate when both were considered covarying across space and time points to a close association. The survey regressions were inverse probability-weighted, which established a pseudo-randomized analytical paradigm in this dataset. The use of E-Values tightly constrained the potent effect demanded of any unmeasured confounding covariate. The median E-value was 13.78, with a median lower bound of 3.31. This establishes the strength of the association some unmeasured confounder must have with both the measured covariates and the outcome variable to explain the close relationship observed. It should also be noted that while E-values can be informative, they do not resolve ecological fallacy, measurement error, selection bias, or outcome misclassification problems.

It should also be noted that several interaction terms produced extremely large or infinite risk ratio estimates, resulting in correspondingly infinite E-values. Such estimates reflect sparse data or near-separation in higher-order interaction terms rather than meaningful causal effects. Thus very elevated E-values have limited practical meaning due to confounding in the context of an ecological study.

Having made these observations, the present study is strictly ecological in design. No individual-level participant data was available to the researchers. This design prevents the drawing of causal inferences in an absolute sense. Results are interpreted as population effects under measured covariate control under the standard assumptions of consistency, positivity, conditional exchangeability and correct model specification. The issue of causality is explicated in detail elsewhere [ref. 2,ref. 17,ref. 71,ref. 91]. This study can therefore be regarded as demonstrating a close association and is thus primarily hypothesis-generating.

4.4. Generalizability

Many features of the present study suggest that the main results are widely generalizable. Confirmatory results were obtained with bivariate and multivariate geospatial regression techniques, Chi squared tests and t-tests of inter-group comparisons and E-values. Thus, the analysis was very internally consistent. Elevated E-values indicate the result is robust to external confounding. External validation is confirmed by previous large studies from many sites that all draw similar conclusions [ref. 2,ref. 3,ref. 4,ref. 5,ref. 6,ref. 7].

4.5. Strengths and Limitations

This study has many strengths, including the use of the NBDPN ASD dataset, the use of space–time analyses, the confirmation of the geospatial analysis results by survey regressions, and the absence of multiple collinearity and E-values, which are one of the principle statistical tools of causal inference. Robustness, sensitivity, temporal lagged, inverse probability weighted survey regression, and non-overlapping spatiotemporal models were also examined and produced uniformly concordant results. Study weaknesses include limitation to only a single race due to the technical limitations of space–time regression, as high missing data rates limited the application of space–time methods to other ethnicities. The lack of individual-level cannabis and ASD data is a limitation shared by many other epidemiological studies. State-level associations do not necessarily reflect individual-level causal effects, and aggregate confounding (e.g., states with liberal cannabis policy may differ systematically in healthcare access, screening intensity, or reporting) introduces possible errors which cannot be ruled out regardless of E-value magnitude. Other ethnicities are addressed in a companion paper [ref. 91]. Earlier single-center retrospective studies have questioned CDC ASDR data integrity based on the finding of 64.3% accuracy in a single-center pediatric population [ref. 92]. As shown in Table 3, for that reference, the accuracy of ASD diagnosis was 19% in patients aged 41–64 years, 23% in age group 21–40 years and 64% < 21 years. This trend suggests increasing rates of diagnostic accuracy in the younger cohorts considered, and it may be reasonably expected that in the very young cohorts reported by the NBDPN, younger than 5 years old, this accuracy is significantly higher again. Whilst this issue relates to this work, it is unlikely to invalidate its main conclusions. This was an ecological study, which limits its ability to formally establish causal relationships. Some covariates were not measured in the present work, including maternal age, obesity, pregestational diabetes, exposure to teratogenic medications (antiepileptics, lithium), folic acid supplementation, environmental pollution, and access to prenatal care. These would be issues which could be taken up in subsequent research. It should also be noted that the NSDUH reflects the general population, not pregnant individuals, and state-level cannabis metrics may be proxies for multiple correlated factors (policy environment, healthcare access, socioeconomic patterns). Having said that, it has been shown several times that cannabis use by pregnant mothers closely parallels that of the general community [ref. 93,ref. 94,ref. 95,ref. 96]. Other covariates which may account for some of the marked variance in ASDRs such as changes in prenatal/postnatal imaging, screening intensity, and reporting completeness over time were not available to the present analysis but may be examined in future research. However it seems unlikely to us that such putative case ascertainment practices would co-vary across space and time coincident and concurrent with spatiotemporal cannabis distribution.

5. Conclusions

Many important conclusions stem from this study. The data demonstrate a close relationship between population cannabis and cannabinoid exposure with the current rising US ASD trend which, in some places, is rising to extremely high levels, approaching almost one in eleven births. Elevated cannabis-associated ASDRs may be an example of intergenerational teratogenicity and thus the transgenerational transmission of genotoxic effects. Several technical aspects of the present report support the strength of this association, including the close co-variation of cannabis exposure across space and time, the concordance between two different forms of multiple regression, the use of pseudo-randomization in the survey regressions, the low p-values and very high E-values, numerous cellular and molecular mechanistic pathways implicating cannabinoid genotoxicity and epigenotoxicity, obvious temporality, and concordance with many other large recent reports from other nations. Several cannabinoids may be implicated, including cannabidiol, Δ9THC, and cannabigerol, consistent with cannabinoid genotoxicity being a class effect across all cannabinoid species. Concerningly, findings in several Midwest States do not agree with local cannabis smoking trends, as documented in national surveys, but may relate to large cannabis crops grown in those vicinities, thereby potentially implicating either other exposure routes or possibly local marketing practices. This study extends earlier work by showing that the population cannabis use–ASDR association is much stronger than previously described. Overall, the findings strongly suggest that communities need to give careful consideration to the genotoxic and epigenotoxic activities of cannabinoids in designing environmentally prudent and sustainable public health policies to protect the health of current generations and the epigenetic health of generations yet to come. Further research on these important topics is strongly indicated.

References

1. 1. National Birth Defects Prevention Network Major Birth Defects Data from Populatoin-Based Birth Defects Surveillance Programs in the United States, 2016–2020 National Birth Defects Prevention Network National Birth Defects Prevention Network, NBDPNHouston, TX, USA2024Volume 2020
2. A.S. Reece, G.K. Hulse. Contemporary epidemiology of rising atrial septal defect trends across USA 1991–2016: A combined ecological geospatiotemporal and causal inferential study. BMC Pediatr., 2020. [DOI | PubMed]
3. M.B. Forrester, R.D. Merz. Risk of selected birth defects with prenatal illicit drug use, Hawaii, 1986–2002. J. Toxicol. Environ. Health, 2007. [DOI]
4. A.S. Reece, G.K. Hulse. Cannabis Teratology Explains Current Patterns of Coloradan Congenital Defects: The Contribution of Increased Cannabinoid Exposure to Rising Teratological Trends. Clin. Pediatr., 2019. [DOI]
5. A.S. Reece, G.K. Hulse. European Epidemiological Patterns of Cannabis- and Substance-Related Congenital Cardiovascular Anomalies: Geospatiotemporal and Causal Inferential Study. Environ. Epigenetics, 2022. [DOI]
6. A.S. Reece, G.K. Hulse. Geotemporospatial and causal inference epidemiological analysis of US survey and overview of cannabis, cannabidiol and cannabinoid genotoxicity in relation to congenital anomalies 2001–2015. BMC Pediatr., 2022. [DOI | PubMed]
7. A.S. Reece, G.K. Hulse. Canadian Cannabis Consumption and Patterns of Congenital Anomalies: An Ecological Geospatial Analysis. J. Addict. Med., 2020. [DOI | PubMed]
8. 8. Colorado: Department of Public Health and the Environment Colorado Responds to Children with Special Needs—Birth Defect Data, Colorado1st ed. Colorado: Department of Public Health and the Environment Colorado: Department of Public Health and the EnvironmentDenver, CO, USA2018Volume 1Available online: https://cohealthviz.dphe.state.co.us/t/HealthInformaticsPublic/views/COHIDBirthDefectsdashboard_15614030383960/BirthDefectCasesStatistics?iframeSizedToWindow=true&%3Aembed=y&%3AshowAppBanner=fals
9. T.J. VanderWeele, P. Ding, M. Mathur. Technical Considerations in the Use of the E-Value. J. Causal Inference, 2019. [DOI]
10. A. Moore, J. Roulean, E. Skarsgard. Chapter 7 Gastroschisis. Congenital Anomalies in Canada, 2013: A Perinatal Health Surveillance Report, 2013
11. A.S. Reece, G.K. Hulse. Cannabinoid- and Substance-Relationships of European Congenital Anomaly Patterns: A Space-Time Panel Regression and Causal Inferential Study. Environ. Epigenetics, 2022. [DOI | PubMed]
12. A.S. Reece, G.K. Hulse. Broad Spectrum epidemiological contribution of cannabis and other substances to the teratological profile of northern New South Wales: Geospatial and causal inference analysis. BMC Pharmacol. Toxicol., 2020. [DOI | PubMed]
13. S. Howell, T. Endo, S. MacLeod, S. Cornes. Congenital Anomalies in Queensland: 1 July 2007 to 30 June 2010. Stat. Anal. Rep., 2011
14. Congenital Heart Defects (CHDs)
15. W.F. Geber, L.C. Schramm. Teratogenicity of marihuana extract as influenced by plant origin and seasonal variation. Arch. Int. Pharmacodyn. Ther., 1969. [PubMed]
16. W.F. Geber, L.C. Schramm. Effect of marihuana extract on fetal hamsters and rabbits. Toxicol. Appl. Pharmacol., 1969. [DOI | PubMed]
17. A.S. Reece, G.K. Hulse. Chapter 3: Geospatiotemporal and Causal Inferential Analysis of United States Congenital Anomalies as a Function of Multiple Cannabinoid- and Substance-Exposures: Phenocopying Thalidomide and Hundred Megabase-Scale Genotoxicity. Epidemiology of Cannabis: Genotoxicity and Neurotoxicity, Epigenomics and Aging, 2025
18. S.N. Grobner, B.C. Worst, J. Weischenfeldt, I. Buchhalter, K. Kleinheinz, V.A. Rudneva, P.D. Johann, G.P. Balasubramanian, M. Segura-Wang, S. Brabetz. The landscape of genomic alterations across childhood cancers. Nature, 2018. [DOI | PubMed]
19. X. Ma, Y. Liu, Y. Liu, L.B. Alexandrov, M.N. Edmonson, C. Gawad, X. Zhou, Y. Li, M.C. Rusch, J. Easton. Pan-cancer genome and transcriptome analyses of 1699 paediatric leukaemias and solid tumours. Nature, 2018. [DOI | PubMed]
20. A.S. Reece, G.K. Hulse. Chapter 5: Multivalent Cannabinoid Epigenotoxicities and Multigenerational Aging. Epidemiology of Cannabis: Genotoxicity and Neurotoxicity, Epigenomics and Aging, 2025
21. A.S. Reece, A. Norman, G.K. Hulse. Cannabis exposure as an interactive cardiovascular risk factor and accelerant of organismal ageing: A longitudinal study. BMJ Open, 2016. [DOI]
22. K.T. Phillips, K.L. Pedula, N.G. Choi, K.K. Tawara, V. Simiola, D.D. Satre, A. Owen-Smith, F.F. Lynch, J. Dickerson. Chronic health conditions, acute health events, and healthcare utilization among adults over age 50 in Hawai’i who use cannabis: A matched cohort study. Drug Alcohol Depend., 2022. [DOI | PubMed]
23. A.S. Reece, G.K. Hulse. Impacts of cannabinoid epigenetics on human development: Reflections on Murphy et. al. ‘cannabinoid exposure and altered DNA methylation in rat and human sperm’ epigenetics 2018; 13: 1208–1221. Epigenetics, 2019. [DOI | PubMed]
24. J.P. Allen, J.S. Danoff, M.A. Costello, G.L. Hunt, A.F. Hellwig, K.M. Krol, S.G. Gregory, S.N. Giamberardino, K. Sugden, J.J. Connelly. Lifetime marijuana use and epigenetic age acceleration: A 17-year prospective examination. Drug Alcohol Depend., 2022. [DOI | PubMed]
25. A.S. Reece, G.K. Hulse. Impact of Converging Sociocultural and Substance-Related Trends on US Autism Rates: Combined Geospatiotemporal and Causal Inferential Analysis. Eur. Arch. Psychiatry Clin. Neurosci., 2022. [DOI]
26. A.S. Reece, G.K. Hulse. Chapter 2: Linked Rise of Cannabis Use and Autism Incidence Demonstrated by Close Three-Level Geospatiotemporal Relationships, USA, 1990–2011. Epidemiology of Cannabis: Genotoxicity and Neurotoxicity, Epigenomics and Aging, 2025
27. S.E. Paul, A.S. Hatoum, J.D. Fine, E.C. Johnson, I. Hansen, N.R. Karcher, A.L. Moreau, E. Bondy, Y. Qu, E.B. Carter. Associations Between Prenatal Cannabis Exposure and Childhood Outcomes: Results from the ABCD Study. JAMA Psychiatry, 2021. [DOI | PubMed]
28. D.A.A. Baranger, S.E. Paul, S.M.C. Colbert, N.R. Karcher, E.C. Johnson, A.S. Hatoum, R. Bogdan. Association of Mental Health Burden with Prenatal Cannabis Exposure From Childhood to Early Adolescence: Longitudinal Findings from the Adolescent Brain Cognitive Development (ABCD) Study. JAMA Pediatr., 2022. [DOI | PubMed]
29. L. Lepow, A. Wagner, S. Peri, F. Adams, S.A. Ramakrishnan, M.A. Alam, R.B. Shaik, N.A. Hubbard, H.W. Koenigsberg, Y. Hurd. Independent and Interactive Impacts of Prenatal Exposure to Legal Substances and Childhood Trauma on Emotion Processing in Pre-Adolescents: Preliminary Findings From the ABCD Study. JAACAP Open, 2024. [DOI | PubMed]
30. National Survey of Drug Use and Health (NSDUH)
31. Cannabis Overview
32. L.O. Gostin, J.G. Hodge, S.A. Wetter. Enforcing Federal Drug Laws in States Where Medical Marijuana Is Lawful. J. Am. Med. Assoc., 2018. [DOI]
33. S. Chandra, M.M. Radwan, C.G. Majumdar, J.C. Church, T.P. Freeman, M.A. ElSohly. New trends in cannabis potency in USA and Europe during the last decade (2008–2017). Eur. Arch. Psychiatry Clin. Neurosci., 2019. [DOI | PubMed]
34. M.A. ElSohly, S. Chandra, M. Radwan, C.G. Majumdar, J.C. Church. A Comprehensive Review of Cannabis Potency in the United States in the Last Decade. Biol. Psychiatry Cogn. Neurosci. Neuroimaging, 2021. [DOI | PubMed]
35. M.A. ElSohly, Z. Mehmedic, S. Foster, C. Gon, S. Chandra, J.C. Church. Changes in Cannabis Potency Over the Last 2 Decades (1995–2014): Analysis of Current Data in the United States. Biol. Psychiatry, 2016. [DOI | PubMed]
36. M.A. ElSohly, S.A. Ross, Z. Mehmedic, R. Arafat, B. Yi, B.F. Banahan. Potency trends of delta9-THC and other cannabinoids in confiscated marijuana from 1980–1997. J. Forensic Sci., 2000. [DOI | PubMed]
37. T.P. Freeman, S. Craft, J. Wilson, S. Stylianou, M. ElSohly, M. Di Forti, M.T. Lynskey. Changes in delta-9-tetrahydrocannabinol (THC) and cannabidiol (CBD) concentrations in cannabis over time: Systematic review and meta-analysis. Addiction, 2021. [DOI | PubMed]
38. A.S. Reece, G.K. Hulse. Epidemiological association of cannabinoid- and drug-exposures and sociodemographic factors with limb reduction defects across USA 1989–2016: A geotemporospatial study. Spat. Spatio-Temporal Epidemiol., 2022. [DOI]
39. A.S. Reece, G.K. Hulse. Epidemiological Associations of Various Substances and Multiple Cannabinoids with Autism in USA. Clin. Pediatr. Open Access, 2019. [DOI]
40. A.S. Reece, G.K. Hulse. Epidemiology of Δ8THC-Related Carcinogenesis in USA: A Panel Regression and Causal Inferential Study. Int. J. Environ. Res. Public Health, 2022. [DOI | PubMed]
41. H. Wickham, M. Averick, J. Bryan, W. Chang, L.D. McGowan, R. Francios, G. Groelmund, A. Hayes, L. Henry, J. Hester. Welcome to the Tidyverse. J. Open Source Softw., 2019. [DOI]
42. E. Pebesma. Simple Features for R: Standardized Suport for Spatial Vector Data. R J., 2018. [DOI]
43. R. Bivand, L. Anselin, O. Berke, A. Bernat, M. Carvalho, Y. Chun, C. Dormann, S. Dray, R. Halbersma, N. Lewis-Koh. The spdep Package. Compr. R Arch. Netw. Version, 2007
44. R. Bivand, G. Millo, G. Piras. A Review of Software for Spatial Econometrics in R. Mathematics, 2021. [DOI]
45. G. Millo, G. Piras. Package ‘splm’. Compr. R Arch. Netw. Version, 2018
46. G. Millo, G. Piras. splm: Spatial Panel Data Models in R. J. Stast. Softw., 2012
47. M.J. Crawley. The R Book, 2007
48. Viridis: Default Color Maps from ‘Matplotlib’
49. Colorplaner: Ggplot2 Extension to Visualize Two Variables Per Color Aesthetic Through Colorspace Projection
50. T.J. VanderWeele, M.B. Mathur. Commentary: Developing best-practice guidelines for the reporting of E-values. Int. J. Epidemiol., 2020. [DOI | PubMed]
51. T.J. Vander Weele, P. Ding. Sensitivity Analysis in Observational Research: Introducing the E-Value. Ann. Intern. Med., 2017. [DOI | PubMed]
52. M. Stevenson, E. Sergeant. epiR: Tools for the Analysis of Epidemiological Data. 2025
53. A. Zeileis, G. Grothendieck. Zoo: S3 Infrastructure for Regular and Irregular Time Series. J. Stat. Softw., 2005. [DOI]
54. “Survey: Analysis of Complex Survey Samples”. R Package Version 4.4
55. W.M. Van der Wal, R.B. Geskus. Ipw: An R Package for Inverse Probabilty Weighting. J. Stat. Softw., 2011
56. D.A. Belsley, E. Kuh, R.E. Welsch. Regression Diagnostics: Identifying Influential Data and Sources of Collinearity, 1980
57. J.P. LeSage, R.K. Pace. Introduction to Spatial Econometrics, 2009
58. D.T. Myran, M. Pugliese, L.D. Harrison, M. Solmi, K.K. Anderson, J.G. Fiedorowicz, Y. Finkelstein, D. Manuel, M. Taljaard, C. Webber. Changes in Incident Schizophrenia Diagnoses Associated with Cannabis Use Disorder After Cannabis Legalization. JAMA Netw. Open, 2025. [DOI | PubMed]
59. D.T. Myran, M. Pugliese, A.J. McDonald, J. Xiao, B. Fischer, Y. Finkelstein, P. Tanuseputro, J. Firth, A. Pakpour, C.W. Hsu. Cannabis Use Disorder Emergency Department Visits and Hospitalizations and 5-Year Mortality. JAMA Netw. Open, 2025. [DOI | PubMed]
60. How Missouri Became a Cannabis Mecca
61. J.J. Coleman. Did the 2018 Farm Bill’s Hemp Provisions Decriminalize Marijuana?. J. Drug Policy Anal., 2021. [DOI]
62. B. Wagner, P. Gerletti, P. Fürst, O. Keuth, T. Bernsmann, A. Martin, B. Schäfer, J. Numata, M.C. Lorenzen, R. Pieper. Transfer of cannabinoids into the milk of dairy cows fed with industrial hemp could lead to Δ9-THC exposure that exceeds acute reference dose. Nat. Food, 2022. [DOI | PubMed]
63. M. Parolini, S. Castiglioni, S. Magni, C. Della Torre, A. Binelli. Increase in cannabis use may indirectly affect the health status of a freshwater species. Environ. Toxicol. Chem., 2017. [DOI | PubMed]
64. M. Parolini, A. Binelli. Oxidative and genetic responses induced by Delta-9-tetrahydrocannabinol (Delta-9-THC) to Dreissena polymorpha. Sci. Total Environ., 2014. [DOI]
65. State Board—Cannabis Cultivation Water Quality
66. Water Boards Cannabis Cultivation
67. B.K. Madras, P.J. Larkin. Rescheduling Cannabis-Medicine or Politics?. JAMA Psychiatry, 2025. [DOI | PubMed]
68. F.W. Busch, D.A. Seid, E.T. Wei. Mutagenic activity of marihuana smoke condensates. Cancer Lett., 1979. [DOI | PubMed]
69. A.M. Zimmerman, A.Y. Raj. Influence of Cannabinoids on Somatic Cells in vivo. Pharmacology, 1980. [DOI | PubMed]
70. S.K. Tahir, A.M. Zimmerman. Influence of marihuana on cellular structures and biochemical activities. Pharmacol. Biochem. Behav., 1991. [DOI | PubMed]
71. A.S. Reece, G.K. Hulse. Impact of Cannabis and Cannabis Legalization on US Atrial Septal Defect Rates. J. Xenobiot., 2026. [DOI | PubMed]
72. B.M. Carlson. Human Embryology and Developmental Biology, 2019
73. C.T. Watson, H. Szutorisz, P. Garg, Q. Martin, J.A. Landry, A.J. Sharp, Y.L. Hurd. Genome-Wide DNA Methylation Profiling Reveals Epigenetic Changes in the Rat Nucleus Accumbens Associated with Cross-Generational Effects of Adolescent THC Exposure. Neuropsychopharmacology, 2015. [DOI | PubMed]
74. R. Schrott, S.K. Murphy, J.L. Modliszewski, D.E. King, B. Hill, N. Itchon-Ramos, D. Raburn, T. Price, E.D. Levin, R. Vandrey. Refraining from use diminishes cannabis-associated epigenetic changes in human sperm. Environ. Epigenet., 2021. [DOI | PubMed]
75. S.K. Murphy, N. Itchon-Ramos, Z. Visco, Z. Huang, C. Grenier, R. Schrott, K. Acharya, M.H. Boudreau, T.M. Price, D.J. Raburn. Cannabinoid exposure and altered DNA methylation in rat and human sperm. Epigenetics, 2018. [DOI | PubMed]
76. E.L. Dempster, C.C.Y. Wong, J. Burrage, E. Hannon, D. Quattrone, G. Trotta, V. Rodriguez, L. Alameda, E. Spinazzola, G. Tripoli. Methylomic signature of current cannabis use in two first-episode psychosis cohorts. Mol. Psychiatry, 2025. [DOI | PubMed]
77. E.W. Fish, L.B. Murdaugh, C. Zhang, K.E. Boschen, O. Boa-Amponsem, H.N. Mendoza-Romero, M. Tarpley, L. Chdid, S. Mukhopadhyay, G.J. Cole. Cannabinoids Exacerbate Alcohol Teratogenesis by a CB1-Hedgehog Interaction. Sci. Rep., 2019. [DOI | PubMed]
78. P. Tsaytler, J. Liu, G. Blaess, D. Schifferl, J.V. Veenvliet, L. Wittler, B. Timmermann, B.G. Herrmann, F. Koch. BMP4 triggers regulatory circuits specifying the cardiac mesoderm lineage. Development, 2023. [DOI | PubMed]
79. O. Asimaki, G. Leondaritis, G. Lois, N. Sakellaridis, D. Mangoura. Cannabinoid 1 receptor-dependent transactivation of fibroblast growth factor receptor 1 emanates from lipid rafts and amplifies extracellular signal-regulated kinase 1/2 activation in embryonic cortical neurons. J. Neurochem., 2011. [DOI | PubMed]
80. D. Richard, F. Picard. Brown fat biology and thermogenesis. Front. Biosci., 2011. [DOI]
81. D. Fraher, M.K. Ellis, S. Morrison, S.L. McGee, A.C. Ward, K. Walder, Y. Gibert. Lipid Abundance in Zebrafish Embryos Is Regulated by Complementary Actions of the Endocannabinoid System and Retinoic Acid Pathway. Endocrinology, 2015. [DOI | PubMed]
82. R. Tanveer, A. Gowran, J. Noonan, S.E. Keating, A.G. Bowie, V.A. Campbell. The endocannabinoid, anandamide, augments Notch-1 signaling in cultured cortical neurons exposed to amyloid-beta and in the cortex of aged rats. J. Biol. Chem., 2012. [DOI | PubMed]
83. A. Vallée, Y. Lecarpentier, R. Guillevin, J.N. Vallée. Effects of cannabidiol interactions with Wnt/β-catenin pathway and PPARγ on oxidative stress and neuroinflammation in Alzheimer’s disease. Acta Biochim. Biophys. Sin., 2017. [DOI | PubMed]
84. J. Petko, T. Tranchina, G. Patel, R. Levenson, S. Justice-Bitner. Identifying novel members of the Wntless interactome through genetic and candidate gene approaches. Brain Res. Bull., 2018. [DOI | PubMed]
85. E.C. Délot. Control of endocardial cushion and cardiac valve maturation by BMP signaling pathways. Mol. Genet. Metab., 2003. [DOI | PubMed]
86. J. Rivera-Feliciano, C.J. Tabin. Bmp2 instructs cardiac progenitors to form the heart-valve-inducing field. Dev. Biol., 2006. [DOI | PubMed]
87. G.S. Brown, J. Jang, D. Li. Growth factors and their roles in cardiac development and regeneration: A narrative review. Pediatr. Med., 2022. [DOI]
88. R. Ilagan, R. Abu-Issa, D. Brown, Y.-P. Yang, K. Jiao, R.J. Schwartz, J. Klingensmith, E.N. Meyers. Fgf8 is required for anterior heart field development. Development, 2006. [DOI | PubMed]
89. F. Hubert, S.M. Payan, F. Rochais. FGF10 Signaling in Heart Development, Homeostasis, Disease and Repair. Front. Genet., 2018. [DOI | PubMed]
90. A.S. Reece, G.K. Hulse. Epigenomic and Other Evidence for Cannabis-Induced Aging Contextualized in a Synthetic Epidemiologic Overview of Cannabinoid-Related Teratogenesis and Cannabinoid-Related Carcinogenesis. Int. J. Environ. Res. Public Health, 2022. [DOI | PubMed]
91. A.S. Reece, G.K. Hulse. Cannabinoid Impacts on Ethnic Modulation of Atrial Septal Defect Prevalence USA. J. Perinatol., 2026
92. F.H. Rodriguez, G. Ephrem, J.F. Gerardin, C. Raskind-Hood, C. Hogue, W. Book. The 745.5 issue in code-based, adult congenital heart disease population studies: Relevance to current and future ICD-9-CM and ICD-10-CM studies. Congenit. Heart Dis., 2018. [DOI | PubMed]
93. N.D. Volkow, W.M. Compton, E.M. Wargo. The Risks of Marijuana Use During Pregnancy. J. Am. Med. Assoc., 2017. [DOI | PubMed]
94. N.D. Volkow, B. Han, W.M. Compton, C. Blanco. Marijuana Use During Stages of Pregnancy in the United States. Ann. Intern. Med., 2017. [DOI | PubMed]
95. K.C. Young-Wolff, G.T. Ray, S.E. Alexeeff, S.R. Adams, M.B. Does, D. Ansley, L.A. Avalos. Rates of Prenatal Cannabis Use Among Pregnant Women Before and During the COVID-19 Pandemic. J. Am. Med. Assoc., 2021. [DOI | PubMed]
96. K.C. Young-Wolff, L. Tucker, S. Alexeeff, M.A. Armstrong, A. Conway, C. Weisner, N. Goler. Trends in self-reported and biochemically tested marijuana use among pregnant females in california from 2009–2016. J. Am. Med. Assoc., 2017. [DOI]