Congenital anomaly epidemiological correlates of Δ8THC across USA 2003-16: panel regression and causal inferential study.

Δ8-Tetrahydrocannabinol (Δ8THC) is marketed in many US states as 'legal weed'. Concerns exist relating to class-wide genotoxic cannabinoid effects. We conducted an epidemiological investigation of Δ8THC-related genotoxicity expressed as 57 congenital anomaly (CA) rates (CARs) in the USA. CARs were taken from the Centers for Disease Control, Atlanta, Georgia. Drug exposure data were taken from the National Survey of Drug Use and Health, with a response rate of 74.1%. Ethnicity and income data were taken from the US Census Bureau. National cannabinoid exposure was taken from Drug Enforcement Agency publications and multiplied by state cannabis use data to derive state-based estimates of Δ8THC exposure. At bivariate continuous analysis, Δ8THC was associated with 23 CAs on raw CA rates, 33 CARs after correction for early termination for anomaly estimates and 41 on a categorical analysis comparing the highest and lowest exposure quintiles. At inverse probability weighted multivariable additive and interactive models lagged to 0, 2 and 4 years, Δ8THC was linked with 39, 8, 4 and 9 CAs. Chromosomal, cardiovascular, gastrointestinal, genitourinary, limb, central nervous system (CNS) and face systems were particularly affected. The minimum E-values ranged to infinity. Both the number of anomalies implicated and the effect sizes demonstrated were much greater for Δ8THC than for tobacco and alcohol combined. Δ8THC appears epidemiologically to be more strongly associated with many CAs than for tobacco and alcohol and is consistent with a cannabinoid class genotoxic/epigenotoxic effect. Quantitative causality criteria were fulfilled, and causal relationships either for Δ8THC or for cannabinoid/s, for which it is a surrogate marker, may be in operation.

Introduction

Δ8-Tetrahydrocannabinol (Δ8THC) is increasingly being marketed in many parts of USA often under the rubric of ‘legal weed’ [1]. Whilst its effects are generally said to be milder than those of Δ9THC [2], its effects in overdose are virtually identical, with reports of sedation, vomiting, drowsiness, confusion, erratic behaviour and psychotic episodes becoming so common that both the Centers for Disease Control (CDC) and the Food and Drug Administration (FDA) have recently issued cautionary warnings related to its use [3, 4].

Δ8THC is made artificially from Δ9THC usually by acidification and isomerization, so the double bond on carbon 9 of the C-ring shifts across to the eighth carbon atom on this ring [2]. Δ8THC is commonly extracted from low THC hemp residues remaining after hemp processing. In this sense, therefore, the new wave of Δ8THC poisonings coming to the attention of the CDC and the FDA represents a nascent but man-made epidemic, which is therefore fundamentally different from the previous popularization of cannabis use to sweep USA.

One of the more concerning effects of cannabinoid compounds is their genotoxicity, which has been a focus of research efforts for several decades. As long ago as 1969, severe effects of prenatal cannabinoid exposure on a range of birth defects were described, including exomphalos, foetal loss and resorption, encephalocele, spina bifida and phocomelia, in rabbits and rodents [5-7]. This report was later confirmed by a report that 21 congenital anomalies (CAs) were elevated in humans after prenatal cannabis exposure in Hawaii [8], and this has since been confirmed by similar studies from Canada, Australia, Colorado and USA [9-14]. The commonest cancer of early childhood acute lymphoid leukaemia and total pediatric cancers were both shown recently to be related to cannabis, Δ9THC and cannabigerol, thereby strongly implicating intergenerational transmission of inheritable genotoxicity and/or epigenotoxicity [15-18]. Indeed, all five chromosomal anomalies followed serially by the CDC were recently noted to be elevated causally in association with THC and cannabigerol along with several common cancers [19].

Cannabidiol and cannabichromene have been noted to be directly toxic to the bases of DNA, causing mutagenic base oxidation, and to chromosomes, causing single- and double-stranded breaks at low doses [20]. Cannabinol, cannabidiol and Δ9THC were shown long ago to be toxic to chromosomes, inducing the formation of minute fragments and ring and chain translocations [21, 22]. They also induce major structural damage to sperms and the reproductive tract in both sexes [23, 24], and cannabis smoke has been shown to induce a massive 20% oocyte loss with just the first meiotic cell division [25]. Hence, there is evidence for the genotoxic effects of multiple cannabinoids. However, rather than this being a particular property of individual cannabinoids, it was shown long ago that the genotoxic moiety lay in the olivetol nucleus, which forms the core structure of the C-ring of the cannabinoid class [26, 27].

Many in vitro studies clearly demonstrate that cannabinoid cellular toxicity exponentiates in the micromolar dose range [28-32] and it accumulates in fat stores, particularly where they are held behind membrane tight junctions as in the blood–brain barrier and the blood–testis barrier [33]. Moreover, this exponential dose-related cumulative toxicity has since been confirmed in epidemiological studies [9, 14, 15]. This implies that far from merely academic interest, these in vitro dose-ranging effects carry far-reaching implications for patterns of disease and genomic dysfunction at the level of population health.

Prenatal cannabis exposure has also been shown to have a marked effect on the epigenome through large-scale and heritable perturbation of DNA methylation, an effect expressed in increased susceptibility to addictive behaviours in later (rodent) life with analogous changes in the human sperm methylome [34-40]. Δ9THC has also been shown to suppress histone synthesis and acetylation and phosphorylation [41-43]. Reduced histone formation necessarily opens the chromatin for unregulated transcription and is known to predispose to oncogenesis.

Clearly, the prospect of multigeneration inheritance of genotoxic and epigenotoxic mutational burdens is a very serious health issue. However, cancer and CAs are together relatively rare health outcomes. Arguably, most concerning of all is the recent formal demonstration by a group from Harvard Medical School that the epigenome directly controls the expression of ageing gene networks and cassettes and can bidirectionally and powerfully impact cellular age [44]. This implies that food chain contamination by cannabinoids [45] can lead to subtle and lasting effects which are disadvantageous to the health of whole populations and their progeny. The pro-ageing effects of cannabis on biological age in humans have previously been demonstrated [46].

Altogether, these several features make a formal investigation of the potential genotoxicity/epigenotoxicity of the novel artificial cannabinoid Δ8THC a high priority. The present paper engages this research in CAs of childhood. The hypothesis that Δ8THC exposure might be linked with CAs of childhood was formulated prior to study commencement.

Methods

Data

Data on the rates of 62 CAs in each state of the USA were taken from published reports of the National Birth Defects Prevention Network (NBDPN) sponsored by the CDC, Atlanta, Georgia [47]. These reports are a collation of usually 5-year incidence rate data averaged across all races from 2003–07 to 2012–16, which are the most recent reports available at the time of writing. Data are based on reports to state registries made after birth and are based on ICD-9 codes. For each report, the mid-year of the reporting period was taken as the reference year for the purposes of this analysis. Data for Montana, Pennsylvania, Wyoming, Alabama, Idaho and South Dakota were not available. The nomenclature for orofacial clefts (cleft lip alone, cleft lip with and without cleft palate, cleft lip and cleft palate and cleft palate alone) changed during the study period, so all four CAs have been combined together for analytical purposes as orofacial clefts. Similarly, the nomenclature for arm, leg and limb reduction changed across the reporting period, so these were all combined into a single category of limb reduction defects for analysis. These adjustments left 57 CAs for analysis. At present. NBDPN reports are based on 45 anomalies, and these have changed across time.

Average annual drug use data were taken from the Restricted Data Access Scheme of the National Survey of Drug Use and Health conducted annually by the Substance Abuse and Mental Health Services Administration [48]. This is an annual survey conducted nationally, which is well established to be the representative of the non-institutionalized US adult population. The survey has a response rate of 74.1% [48]. The substances of interest were last-month cigarettes use, last-month alcohol abuse or dependence [denoted as Alcohol Use Disorder (AUD)], last-month cannabis use and last-year non-medical use of narcotic analgesics or cocaine. Ethnicity and income data were taken from the US Census Bureau using the R package tidycensus [49]. The concentration of cannabinoids in Federal seizures was taken from published reports from the Drug Enforcement Agency [50, 51]. The main cannabinoid of interest was state-level Δ8THC exposure.

Derived Data

Quintiles of substance use were calculated across the whole period by dividing it into five equal intervals, with the lowest quintile being designated as Quintile 1 and the highest quintile denoted as Quintile 5. State-level cannabinoid exposure was estimated as the product of the amount of cannabis used in each state last month and the concentration of that cannabinoid in Federal seizures by methodology established previously [15, 16, 19, 52, 53].

Statistics

Data were processed in R Studio version 1.4.1717 based on R version 4.1.1 from the Comprehensive R Archive Network. Data were manipulated using dplyr from the tidyverse suite [54]. Maps were drawn using sf [55]. Graphs were drawn using ggplot2 also from tidyverse [54], and multiple graphs were arranged using ggpubr [56]. The ggrepel package was used to position labels on volcano plots. Bivariate linear regression was conducted in Base R. Panel multivariable regression was conducted using the plm package using the two-way (simultaneous space–time) effect [57]. Overall or marginal effects of covariates in interactive multivariable models were calculated using the margins package in units of standard deviations of the dependent variable, which was the defect-specific CA rate (CAR) [58]. Categorical analyses were conducted using the epiR package [59]. Multiple models for each CA were considered simultaneously in a single analytical run using the purrr package from tidyverse to facilitate recursive modelling together with functions from the expected value (E-value), margins, broom [60, 61], epiR and dplyr packages. The analysis was performed in October 2021. All t-tests were two-tailed. P < 0.05 was considered significant.

Causal inference

All panel models were inverse probability weighted for Δ8THC exposure across all other substance exposures, a methodology which has the effect of evening out substance exposures across all groups, creating a pseudorandomized paradigm from which formal causal relationships can properly be assessed. E-values were also calculated for bivariate and multivariable regression equations using the R package E-value [62]. This value shows the degree of association required of some hypothetical unmeasured confounder variable with both the exposure of interest and the outcome of concern in order to explain away the apparently causal nature of the observed association. Inverse probability weights were calculated with the R package ipw and were utilized in all multivariable panel models [63].

Results

Table S1 provides an overview of the CAs studied and median household income and ethnicity covariate details. In total, 57 CAs were studied after adjusting facial and limb anomalies as described in the Methods section. An important determinant of the total CAR is the early termination for anomaly (ETOPFA) rate, which varies by CA. Table S2 shows a selection of surveys of ETOPFA rates for selected CAs from Australian and American series [64-71]. The final ETOPFA rate is the composite value shown in the right-hand column. It is likely that the ETOPFA rate has changed over time. The only detailed longitudinal time series we were able to locate is from the Western Australian Registry of Developmental Anomalies for Down’s syndrome, which is shown in Table S3 and runs from 1980 to 2015. The fraction of maximal termination rate in the right-hand column was applied to all CAs by year as an estimate of the temporal effect of rising ETOPFA rates in order to calculate estimates of the ETOPFA-corrected CARs (ETOPFACARs).

Figure S1 shows the concentration of various cannabinoids reported in assays of Federal seizures across the study period. Many cannabinoids including Δ8THC are noted to have risen. Figure S2 represents map-graphically estimates of the Δ8THC exposures across the country in this period.

Continuous Bivariate Analysis

Figure S3 shows the unadjusted rate of the various anomalies as a function of tobacco exposure for the raw reported CARs. Many CAs in the first three lines of this graph appear to demonstrate a positive relationship with tobacco exposure. When the exercise is repeated for AUD exposure, some of the CAs in the first line appear to be positively related (Fig. S4). When the exercise is repeated for Δ9THC, most of the CAs in the first five lines are noted to be positively related (Fig. S5). When this is done for Δ8THC, the CAs in the first three to four lines demonstrate positive regression slopes (Fig. S6). When this is done for analgesics and cocaine exposures, the results are more modest, with 3.5 and 1.5 lines of CAs demonstrating positive regression slopes (Figs S7 and S8).

The regression results of these efigures are summarized formally by recursive linear models and are presented for these substances and others in Tables S4–S11. These tables show that for cigarettes, AUD, cannabis, Δ9THC, Δ8THC, cannabidiol, analgesics and cocaine, 14, 6, 29, 31, 23, 11, 13 and 8 CAs demonstrated elevated minimum E-values (mEVs), respectively.

This exercise is then repeated for the ETOPFACARs. These results are presented graphically in Figs S9 and S10, Figs 1 and 2, Figs S11 and S12 and Tables S12–S17 for cigarettes, AUD, Δ9THC, Δ8THC, analgesics and cocaine, respectively. The results show that each substance is associated with 12, 3, 36, 33, 13 and 6 ETOPFACARs, respectively. This result makes it clear that the associations of ETOPFACARs with cannabinoids are much stronger than for other widely recognized clinical teratogens.

Trends of ETOPFA-corrected CARs as a function of Δ9THC exposure

Trends of ETOPFA-corrected CARs as a function of Δ8THC exposure

Since the above tables present rather a lot of information, these data for significant CAs are summarized in volcano plots designed by analogy with plots common in the genomic and epigenomic literature for charting the significance level of gene expression changes against its fold change. A similar exercise has been done here for significantly positively related CAs with the charting of the (negative) log of the P-value (on the y-axis) against the mEV (on the x-axis). Figure S13 performs this function for the 12 tobacco-related ETOPFACAR CAs identified in Table S12. Similarly Fig. S14 does this for three AUD-related ETOPFACARs identified in Table S13.

Figure S15 is a similar volcano plot for the 36 ETOPFACARs noted to be linked with Δ9THC exposure. Figure 3 performs a similar function for the 33 ETOPFACARs linked with Δ8THC. Several features in this figure are noteworthy. There are clearly many more CAs on this plot than for both cigarettes and AUD combined. Moreover, the scale of this plot is much greater than that of the earlier plots with log mEVs of >200 and negative log P-values of >60. Four of the five chromosomal anomalies (trisomies 21, 18 and 13 and Turner syndrome) are in or near the top-right quadrant. The fifth chromosomal anomaly also features on this plot (Deletion 22aq11.2), which although it has a P-value of 0.0034 and mEV of 5.68 × 1014 (Table S15) features lower in the plot.

Figure 3:

Volcano plot of negative log of P-value against log of minimal E-value for CAs significantly positively related to Δ8THC on bivariate testing

Finally, 11 ETOPFACARs which are similarly raised for cannabidiol are shown in Fig. S16. Again the magnitude of the represented effects is much greater for cannabidiol than for tobacco or alcohol.

Categorical Bivariate Analysis

Substance exposure metrics can be quantized into categorical quintiles of exposure in order to calculate key epidemiological indices, such as relative risk (RR), attributable fraction among the exposed (AFE) and population attributable risk (PAR). Panelled boxplots representing comparisons of the extreme quintiles are shown for tobacco, AUD, cannabis, Δ9THC, Δ8THC, analgesics and cocaine in Figs S17–S22 and Fig. 4, respectively. Reports of the RR, AFE, PAR and chi-squared tests including P-values and mEVs are shown in tabular form in Tables S18–S23 where 18, 3, 34, 41, 17 and 7 CAs are found to occur with increased mEV in the highest quintiles.

Figure 4:

Categorical panelled boxplot contrasting the ETOPFA-adjusted CAR by highest and lowest exposure quintiles of Δ8THC exposure

These results are summarized across substances in Table S24. This table shows the number of CAs having elevated RRs, AFEs and PARs and their lower bounds, along with summary metrics for the total of the negative P-value exponents, a count of the number of CAs with elevated mEVs and the total of the mEVs.

Six of these metrics are summarized in Fig. S23. It is noted that across all the domains examined, Δ8THC consistently outperforms other substance exposures and appears at the extreme right-hand end of each of the six graphs.

These results from categorical bivariate analysis therefore confirm the results obtained in continuous bivariate analysis, indicating that Δ8THC and Δ9THC are implicated with numerically more CAs and with a much higher effect size than traditional teratogens such as tobacco and alcohol.

Multivariable Panel Regression

Panel regression was the multivariable technique chosen to adjust the above bivariate results for common socio-demographic covariates including other substances, ethnicity and median household income. The substances studied were cigarettes, AUD, Δ8THC, analgesics and cocaine. All models were inverse probability weighted to allow formal causal inferences to be drawn. E-values were calculated in each case.

Additive Panel Model

Table S25 shows the 164 terms with positive and significant coefficients to emerge from additive panel modelling incorporating all socio-demographic substance and income covariates. Table S26 extracts 39 terms for Δ8THC. Urological, cardiovascular, chromosomal and gastrointestinal anomalies feature prominently.

Table S27 summarizes these findings by covariate listed in descending order of total mEV exponents. In this table, Δ8THC is noted to have higher results across all domains including the number of CAs implicated [39], mean mEV (2.789 × 10124), median mEV (2.09 × 1027), total negative P-value exponents (244) and total mEV exponents [1, 119].

These results are shown graphically in Fig. S24, where again Δ8THC is noted to occupy the extreme right-hand position on all four graphs for the number of implicated CAs, log of mean mEVs, log of total mEVs and sum of total negative P-value exponents.

Table 1 shows these results analysed by organ systems and listed in descending order of the total percentage of CAs in each organ group. This list is headed by chromosomal and body wall anomalies, each with 100% CAs affected. Cardiovascular, genitourinary, limb defects (66.7%) and gastrointestinal anomalies all have in excess of 50% of their CAs listed. central nervous system (CNS) and facial anomalies have 50% and 43% listed, respectively.

Table 1:

Summary table by system—additive multivariable panel model

System	Number of Anomalies Implicated	Total Number of Anomalies in System	% of Anomalies in System	Mean Lower Bound E-Value	Median Lower Bound E-Value	Total of P-Value Negative Exponents	Total of E-Value Exponents
Chromosomal	5	5	100.0	2.85E + 46	2.85E + 46	55	233
Body wall	3	3	100.0	8.78E + 28	8.78E + 28	18	82
CVS	15	19	78.9	1.90E + 21	1.90E + 21	78	404
GUT	4	6	66.7	1.68E + 27	1.68E + 27	34	86
Limb	2	3	66.7	1.50E + 28	1.50E + 28	7	45
GIT	4	7	57.1	7.08E + 40	7.08E + 40	24	131
CNS	3	6	50.0	7.00E + 23	7.00E + 23	11	59
Face	3	7	42.9	4.46E + 32	4.46E + 32	17	79

These results are shown graphically in Fig. 5, where again chromosomal, cardiovascular, limb, genitourinary and gastrointestinal anomalies feature prominently.

Figure 5:

Summary of numbers of affected anomalies and P- and E-values by organ system from additive multivariable panel model of Table S25: (a) number of CAs affected, (b) percentage of affected anomalies in each organ system, (c) total of mEV exponents by the organ system and (d) total of negative P-value exponents by the organ system

Interactive Panel Model

A similar model including all covariates was considered incorporating a three-way interaction between cigarettes: AUD and Δ8THC. It was also inverse probability weighted. Table S28 shows the 126 terms with positive and significant coefficients to emerge from this model. From this model, 13 terms including Δ8THC can be extracted (Table S29). Five of these mention the same CA twice. Extracting the most significant of these terms leaves the results shown in Table 2, which lists eight CAs. The Table also lists various metrics for the overall or “marginal” effect of covariates at the average of model covariates which is very useful for interactive models where different interactions may often result in opposing signs of model terms. This parameter is known as the average marginal effect (AME). In this case, it is measured in units of the standard deviation of the CAR for each CA. Facial (including holoprosencephaly), chromosomal, body wall, cardiovascular and gastrointestinal syndromes feature in the list of CAs shown in Table 2.

Table 2:

Summary table for significant anomalies with positive coefficients from interactive multivariable panel model P- and E-values, AME and percent AMEs from Table S28

Defect	Term	Mean CAR	SD CAR	Estimate	Std.Error	Model SD	T-Statistic	Adj.R2	P-Value	Covariate	AME	P-Value of AME	AME Lower Bound	AME Upper Bound	E-Value Estimate	E-Value Lower Bound
Anotia/microtia	Cigarettes: AUD: Δ8THC	2.23	0.76	25 618.243	9089.476	0.511	2.818	0.275	0.005	Δ8THC	0.618	3.25E–13	0.452	0.785	Infinity	Infinity
Atrioventricular septal defect	Cigarettes: AUD: Δ8THC	2.34	0.58	21 348.588	7648.960	0.429	2.791	0.034	0.006	Δ8THC	0.407	1.08E–08	0.267	0.546	Infinity	Infinity
Craniosynostosis	Δ8THC	2.05	0.89	4227.152	1813.896	0.541	2.330	0.239	0.022	Δ8THC	0.329	0.112328064	−0.077	0.734	Infinity	Infinity
Diaphragmatic hernia	Cigarettes: AUD: Δ8THC	1.66	0.37	15 896.251	5238.978	0.294	3.034	0.066	0.003	Δ8THC	0.281	9.27E–09	0.185	0.377	Infinity	Infinity
Holoprosencephaly	Cigarettes: AUD: Δ8THC	2.76	0.90	41 386.379	15 218.314	0.612	2.720	0.201	0.007	Δ8THC	0.317	0.034346152	0.023	0.611	Infinity	Infinity
Single ventricle	Δ8THC	0.98	0.53	2306.041	1152.334	0.384	2.001	0.049	0.047	Δ8THC	0.504	7.85E–08	0.320	0.688	Infinity	6.34E + 53
Small intestinal atresia/stenosis	Cigarettes: AUD: Δ8THC	1.85	0.45	16 720.952	7689.631	0.315	2.174	0.236	0.031	Δ8THC	0.370	1.29E–07	0.233	0.508	Infinity	Infinity
Trisomy 13	Δ8THC	2.45	0.69	2192.900	974.297	0.488	2.251	0.311	0.025	Δ8THC	0.381	2.76E–06	0.222	0.540	Infinity	Infinity

Note that the units for the AME are standard deviations of the applicable CA.

Table S28 may be summarized by covariate, as shown in Table 3. This table is also listed in the descending order of total mEV exponents. Here again, Δ8THC is noted to head the list of mean, median and total mEV exponents. These results are shown graphically in Figs 6 and 7. As judged by the marginal effect (Fig. 7), Δ8THC holds a mostly mid-rank position amongst the other covariates and demonstrates a lower degree of variability (highest mean to SEM of AME, panel d).

Table 3:

Summary table for key parameters from interactive panel model of Table S28

Covariate	Number of Anomalies	Total of Negative P-Value Exponents	Total of mEV exponents	Mean mEV exponents	Median mEV exponents	Mean AME	Median AME	SEM of AME	Total AME	Mean to SEM of AME Ratio	AME Lower Bound	AME Upper Bound
Δ8THC	9	10	2216	246.22	307	0.39	0.38	0.04	3.55	9.75	0.35	0.43
Cigarettes	15	27	1946	129.73	21	0.16	0.05	0.1	2.45	1.6	0.06	0.26
AUD	12	16	1391	115.92	32.5	−0.18	−0.07	0.07	−2.21	−2.57	−0.25	−0.11
NHPI	15	36	200	13.33	5	0.34	0.29	0.06	5.17	5.67	0.28	0.4
Cocaine	11	39	159	14.45	12	0.26	0.23	0.04	2.89	6.5	0.22	0.3
Analgesics	12	29	75	6.25	7.5	0.23	0.2	0.05	2.72	4.6	0.18	0.28
AIAN	17	55	32	1.88	2	0.24	0.22	0.03	4.1	8	0.21	0.27
Caucasian	14	35	7	0.5	0.5	0.52	0.43	0.07	7.25	7.43	0.45	0.59
African	10	24	0	0	0	0.36	0.37	0.06	3.6	6	0.3	0.42
Hispanic	17	106	0	0	0	0.47	0.39	0.07	7.97	6.71	0.4	0.54
Income	4	7	0	0	0	0.29	0.28	0.06	1.18	4.83	0.23	0.35

Figure 7:

Summary of AMEs by covariate from the interactive panel model summarized in Table S28. (a) Mean AME, (b) median AME, (c) total AME and (d) ratio of the mean to the standard error of the mean AME as a measure of AME variability by covariate

Summary of numbers of affected anomalies and P- and E-values by covariate from interactive multivariable panel model summarized in Table S28: (a) number of significant positive anomalies by covariate, (b) median mEV by covariate, (c) log of total mEV exponents and (d) log (base 10) of sum of negative P-value exponents

When body systems are considered chromosomal, body wall and CNS systems predominate (Table S30 and Fig. S25).

Interactive Panel Model at Two Temporal Lags

The above interactive model was again considered with all independent covariates lagged by 2 years. In this instance, 90 terms with positive and significant coefficients emerged, as shown in Table S31. It should be noted here that in R, the largest integer which is routinely considered is 1.7 × 10307 so that numbers larger than that are usually denoted as infinity. To make this analysable, this infinite number has been replaced by its minimum exponent (307) in all summary indices. Δ8THC was mentioned in eight of the terms mentioned in Table S31 (shown in Table S32) and for each of the four anomalies identified duplicate entries were recorded. When the most significant of these was retained (which was Δ8THC itself in each case), the results shown in Table S34 were obtained. Interestingly, the percentage AME in each case ranged from 2.1% to 10.5% SDs for each CA.

The results from Table S31 are summarized in Table S34, which is again listed in the descending order of total mEV exponents. Δ8THC is noted to have the highest position of median and total mEVs and occupies a mid-position for the number of CAs implicated and the sum of the P-value exponents (Fig. S26). Δ8THC has a lower position for summary indices of AMEs and a smaller range of variability (Fig. S27).

Table S31 may also be summarized by body organ systems for Δ8THC. These results are shown in tabular form in Table S35 and graphically in Fig. S28. Once again limb, genitourinary, cardiovascular and gastrointestinal systems feature in these graphs of percentage of anomalies affects, mEV exponents, P-value exponents and total AME.

Interactive Panel Model at Four Temporal Lags

This interactive model was again investigated with all independent covariates lagged by 4 years. In this instance, 78 terms with positive and significant regression coefficients were obtained which are presented in Table S36. In 18 cases terms included Δ8THC (Table S37) and when the most significant were selected out, nine terms remained (Table S38). Here, the percentage AME ranged up to 21.7% for biliary atresia.

Table S36 may be summarized by covariate as shown in Table S39, which is again listed in descending order of total mEV exponent. Once again, Δ8THC heads this list across several columns. As shown in Fig. S29, Δ8THC is at the extreme right-hand position for the number of anomalies implicated and the median and total mEV exponents. At four lags its marginal effect wanes (Fig. S30).

When the terms including Δ8THC from Table S36 are summarized by the system, the highest number of affected anomalies occurs in the limb, genitourinary and face systems, the highest mEVs are in the cardiovascular and genitourinary systems and the most affected systems are the gastrointestinal and cardiovascular systems (Table S40 and Fig. S31).

Discussion

Main Results

The data show that Δ8THC is associated with more CAs than tobacco and alcohol combined on bivariate analysis both before and after correction for estimates of ETOPFAs and also on categorical analysis. These changes were increased by adjustment in an additive multivariable model. Investigation in interactive models lagged to 0, 2 and 4 years reduced the numbers of CAs implicated but greatly increased the mEVs often to infinity (Table 2; Tables S33 and S38). Importantly since all multivariable models were inverse probability weighted, this analytical context allows causal inferences to be drawn from these results. Together, the use of inverse probability weighting and extremely elevated mEVs indicates that groups are broadly comparable and that uncontrolled confounding is effectively excluded on quantitative criteria.

Amongst continuous bivariate linear models corrected for estimates of ETOPFAs, the number of CAs was Δ9THC > Δ8THC > analgesics > tobacco > cocaine > AUD with 36, 33, 13, 12, 6 and 3 of 57 CAs affected (Tables S18–S23). In additive fully adjusted inverse probability-weighted panel models, the number of CAs implicated was Δ8THC > cocaine = cigarettes > analgesics > AUD > with 39, 10, 10, 12 and 3 anomalies affected (Table S27).

Interpretation

Our interpretation of these obviously serious results is both concerned and guarded. We are concerned that the effects demonstrated are significant and that a number of CAs along with most organ systems are strongly implicated. We are concerned that the effects we are seeing are reminiscent of genotoxic effects recently described for cannabis and other cannabinoids and that are consistent with earlier cellular studies [20, 72, 73], the effect appears to be a class effect implicating many cannabinoids [26, 27, 74]. We are concerned that much experimental and epidemiological evidence points towards a clear threshold effect beyond which exponential effects become commonplace [28-32]. And we are concerned that Δ8THC marketing is being aggressively pursued despite its class membership in a group of well-established genotoxic compounds. We are also very concerned that the quantitative criteria of causality are clearly fulfilled by the present findings.

Having said that the concentrations of Δ8THC in the samples analysed from Federal seizures during the period of this study were very low at about 0.1% of the population (Fig. S1). It thus seems mechanistically implausible to us that so low a population exposure could be driving birth anomaly patterns generally across the whole population. As Δ8THC generally occurs with other cannabinoids, it may be a marker for cannabinoid genotoxicity more broadly defined. Or perhaps the manufacturing processes that are used to generate it (strong acid, heating, etc.) also transform other cannabinoids, and Δ8THC is a marker for these products. More research will be required on these points. In this respect, therefore, our interpretation must remain guarded.

However, it is also clear that secular cannabinoid trends have moved a long way from 2014 and the penetration of this genotoxin into the community may now be much more significant, possibly even to the point where population effects might become feasible. For these reasons, it seems to us that further research into these issues is urgently required both in epidemiology and in laboratory studies of Δ8THC- and cannabinoid-related genotoxicity generally.

It is of interest that the findings in the present study contrast with those in many others which have not demonstrated a link between cannabinoid exposure and significant CAs [75-79]. Most such studies used much smaller samples and generally different methodologies including inpatient samples and self-report questionnaires. Some studies, for example, found a lack of association with gastroschisis [80], but it is clear from CDC–NBDPN studies [81] and many other reports including some from the same authors [8] and careful multivariate regression studies [82] that gastroschisis is indeed linked to cannabis use [8, 82–87]. Similarly, some studies have not demonstrated a link between Ebstein's anomaly and cannabis use [88], but reviews from the American Academy of Pediatrics and the American Heart Association have noted that it is indeed associated [89]. Therefore, a certain amount of disagreement in the literature seems not unusual.

It may be that the use of large national samples in the present study provides a more sensitive study methodology for these rare anomalies, particularly in an epidemiological context where cannabis use was formerly at a relatively low level. We also note the close concordance of the present results with recently published results from Europe, which similarly identified over 80 CAs as being linked with various metrics of cannabis exposure [90]. As noted in the Introduction section, it is not unusual in large regional or national studies to find many CAs linked with cannabis exposure [8–13, 19, 53, 91].

Causal Inference

It is of interest to consider how the present findings address the nine hallmarks of causality set out by Hill [92]. The findings reported herein demonstrate the strength of association, specificity to cannabinoids and not other substances, share appropriate temporality, are coherent with known data, are biologically plausible, demonstrate an appropriate dose–response curve and are consistent with experimental findings. As this is the first such report, it has not been replicated elsewhere and so cannot demonstrate consistency amongst studies.

One of the major pitfalls of observational studies is that the study groups may not be comparable. This aberration is corrected by the use of inverse probability weighting, which is the adjustment of choice for causally relevant modelling.

Another major hypothetical concern is that some unmeasured extraneous confounder covariate might exist, which explains away an apparently causal effect. However, such covariates cannot simply be fantasized. The E-value defines the quantitative correlation required of this unknown confounder covariate with both the exposure of concern and the outcome of interest in order to obviate the reported association. The usual mEV required for reporting causal relationships is 1.25 [93]. The E-value for the tobacco–lung cancer relationship is 9, which is regarded as high [94]. Clearly, our study where inverse probability-adjusted multivariable mEVs were infinite reports unusually high and in fact prohibitive E-values and therefore excludes unmeasured confounding on quantitative criteria.

Mechanisms

Consideration of the basic cellular mechanisms which might underpin the observed epidemiological findings is central to any consideration of a possible causal effect. This subject has been reviewed elsewhere [33, 39, 95]. Δ8THC is known to interfere with the synthesis of DNA, RNA and proteins [96-100] and thus necessarily interferes with nuclear and DNA physiology and cell growth like other cannabinoids [28, 101, 102].

Importantly, Δ8THC binds the type 1 cannabinoid receptor [2, 99] and inhibits mitochondrial activity [103-107]. Mitochondria are in close communication with the cell nucleus to which they supply ATP and numerous small molecular substrates for epigenetic reactions and DNA maintenance. Mitochondria also maintain close indirect communication with the nucleus through several small-molecule shuttles [108]. Hence, the disruption of mitochondria directly and indirectly perturbs and destabilizes genomic maintenance and normal epigenomic regulation.

The epigenome was recently shown to regulate a stunning array of fundamental cell functions including brain development [109], neuronal specification [109], neural network formation [110], learning [111], neocortical expansion [112, 113], glutamatergic neuronal specification [112] and long-term immune and metabolic memory held on lymphoid stem cells’ epigenome [114, 115].

It is also of relevance to observe that many cannabinoids interfere with key body morphogens including sonic hedgehog [32], retinoic acid [116-118], bone morphogenetic proteins [119-121], fibroblast growth factor [122] and many others.

For example, it is easy to see how cannabinoid exposure during a critical window of gestation at Days 24–26 could disrupt limb bud formation or digital development or, at different times, formation of the facial forebrain or central cardiovascular structures since all are dependent on sonic hedgehog signalling [123].

It has recently been formally demonstrated in Europe that the continent is undergoing a triple convergence of rising prevalence of cannabis use, increased intensity of cannabis daily use and increased cannabis potency [124, 125]. The same is true in USA and in many other places [15, 16, 19, 51, 52, 126].

Pathophysiologically, another very serious and extremely concerning convergence is also occurring as the exponential dose–response effects of cannabinoids [28-32], their entry into the food chain without consumer consent, damage to the genetic material at chromosomal, genetic and epigenomic–mitochondriopathic levels, their effects on CAs, cancerogenesis and population genomic ageing and multigenerational impacts increasingly come to bear.

The documentation in the present report that all the chromosomal anomalies affecting chromosomes 13, 18, 21, 22 and X are elevated in association with Δ8THC exposure is of particular concern, particularly in the light of previous findings implicating cannabinoids with elevated rates of acute lymphoid leukaemia of childhood [15], which entails major damage to chromosome 9 and testicular cancer [52] and which further entails damage to chromosome 12. Together the length of these chromosomes comprises 585 megabases of the human genome or ∼19.5%. Clearly, these findings directly implicate diverse cannabinoids in major human genotoxicity at the hundred megabase scale.

Also concerningly it appears that as the major CA and cancer databases lag behind current developments by several years, our public health systems are presently ill-prepared to monitor and track pertinent developments in real time.

Generalizability

Since this report is based on some of the best available data on this subject in the world, given the very strong effect sizes reported, given its detailed and comprehensive inclusion of the techniques of causal inference (inverse probability weighting and E-values) at many points and the consequent demonstration of putatively casual relationships, it seems to us that these results are likely to be widely generalizable wherever data of sufficient quality exist to investigate these relationships. However, given that this is an initial report, we are keen to see these results replicated elsewhere geographically and also in a more up-to-date temporal context.

Strengths and Limitations

This study has a number of strengths relating to data inputs including the use of a large national registry population-based database of CAs, the use of a large nationally representative sample of drug use and exposure and the use of cannabinoid exposure data from a major federally contracted analytical laboratory. It has several strengths relating to data analysis including the use of simultaneous multiple model investigation, the use of the quantitative tools of causal inference, particularly inverse probability weighting and E-values, and the use of multi-panelled graphs and maps for broad-ranging data visualization. Study weaknesses include the unavailability of solid data for ETOPFA rates and the unavailability of state-level cannabinoid concentration data. In common with many epidemiological studies, individual participant exposure data were not available to the present investigators. Also, we were not able to account for hereditary factors which might apply to anomalies such as orofacial clefts and neural tube defects, nutritional factors, paternal factors (contributing to factors such as non-disjunction) or the use of teratogenic medications taken during pregnancy (such as statins, warfarin, valproate, radiation, anticonvulsants and antidepressants).

Conclusion

This report documents that rising Δ8THC levels are closely associated with many CAs, particularly in the cardiovascular, central nervous, limb, gastrointestinal, genitourinary and face systems with high effect sizes and mEVs. Moreover, the results indicate that causal relationships are possible. Given the low concentration of Δ8THC in community samples, the present results do not exclude the possibility that Δ8THC is a marker for cannabinoids generally or another more toxic cannabinoid or combination thereof. The results are consistent with genotoxic effects reported for several other cannabinoids, and together these data indicate that genotoxicity is likely a class effect extending to many compounds in the cannabinoid class. Whilst the key quantitative metrics of formal causal inference are amply fulfilled by these results, we are keen to see these studies replicated elsewhere and in a more contemporaneous temporal context before formal causal attribution is assigned. Clearly, more laboratory research is required on Δ8THC genotoxicity, epigenotoxicity and embryonal and developmental toxicity. Meanwhile, we feel that these results clearly raise a high index of suspicion implicating Δ8THC along with many other cannabinoids as significantly genotoxic compounds. Of particular concern is that many cannabinoids demonstrate prominent threshold effects in the micromolar range beyond which genotoxicity, mitochondriopathy and epigenotoxicity exponentiate [28-32]. The presence of this threshold effect in the context of rapidly rising community penetration of a genotoxic agent is a grave concern indeed. We therefore feel that this initial report is a powerful stimulus to further research on Δ8THC in the laboratory, in the clinic and amongst epidemiological investigators. Whilst further work is awaited, it would seem prudent and precautionary to tightly control community penetration and marketing of Δ8THC in line with restrictive policies universally imposed on other compounds with potently genotoxic effects. The very serious intergenerational implications of the present findings for the health of subsequent generations across multiple major body systems are of great concern indeed.

Click here for additional data file.