Cannabinoid and substance relationships of European congenital anomaly patterns: a space-time panel regression and causal inferential study.

With reports from Australia, Canada, USA, Hawaii and Colorado documenting a link between cannabis and congenital anomalies (CAs), this relationship was investigated in Europe. Data on 90 CAs were accessed from Eurocat. Tobacco and alcohol consumption and median household income data were from the World Bank. Amphetamine, cocaine and last month and daily use of cannabis from the European Monitoring Centre for Drugs and Drug Addiction. Cannabis herb and resin Δ9-tetrahydrocannabinol concentrations were from published reports. Data were processed in R. Twelve thousand three hundred sixty CA rates were sourced across 16 nations of Europe. Nations with a higher or increasing rate of daily cannabis use had a 71.77% higher median CA rates than others [median ± interquartile range 2.13 (0.59, 6.30) v. 1.24 (0.15, 5.14)/10 000 live births (P = 4.74 × 10-17; minimum E-value (mEV) = 1.52]. Eighty-nine out of 90 CAs in bivariate association and 74/90 CAs in additive panel inverse probability weighted space-time regression were cannabis related. In inverse probability weighted interactive panel models lagged to zero, two, four and six years, 76, 31, 50 and 29 CAs had elevated mEVs (< 2.46 × 1039) for cannabis metrics. Cardiovascular, central nervous, gastrointestinal, genital, uronephrology, limb, face and chromosomalgenetic systems along with the multisystem VACTERL syndrome were particularly vulnerable targets. Data reveal that cannabis is related to many CAs and fulfil epidemiological criteria of causality. The triple convergence of rising cannabis use prevalence, intensity of daily use and Δ9-tetrahydrocannabinol concentration in herb and resin is powerfully implicated as a primary driver of European teratogenicity, confirming results from elsewhere.

Introduction

Whilst it is often said that prenatal cannabis exposure has relatively benign implications in postnatal life [1-3], recent independent reports from Hawaii [4], Colorado [5], Canada [6, 7], Australia [8] and USA [9-11] indicate that dozens of congenital anomalies (CAs; birth defects) are likely epidemiologically causally associated with rising rates of community cannabis consumption. Systems that are particularly affected include the cardiovascular, gastrointestinal, chromosomal, genitourinary, limb and body wall systems. Concerns regarding prenatal exposure were provided with heightened salience by reports from many places indicating increased use of cannabis and cannabinoid products by pregnant women in recent times [12], by rates of cannabis use in pregnancy amongst teenagers as high as 25% [13], by increased use of cannabis in pregnancy since the COVID-19 pandemic [12] and by reports that 69% of cannabis dispensaries positively recommend cannabis use to women whilst pregnant [14]. Moreover, recent reports note a quadruple convergence of rising rates of cannabis use, Δ9-tetrahydrocannabinol (THC) potency, intensity of daily use and cannabis use disorder in Europe, suggesting that the modern era is actually experiencing a confluence of concerning teratogenic trends [15, 16].

The implications of cannabinoid genotoxicity are further highlighted with the recent data suggesting that multiple cancers (of breast, pancreas, thyroid, liver and acute lymphoid and myeloid leukaemias) are also epidemiologically causally related to cannabis use [17, 18] and, with the formal experimental demonstration in mice, that epigenomic programming actually controls the organism-wide ageing epigenomic cassettes [19]. The recent demonstration that cannabis is a major driver of the rise in USA paediatric cancer rates underscores the transgenerational nature of this mutagenesis [17, 20]. These data together indicate that cannabinoid-related epigenomic disturbances likely have broad public health implications for diverse communities extending to cancerogenesis on the one hand and pan-systemic ageing on the other and including transgenerational effects.

Key to any consideration of the possible causal relationships of cannabis with mutagenesis and teratogenesis is the elucidation of the biological pathways, which may underlie any apparently causal relationship. Multiple cannabinoids have long been known to be toxic to chromosomes, genes, DNA strands, DNA nucleosides, the epigenome, sperm, oocytes, mitosis and meiosis and the male and female reproductive tracts in multiple respects [21-38]. Several studies have also shown cannabis to have a major effect perturbing DNA methylation [31-37], with these changes shown to be inheritable to subsequent generations [31-37], to perturb DNA methylation in the nucleus accumbens of offspring and affect behaviour [33, 34], to be seen in human sperm [31, 32] and to improve in both rats and humans after cessation of cannabis exposure [32]. Cannabis has an adverse effect on protein synthesis including histone formation, a change which necessarily opens up chromatin for aberrant gene expression [40, 41]. Thus cannabinoids derange the ‘histone code’. They also adversely affect tubulin synthesis and the post-translational modifications (particularly glycosylation) of tubulin, which have been collectively referred to as the ‘tubulin code’ [42], which adversely affects both the microtubules of the mitotic spindle and the anaphase separation of chromosomes, and also the motility of sperm flagella and their ability to maintain linear forward progression in fertilization assays [42]. Numerous cannabinoids have adverse effects on mitochondrial metabolism in neurons, sperm, lymphocytes and pulmonocytes [22, 23, 43–51], which necessarily impairs the supply of methyl, acetyl, ubiquinyl, propyl, adenosine-ribose and many other groups for the epigenomic machinery, impairs ATP energy supply for the numerous energy-dependent genomic and epigenomic reactions required for normal genome maintenance and also deranges the delicate mitonuclear balance [52, 53].

Ready acces to various European metrics of cannabis exposure indicating the quadruple convergence of rising cannabinoid exposure [15, 16] along with access to comprehensive congenital anomaly rates from multiple national European registries together with newer statistical techniques allowing the examination of multiple models in a single analytical run for all anomalies considered together in a space-time context have provided an ideal opportunity to investigate these relationships in the contemporary European context. The hypothesis to be tested in this investigation was that the well-described genotoxicity and epigenotoxicity of cannabinoids seen in vitro may be manifested clinically in vivo at the level of child population health with various of the described metrics for cannabis use. Furthermore, we sought to employ statistical techniques of formal causal inference to allow epidemiologically causal relationships to be investigated beyond merely those of simple association. It was also relevant to compare links described in other jurisdictions with the European findings and to compare the relative effects of the known teratogens tobacco and alcohol.

In terms of anomaly classes of special interest, we were particularly interested to study those that had been previously identified in the literature as being cannabis related, such as chromosomal and genetic, cardiovascular, central nervous, gastrointestinal, urogenital and nephrological and limb anomalies [4–9, 15–18, 20, 54, 55]. Interestingly, the presence in the European data of a rare multisystem anomaly known as vertebral, anorectal, cardiac, tracheo-esophageal fistula ± oesophageal atresia, renal anomalies and limb abnormalities (VATER/VACTERL), which was described from Great Ormond St Hospital as a group of co-occurring anomalies [56] whose aetiology was recently ascribed to inhibition of sonic hedgehog signalling in utero, which is a known target of many cannabinoids [57], was of particular interest. All hypotheses were formulated prior to study commencement.

Methods

Data

Data on total congenital anomaly rates per 10 000 live births was downloaded from the Eurocat website for each nation and for each year separately for all available CAs [58]. Data on national birth rates and populations were taken from the World Bank [59]. Data on tobacco (percentage smoking) and alcohol (per capita annual litres of alcohol consumption) use were taken from the World Health Organization (WHO) Global Health Observatory [60]. Data on drug use were taken from the European Monitoring Centre for Drugs and Drug Addiction (EMCDDA) [61]. EMCDDA data on cannabis use and potency were supplemented by data provided in the recent report from the Manthey group relating to monthly and daily cannabis use, which was itself derived from EMCDDA [16]. Median household income data were derived from the World Bank [62]. Nations were chosen on the basis of base population, the availability of comprehensive data for 2010–2019, and for their place in the Supplementary Fig. 4 of Manthey and colleagues relating to rising/high or falling/low daily cannabis exposure [16]. Nations were divided dichotomously into rising or falling groups based on the findings presented in this figure.

Derived Data

As noted in the Introduction, Europe has been subject to a convergence of rising indices of cannabinoid exposure in the past decade. It was therefore of interest to see if a combination of these variables may have more explanatory power than the simple covariates themselves. Hence, last month cannabis use was multiplied by the THC potency of cannabis herb and cannabis resin to form compound covariates. Similarly these metrics were also multiplied by daily use rates to gain compound indices of use-intensity-potency for each nation. Quintiles of substance exposure were calculated by dividing the total range across the whole period into five equal parts with the ggplot2 function cut_number.

Data Interpolation

Data interpolation was undertaken for drug use for years for which data were missing. For nations with no data relating to drug exposure in any year, their data field was allowed to remain entirely missing. Both raw data and data after interpolation are provided in the online files in the Mendeley data repository.

Statistics

Data were processed in R studio 1.4.1717 based on R 4.1.1 from the comprehensive R Archive Network. Data were log transformed based on the Shapiro–Wilk test. Negative or zero rates were arcsinh transformed where required. Normally distributed data are quoted as mean ± SEM. Non-normally distributed data are quoted as median and interquartile range. Data were manipulated in dplyr, and graphs were drawn in ggplot2, both from the tidyverse [63]. All artwork is original and is not under pre-existing copyright. Linear regression was conducted in Base R. Linear models were reduced by the classical method of manual deletion of the least significant term [64]. Overall or marginal effects of additive or interactive models were calculated from the margins package [65]. Point estimates for the E-value and its 95% lower bound were calculated using the R package E-value [66-69]. Relative risk (RR), attributable fraction in the exposed (AFE) and population attributable risk (PAR) were calculated in the package epiR version 2.0.38 [70]. The R package ranger was used to conduct random forest regressions [71], and the package vip was used to construct variable importance plots [72]. Heatmaps were drawn using the R Package gplots [73]. Multiple models were analysed simultaneously as described below using purrr-map pipelines (from tidyverse [63]) incorporating functions from the R packages broom [74, 75], dplyr, margins and E-value. P-values were corrected for multiple testing by the algorithms of Bonferroni [76], Holm [77] and false discovery rate [78].

Multivariable regression was conducted using panel regression using the R package plm [79]. Panel regression was chosen for several reasons, including the fact that space and time could be considered simultaneously without consuming degrees of freedom, that models could be inverse probability weighted (IPW), that temporal lagging could be conducted, that models can be incorporated into purrr-map analytical pipelines and that model objects contained a standard deviation that allowed E-values to be calculated. All panel models contained all drug and income covariates. All panel models were IPW. Panel models were conducted using the two-ways method, which allows space and time to be studied simultaneously.

Causal Inference

The use of inverse probability weighting makes all groups in observational studies comparable, effectively pseudo-randomizing the study groups and transforms findings from mere associations into a formal causal paradigm. For this, the R package ipw was employed [80]. One of the classic issues faced by observational studies is that low-level findings may be due to an unidentified confounder variable, which is not controlled in the analysis. This issue is addressed by the use of the E-value (or expected value) [66], which is the degree of association required of some unobserved extraneous covariate with both the exposure of concern and the outcome of interest to explain away the observed effect. Both a point estimate and the 95% lower confidence interval are calculated. In the published literature, minimum E-values (mEV) >1.25 are said to indicate potentially causal effects [81].

Data Availability

Input and output data have been provided online through the Mendeley data repository doi: 10.17632/vd6mt5r5jm.1. Four files with R source code running to 18 830 lines of code are also supplied.

Ethics

Ethical approval for this study was provided from the Human Research Ethics Committee of the University of Western Australia number RA/4/20/4724 on 24 September 2021.

Results

Introduction

As shown in eTable 1, data on 90 CAs were downloaded from the Eurocat website. The Eurocat dataset has a major advantage in that it provides a total rate of anomalies so that foetuses that are not live born, either due to the severity of their condition or due to early termination for anomaly, are included in the total rates. These anomalies are listed in abbreviated form in this table. The key to their full name is shown in eTable 2. Data were derived from the 16 nations indicated. Measures for the compound derived indices of cannabinoid exposure are also shown in eTable 1. eTable 2 also provides the system assignment used. In most cases, this is self-evident. However, the eye is derived from both face structures in its anterior segments, and the retina and optic nerve are derived from outgrowths of the forebrain [82]. Whilst lens and glaucoma abnormalities have been assigned alongside face structures, eye anomalies overall have been assigned to the central nervous system (CNS). eTable 2B lists the anomalies by organ system. A summary of the numbers of anomalies in each system is provided in eTable 3. eTable 4 provides a breakdown of the numbers of anomalies in each system by system both as numbers and as percentages of the totals.

Overall Picture

eTable 5 indicates the assignment of nations into those in which daily cannabis use was either high/increasing or low/decreasing as documented in eFig. 4 of Manthey [16]. When the log (as arcsinh) of the anomaly rate is graphed against time, the results depicted in Fig. 1 are shown. The median (interquartile range) in the decreasing group is 1.041 (0.149, 2.338) and in the increasing group is 1.50 (0.56, 2.54). These results correspond to raw (sinh) anomaly rates of 1.24 (0.149, 5.136) and 2.13 (0.589, 6.300)/10 000 live births, respectively, indicating a 71.77% elevation in the increasing cannabis use intensity countries (t = 8.204, df = 4660, P = 2.99 × 10−16). At linear regression against time, this finding is also significant (β-estimate = 0.2506 (0.1092, 0.3089), P = 4.74 × 10−17), which correspond to point estimates for the E-value of 1.63 and 1.52 for its lower bound (mEV). E-values exceeding 1.25 are said in the literature to indicate likely causal relationships [81].

Figure 1:

Overall trend of all CAs after log transformation

Substance and Time Trends—Continuous Analysis

eFigure 1A and B shows the rates of the 90 CAs across time. The figure has been split into two as there are so many anomalies to aid with presentation and readability. Genetic/chromosomal, cardiovascular and CNS anomalies are noted to feature amongst those which are rising. This list includes holoprosencephaly (a severe facial deformity) and VACTERL (a complex multisystem series of anomalies). eFigures 2A and 2B show the CA rate as a function of tobacco exposure. Only a few CAs are noted to rise. Again when alcohol is considered in eFig. 3, only a few anomalies are noted to rise. eFigures 4 and 5 perform similar roles for amphetamine and cocaine exposures, and it is noted that many more CAs are associated with positive gradients from these agents.

A similar exercise may be done for the THC concentration of cannabis herb. This is shown in eFig. 6. Here, all the CAs in eFig. 6A and half those in eFig. 6B are noted to demonstrate a positive relationship with rising herb THC concentration. This pattern is continued when last month cannabis exposure is considered in Fig. 2A and B, where all the CAs in Fig. 2A are noted to be rising with cannabis exposure. When daily cannabis use interpolated is considered, this pattern is again repeated as shown in Fig. 3A. The pattern is again repeated in the compound indices of last month herb and resin THC concentrations × interpolated daily use shown in eFigs 7 and 8.

Figure 2:

Time trends of CAs

Figure 3:

Trends of CAs with daily cannabis use (interpolated)

These relationships are formally analysed by linear regression in eTables 6–15 for each substance, respectively. These eTables list the usual metrics for linear models along with the applicable E-value point estimates and lower bounds. For the series of substances tobacco, alcohol, amphetamines and cocaine, 3, 12, 23 and 68 anomalies had elevated mEV, respectively. For the series of substances last month cannabis use, herbal THC concentration, resin THC concentration, daily use interpolated, last month cannabis use × herbal THC content × daily use interpolated (LMC_Herb_Daily) and last month cannabis use × resin THC content × daily use interpolated (LMC_Resin_Daily), the applicable numbers were 23, 45, 34, 41 and 42, respectively.

These data are summarized in Table 1. As the table is rather dense with information, these results are illustrated graphically for comparison in eFig. 9. Here, the number of anomalies for each substance is shown in panel A, and the cumulative exponents of the mEV in Panel B and the cumulative negative exponents of the P-values in Panel C. In each case, indices for cannabis exposure outperform teratogenic indices for tobacco, alcohol and amphetamines. eFigure 10 presents a study of the marginal or overall effects. Panel A shows the cumulative percentage average marginal effect, panel B the log of the mean percentage change, panel C the log of the standard error of the percentage change and panel D presents the SEM/average marginal effect ratio as a measure of the variability of the indices. In the first three cases, cannabis indices are noted to be higher than those of the other substances. The variability of the cannabinoid indices is lower than that of tobacco, alcohol and amphetamines (Panel D).

Table 1:

Summary table of bivariate continuous relationships by substance

Substance	Number of terms with elevated E-values	Sum of P-value exponents	Sum of mEV exponents	Total % change	Mean % change	S.E. % change	Median % change	First quartile % change	Third quartile % change
Tobacco	3	3	0	6.45	2.15	1.56	1.18	0.63	3.19
Alcohol	12	24	0	99.40	8.28	2.79	6.88	2.21	8.43
Amphetamine	23	45	0	185.25	8.05	2.51	3.21	0.72	10.20
Cocaine	68	268	0	991.55	14.58	1.65	10.71	3.19	20.85
LM_Cannabis	23	48	35	15 019.46	653.02	111.64	468.06	250.90	959.87
Herb_THC	45	107	55	9780.06	217.33	40.28	112.64	26.28	295.60
Resin_THC	34	75	8	2261.72	66.52	12.01	47.18	12.18	86.21
LM_Cann_×_Herb	42	107	68	20 370.06	485.00	73.45	325.67	146.58	653.78
LM_Cann_×_Resin	38	72	3	3500.09	92.11	15.32	61.13	15.19	132.36
Daily_Interpol	41	116	241	61 794.77	1507.19	226.63	1145.78	467.88	2180.68
Herb_×_Day_Int	41	111	18	6736.93	164.32	24.53	131.28	59.70	230.31
Resin_×_Day_Int	42	90	NA	2860.11	68.10	11.06	50.12	11.35	92.23

Since some marginal effect data are not distributed normally, the median and first and third quartile data are shown in eFig. 11. In all cases, the cannabis indices are substantially higher than those of the other substances. eTable 16 presents a categorization of the organ systems affected by their substance exposures by numbers of anomalies. These data are also presented as a heatmap in Fig. 4. Cardiovascular, chromosomal, gastrointestinal and uronephrological anomalies are noted to be prominent. eFigure 12 presents the number of systems impacted by each substance.

Figure 4:

Heatmap of systems by substance exposure

In the genomic and epigenomic literature, volcano plots are commonly used to represent the significance levels against the fold change in gene expression. The equivalent in this work might be to chart the significance level against the mEV as a measure of fold change implied in the data. eFigures 13, 14 and 15 do this for tobacco, alcohol and amphetamines. The eFigures are noted to be rather lean and to have relatively low levels of significance and mEV. eFigure 16 performs a similar role for cocaine, and this eFigure is noted to be heavily populated and quantitatively much greater. eFigure 17, Figs 5 and 6 and eFigs 18 and 19 perform this role for herbal THC content, last month cannabis use, daily cannabis use interpolated, LMC_Herb_Daily and LMC_Resin_Daily, respectively. These figures are noted to be more densely populated and quantitatively much higher than those for the other substances.

Volcano plot of significance (negative log (P-value)) against log (mEV) for past month cannabis use

Categorical Analysis

The data also lend itself to categorization for the purposes of calculating key epidemiological indices, including RR, AFE and PAR. For this purpose, substance exposure was divided into quintiles and boxplots as shown in eFigs 20–24 for tobacco, alcohol, amphetamines, cocaine and LMC_Resin_Daily, which contrast the highest and lowest quintiles of substance exposures. Where notches do not overlap, this indicates a statistically significant difference. This method also allows the calculation of highest:lowest quintile ratios for each anomaly by substance.

As shown in eTable 17, many substances demonstrate higher anomaly rates in the highest quintiles. However, the indices of cannabis exposure, which include daily exposure, have the greatest number of elevated ratios. These data are illustrated graphically in eFig. 25. eTable 18 shows the number of organ systems affected for LMC_Resin_Daily as a function of the total number of anomalies in each organ system. The table is ordered from those with the greatest percentage of anomalies per system. In fact, all 11 measured body systems are represented at levels above 50% in these data with general, chromosomal, uronephrolgical, limb, body wall, respiratory, cardiovascular, face, gastrointestinal and CNS anomalies more than 80% represented. These data are presented graphically in eFig. 26, which highlights these findings. eTables 19–29 and Tables 2 and 3 show the formal quantitative analysis of this data concatenated as a series of two-by-two epidemiological tables in long format. Many interesting features emerge from these tables, including that foetal alcohol syndrome is at the top of the tobacco, alcohol and herb THC, last month herb THC, last month Resin THC and LMC_Resin_Daily tables with mEVs of 8.67, 34.24, 13.44, 9.69, 7.37 and 43.19, respectively. Interestingly, VACTERL syndrome is near the top of the alcohol, amphetamine, cocaine, last month cannabis, last month cannabis × daily cannabis interpolated, last month cannabis herb, LMC_Herb_Daily and LMC_Resin_Daily lists, with mEVs of 5.93, 4.64, 24.90, 11.35, 26.43, 4.92, 34.61 and 43.19, respectively. t-test results are listed along with their accompanying P-values in eTable 29.

Volcano plot of significance (negative log (P-value)) against log (mEV) for daily cannabis use interpolated

The numbers of anomalies with elevated mEVs from this analysis by substance are shown in tabular and graphical formats in eTable 30 and eFig. 27.

If one studies the lists for the cannabis metrics closely, it is noted that 89 of the 90 (98.9%) CAs demonstrate elevated mEVs by one or more indices of cannabinoid exposure. The sole exception is indeterminate sex. This anomaly, however, is mentioned positively in the tables shown below. Complete coverage of this set of 89 CAs can be achieved by considering together the three covariates cannabis herb THC concentration, daily cannabis use and LMC_Resin_Daily.

Forest Regression

Since the analysis clearly needs to move from bivariate into multivariable regression, a salient issue is which variable/s should be used as the key metric of cannabis exposure moving forwards? This issue is not immediately apparent in the rich European dataset.

Additive and interactive linear model of all covariates against (log) CA rates were constructed, and the final forms of these linear models after model reduction are shown in eTable 31. A three-way interaction between tobacco, alcohol and LMC_Herb_Daily was used in the interactive model. Forest regression was conducted on these models using the Gini (‘impurity’) index as the measure of variable importance. This procedure derived the results shown in tabular form in eTable 32 and in graphical form as variable importance plots in Fig. 7. It is clear from these analyses that for both models compound indices of daily cannabis exposure and cannabis herb THC concentration are the most powerful covariates for entry into multivariable regression.

Figure 7:

Variable importance plots for (A) additive and (B) interactive random forest models (at zero lags)

Multivariable Regression

Additive Models

For these reasons, the three covariates cannabis herb THC concentration, LMC_Herb_Daily and LMC_Resin_Daily were selected for use in multivariable regression. They were first applied together with the other substances and median household income in an additive panel regression model. All panel regression models were IPW. eTable 33 shows the result of these regressions for 196 statistically significant positive terms From this eTable, 86 terms were extracted, which included a term for the metrics of cannabis (eTable 34). The most significant terms for each separate CA were then selected out, leaving 76 anomalies, which appear in Table 4. Interestingly, this table is headed by cardiovascular disorders, genital disorders, microphthalmos, polydactyly, nervous system disorders, patent ductus arteriosus (PDA), atrial septal defect (ASD) and genetic and facial anomalies, which have all been previously reported to be cannabis associated [4–8, 18, 55]. It is also noted that the P-values in this table ascend from 1.5 × 10−23 (for teratogenic syndromes) and the mEV descend from 8.2 × 109.

Table 4:

Multivariate additive IPW panel regression results

Anomaly	Term	Mean rate	β-estimate	Standard error	Sigma	T-statistic	Adj. R2	P-value	E-value estimate	E-value (95% lower bound)	P-value exponent	E-value exponent	% Increment	P—Bonferroni	P—false discovery rate	P—Holm
Congenital heart	Herb	82.6649	8.8588	0.9896	0.2846	8.9516	0.5218	5.16E-14	3.99E + 12	8.19E + 09	14	9	3.07	3.82E-12	9.55E-13	3.67E-12
Genital	Herb	22.6289	7.7667	1.3378	0.3848	5.8055	0.5099	9.98E-08	1.90E + 08	3.89E + 05	8	5	11.20	7.39E-06	8.21E-07	6.59E-06
An/microphthalmos	LpmHerbDailyInt	0.9004	12.2684	1.3776	0.7287	8.9056	−0.0865	4.98E-14	9.00E + 06	3.11E + 05	14	5	281.40	3.68E-12	9.55E-13	3.58E-12
Aortic valve S/A	LpmHerbDailyInt	1.3094	14.1512	1.6110	0.8522	8.7841	−0.0520	8.93E-14	7.31E + 06	2.53E + 05	14	5	193.51	6.61E-12	1.32E-12	6.25E-12
Polydactyly	Herb	10.3496	9.1270	1.6788	0.4957	5.4367	0.3313	4.52E-07	3.79E + 07	9.12E + 04	7	4	24.48	3.34E-05	2.57E-06	2.80E-05
Nervous system	Herb	26.1182	7.7253	1.4583	0.4194	5.2975	0.3384	8.58E-07	3.80E + 07	7.80E + 04	7	4	9.70	6.35E-05	4.53E-06	5.23E-05
PDA	LpmHerbDailyInt	2.9086	11.1711	1.5486	0.8192	7.2137	0.3160	1.57E-10	4.90E + 05	1.69E + 04	10	4	87.11	1.16E-08	1.94E-09	1.08E-08
ASD	Herb	21.4739	9.2217	1.9147	0.5653	4.8163	0.1194	5.80E-06	5.59E + 06	1.35E + 04	5	4	11.80	4.29E-04	2.42E-05	3.37E-04
Patau syndrome	LpmHerbDailyInt	1.7507	9.5652	1.3862	0.7333	6.9003	0.0513	6.74E-10	2.86E + 05	9.88E + 03	10	3	144.73	4.99E-08	7.13E-09	4.59E-08
Ear, face and neck	Herb	2.4050	18.1635	4.1690	1.2309	4.3568	0.0233	3.46E-05	1.36E + 06	3.27E + 03	4	3	105.35	0.00255941	1.07E-04	0.00179197
Teratogenic synds	LpmResinDailyInt	1.2227	5.9944	0.4423	0.6403	13.5526	0.4714	1.50E-23	1.00E + 04	2.93E + 03	23	3	178.34	1.11E-21	1.11E-21	1.11E-21
Arterial truncus	Herb	1.0030	10.0572	2.3589	0.6965	4.2636	0.0782	4.91E-05	1.02E + 06	2.45E + 03	4	3	252.62	0.00363094	1.40E-04	0.00240427
Hypospadias	Herb	19.9122	7.3602	1.7734	0.5101	4.1502	0.4986	7.65E-05	1.01E + 06	2.07E + 03	4	3	12.72	0.00566086	2.02E-04	0.00359541
Orofacial clefts	Herb	14.1998	8.6815	2.1057	0.6057	4.1229	0.1253	8.45E-05	9.25E + 05	1.90E + 03	4	3	17.84	0.00625453	2.09E-04	0.00388795
Posterior urethral valve	Herb	1.2952	8.8449	2.2296	0.6583	3.9670	0.2628	1.45E-04	4.08E + 05	982.80	3	2	195.63	0.01072607	3.15E-04	0.00594282
Matern infect malform	LpmResinDailyInt	0.8059	5.4136	0.4771	0.6906	11.3481	0.3104	4.07E-19	2.51E + 03	732.22	19	2	270.57	3.02E-17	1.51E-17	2.97E-17
Urinary	Herb	36.0821	6.0183	1.5679	0.4629	3.8386	0.1964	0.0002	2.75E + 05	661.40	3	2	7.02	0.016893	4.57E-04	0.00867478
Limb	Herb	37.7465	5.5628	1.5229	0.4380	3.6527	0.5350	0.0004	2.09E + 05	428.32	3	2	6.71	0.0325785	7.97E-04	0.01540875
Tetralogy of Fallot	LpmHerbDailyInt	3.1167	8.3069	1.7397	0.9203	4.7748	−0.0571	6.85E-06	7.38E + 03	254.72	5	2	81.30	5.07E-04	2.55E-05	3.83E-04
Hydrocephalus	LpmHerbDailyInt	5.3737	7.0000	1.4664	0.7757	4.7735	−0.0047	6.88E-06	7.37E + 03	254.19	5	2	47.15	5.09E-04	2.55E-05	3.83E-04
Multicystic renal dys	Herb	3.5369	9.7517	2.7887	0.8234	3.4969	0.3369	0.0007	9.59E + 04	230.40	3	2	71.64	0.05402309	0.00125635	0.02336134
Syndactyly	Herb	4.2302	7.3389	2.1361	0.6307	3.4357	0.0478	0.0009	7.94E + 04	190.72	3	2	59.90	0.0660369	0.00150084	0.02766411
Encephalocele	LpmHerbDailyInt	1.1446	6.2213	1.4122	0.7470	4.4055	−0.0748	2.88E-05	3.91E + 03	134.71	4	2	221.37	0.00212862	1.01E-04	0.00155332
Annular pancreas	LpmHerbDailyInt	0.2372	3.0116	0.6855	0.3626	4.3930	−0.0313	3.02E-05	3.83E + 03	131.85	4	2	1068.20	0.00223154	1.01E-04	0.00159826
Severe microcephaly	LpmHerbDailyInt	2.6832	7.5642	1.7358	0.9182	4.3578	0.1977	3.45E-05	3.60E + 03	124.06	4	2	94.43	0.00255012	1.07E-04	0.00179197
Spina bifida	Herb	4.3705	5.6604	1.7243	0.5091	3.2828	0.1039	0.0015	4.96E + 04	118.87	2	2	57.97	0.10797166	0.00234721	0.04231322
Down syndrome	LpmHerbDailyInt	20.7996	3.3573	0.7784	0.4117	4.3132	0.2441	4.08E-05	3.34E + 03	114.86	4	2	12.18	0.00301611	1.21E-04	0.00203791
VATER/VACTERL	Herb	0.4798	6.1281	1.9199	0.5522	3.1919	0.3359	0.0020	4.86E + 04	99.31	2	1	528.09	0.14505544	0.00302199	0.05315528
Situs inversus	LpmHerbDailyInt	0.6187	3.5795	0.8573	0.4535	4.1752	0.2135	6.81E-05	2.63E + 03	90.49	4	1	409.53	0.00503876	1.87E-04	0.00326839
Diaphragmatic hernia	LpmHerbDailyInt	2.5089	7.1983	1.7493	0.9253	4.1151	−0.0659	8.49E-05	2.37E + 03	81.54	4	1	100.99	0.00628208	2.09E-04	0.00388795
Cleft lip ± palate	LpmHerbDailyInt	8.6449	5.8082	1.4199	0.7511	4.0904	0.0181	9.29E-05	2.27E + 03	78.14	4	1	29.31	0.00687259	2.22E-04	0.00408641
Anotia	LpmHerbDailyInt	0.2330	1.8699	0.4817	0.2548	3.8816	0.0073	0.0002	1.59E + 03	54.40	3	1	1087.45	0.01452356	4.03E-04	0.00765431
Severe CHD	LpmHerbDailyInt	21.9866	4.5125	1.1816	0.6251	3.8189	−0.0172	0.0002	1.43E + 03	48.79	3	1	11.52	0.01809651	4.66E-04	0.00904826
Hydronephrosis	LpmHerbDailyInt	13.5286	5.9134	1.5490	0.8194	3.8176	0.2531	0.0002	1.42E + 03	48.68	3	1	18.73	0.01817657	4.66E-04	0.00904826
Vascular disruptions	LpmResinDailyInt	6.9668	2.9516	0.4384	0.6346	6.7330	0.1494	1.64E-09	137.29	39.78	9	1	31.30	1.21E-07	1.51E-08	1.10E-07
Genetic syndromes	LpmHerbDailyInt	6.5208	6.4410	1.7660	0.9342	3.6473	−0.0163	0.0004	1.06E + 03	36.18	3	1	38.86	0.03265865	7.97E-04	0.01540875
Hypoplastic left ht	LpmHerbDailyInt	2.3257	5.8086	1.6505	0.8731	3.5193	0.0102	0.0007	851.29	28.93	3	1	108.95	0.05015837	0.00119425	0.02236792
VSD	Herb	37.2628	5.3939	1.9550	0.5772	2.7591	0.1610	0.0070	9.87E + 03	23.25	2	1	6.80	0.5185063	0.00925904	0.13313
Gastroschisis	LpmResinDailyInt	2.3468	3.5118	0.6263	0.9067	5.6071	0.1511	2.19E-07	67.37	19.34	7	1	92.92	1.62E-05	1.52E-06	1.42E-05
Blader extr/epispad	LpmResinDailyInt	0.7056	1.9668	0.3513	0.5085	5.5991	−0.0206	2.26E-07	67.03	19.24	7	1	309.03	1.67E-05	1.52E-06	1.45E-05
Indeterminate sex	Herb	0.5005	5.2208	1.9503	0.5758	2.6769	−0.0617	0.0088	7.66E + 03	17.93	2	1	506.25	0.65226722	0.01105538	0.14179341
Cong. glaucoma	LpmResinDailyInt	0.2853	2.2735	0.4156	0.6016	5.4708	0.0968	3.91E-07	61.80	17.70	7	1	764.30	2.89E-05	2.41E-06	2.46E-05
Ebstein’s anomaly	Herb	0.4111	3.4408	1.3011	0.3842	2.6445	0.0493	0.0096	6.93E + 03	16.18	2	1	616.34	0.71306363	0.01188439	0.14453992
All anomalies	LpmHerbDailyInt	249.8954	2.1805	0.6920	0.3661	3.1510	0.2718	0.0022	451.56	15.11	2	1	1.01	0.16299042	0.00332634	0.0572669
Chromosomal	LpmHerbDailyInt	35.9025	2.9271	0.9356	0.4949	3.1286	0.1256	0.0024	434.41	14.51	2	1	7.06	0.17465658	0.00349313	0.0590056
Cong. cataract	LpmResinDailyInt	1.1682	3.0794	0.5990	0.8672	5.1406	0.0874	1.56E-06	50.12	14.29	5	1	186.66	1.15E-04	7.68E-06	9.35E-05
Abdominal wall defx	LpmResinDailyInt	5.7704	2.2561	0.4415	0.6392	5.1100	0.3983	1.77E-06	49.16	14.00	5	1	37.79	1.31E-04	8.17E-06	1.04E-04
Craniosynostosis	LpmHerbDailyInt	2.9964	4.9570	1.6039	0.8485	3.0905	0.2027	0.0027	406.88	13.56	2	1	84.56	0.19616883	0.00384645	0.06362232
AVSD	Herb	4.0854	6.2821	2.4410	0.7207	2.5735	−0.0096	0.0117	5.57E + 03	12.89	1	1	62.02	0.86452845	0.0141726	0.16355944
Lateral anomalies	Herb	1.8149	6.2780	2.4679	0.7098	2.5438	0.1470	0.0127	6.26E + 03	12.32	1	1	139.61	0.94057227	0.01517052	0.16523567
Limb reductions	LpmResinDailyInt	5.3892	2.1126	0.4390	0.6355	4.8126	0.1210	5.89E-06	40.69	11.52	5	1	40.46	4.36E-04	2.42E-05	3.37E-04
Mitral valve anomalies	LpmHerbDailyInt	1.6067	5.6814	1.9368	1.0205	2.9334	0.3100	0.0043	316.61	10.29	2	1	157.70	0.31625901	0.00596715	0.09614151
Conjoined twins	Herb	0.1128	1.9128	0.7687	0.2270	2.4882	−0.0424	0.0147	4.28E + 03	9.78	1	0	2246.24	1	0.01721843	0.17590719
Skeletal dysplasias	LpmHerbDailyInt	1.8518	4.7853	1.6649	0.8807	2.8743	−0.0040	0.0050	280.35	9.18	2	0	136.83	0.37296464	0.00678118	0.10080125
Transpos grt vess	Herb	3.4440	5.7348	2.3746	0.7011	2.4151	0.1216	0.0177	3.42E + 03	7.70	1	0	73.57	1	0.0201855	0.17730503
Anorectal S/A	LpmResinDailyInt	3.1853	1.9684	0.4825	0.6985	4.0795	0.1261	9.67E-05	25.47	7.06	4	0	68.46	0.00715222	2.24E-04	0.00415602
Choanal atresia	LpmResinDailyInt	0.9191	1.9449	0.4883	0.7070	3.9826	0.0914	1.37E-04	23.94	6.61	3	0	237.25	0.01014212	3.07E-04	0.00575634
Foetal alcohol	LpmResinDailyInt	0.2577	0.8498	0.2158	0.3124	3.9386	0.3406	0.0002	23.27	6.41	3	0	846.15	0.0118676	3.39E-04	0.00641492
Single ventricle	LpmResinDailyInt	0.7525	1.2747	0.3735	0.5407	3.4128	0.0451	0.0010	16.57	4.43	3	0	289.77	0.07114277	0.00158095	0.02884166
Hip dysplasia	LpmResinDailyInt	6.3491	2.3960	0.7482	1.0831	3.2023	0.2890	0.0019	14.45	3.80	2	0	34.34	0.1404818	0.00298897	0.05315528
Cleft palate	Herb	5.4928	6.5824	2.9987	0.8625	2.1951	−0.0213	0.0308	2.07E + 03	3.68	1	0	46.13	1	0.03302199	0.21425568
Eye	LpmResinDailyInt	3.6805	1.3682	0.4656	0.6741	2.9384	0.2060	0.0042	12.16	3.11	2	0	59.25	0.30932487	0.00594856	0.09614151
Respiratory	LpmResinDailyInt	3.2843	1.7470	0.5971	0.8644	2.9258	0.5774	0.0044	12.06	3.08	2	0	66.39	0.32335886	0.00598813	0.09614151
Bilat renal agenesis	LpmHerbDailyInt	1.5241	3.3242	1.4461	0.7649	2.2988	−0.0651	0.0238	103.84	3.01	1	0	166.25	1	0.02669179	0.21425568
PV atresia	LpmHerbDailyInt	1.1589	3.3990	1.4839	0.7850	2.2905	−0.0789	0.0243	102.37	2.95	1	0	218.63	1	0.02684011	0.21425568
Bile duct A	LpmResinDailyInt	0.2854	0.9697	0.3556	0.5147	2.7273	0.1397	0.0077	10.58	2.63	2	0	764.03	0.56702174	0.00994775	0.13792421
Edward syndrome	LpmResinDailyInt	4.9591	1.4804	0.5490	0.7947	2.6968	0.2853	0.0083	10.37	2.57	2	0	43.97	0.61721838	0.0106417	0.14179341
Neural tube defects	LpmHerbDailyInt	9.0634	2.2506	1.0304	0.5451	2.1841	0.1583	0.0315	85.16	2.32	1	0	27.96	1	0.03303907	0.21425568
Holoprosencephaly ∼	LpmResinDailyInt	1.6371	1.0551	0.4262	0.6170	2.4756	0.1276	0.0152	8.95	2.12	1	0	133.20	1	0.01751788	0.17590719
Hypoplastic right ht	LpmHerbDailyInt	0.6594	2.3160	1.0852	0.5741	2.1342	−0.0720	0.0355	78.11	2.06	1	0	384.25	1	0.03650427	0.21425568
Aortic atresia ∼	LpmResinDailyInt	0.5426	0.8590	0.3779	0.5470	2.2732	−0.0431	0.0254	7.81	1.74	1	0	401.89	1	0.02769426	0.21425568
Digestive system	Herb	17.0099	3.1511	1.5699	0.4515	2.0072	0.1302	0.0478	1.15E + 03	1.63	1	0	14.90	1	0.04779159	0.21425568
Club foot	LpmResinDailyInt	10.6471	1.0225	0.4687	0.6785	2.1818	0.3969	0.0317	7.35	1.57	1	0	20.48	1	0.03303907	0.21425568
Small intestine S/A	LpmResinDailyInt	1.0797	1.1091	0.5432	0.7863	2.0420	−0.0879	0.0440	6.68	1.30	1	0	201.96	1	0.04465153	0.21425568

These data are summarized in Table 5, which summarizes the key metrics for these results. mEV exponents are illustrated graphically in Fig. 8 and negative P-value exponents are shown graphically in Fig. 9. Figure 10 shows the number of significant associated anomalies (Panel A), the cumulative E-value exponents (Panel B) and the cumulative P-value negative exponents (Panel C). Figure 11 portrays the overall marginal effects as (A) sum, (B) mean and (C) median values.

Table 5:

Summary table of multivariate additive IPW panel regression results by substance

Covariate	Number of positive terms	Sum mEV exponents	Sum P-value exponents	Mean % increment	Median % increment	Sum % increment
Alcohol	16	0	35	−11.22	−6.37	−179.46
Herb_THC	28	57	92	184.08	41.25	5154.28
Amphetamines	13	0	22	8.19	2.75	106.52
Cocaine	28	0	117	−31.23	−14.19	−874.33
LM.Cann_Resin_THC_Daily.Int	28	41	116	187.10	97.71	5238.80
Median_Income	30	16	156	211.34	113.06	6340.35
LM.Cann_Daily.Int	37	0	77	−16.78	−8.99	−620.91
Tobacco	22	0	40	2.09	0.68	45.91

Figure 8:

Scatterplot of log (base 10) mEV by congenital anomaly type from additive multivariable IPW panel model

Figure 9:

Scatterplot of negative log (base 10) P-value by congenital anomaly type from additive multivariable IPW panel model

Figure 10:

Summaries of E- and P-values by substance type for (A) number of anomalies with elevated mEV, (B) the sum of the mEV exponents and (C) the sum of the negative exponents of the significant P-values for the additive IPW multivariable panel model

Figure 11:

Summaries of marginal (overall) effects by substance type for (A) total percentage change at average marginal effect (AME), (B) the mean percentage change at AME and (C) the median percentage change at AME for the additive IPW multivariable panel model

Table 6 shows these data by organ system. It is noted that the table is headed by CAs affecting the face, genitalia, limbs and uronephrological systems, each of which show 100% of anomalies affected. Summary metrics relating to E- and P-values are shown in Fig. 12, and studies relating to marginal effects are illustrated in Fig. 13.

Table 6:

Summary table of multivariate additive IPW panel regression results by organ system

System	Number positive terms	Total system count	% Anomalies in system	Sum mEV exponents	Sum P-value exponents	Sum % increment	Mean % increment	Median % increment
Face	9	9	100.00	10	32	2607.48	289.72	133.20
Genital	2	2	100.00	6	10	517.44	258.72	258.72
Limb	6	6	100.00	9	21	186.38	31.06	29.41
Uronephrology	7	7	100.00	11	24	781.02	111.57	71.64
Cardiovascular	19	23	82.61	36	85	3004.32	158.12	87.11
Central nervous	9	11	81.82	18	41	883.79	98.20	59.25
General	9	11	81.82	9	57	4756.37	528.49	270.57
Body wall	3	4	75.00	3	16	231.69	77.23	92.92
Gastrointestinal	5	8	62.50	2	12	2117.53	423.51	201.96
Chromosomal	4	7	57.14	5	17	234.61	58.65	41.41
Respiratory	1	2	50.00	0	2	66.39	66.39	66.39

Figure 12:

Summaries of E- and P-values by organ system. (A) Number of anomalies with elevated mEV. (B) Percentage of anomalies with elevated mEV. (C) The sum of the mEV exponents. (D) The sum of the negative exponents of the significant P-values for the additive IPW multivariable panel model

Figure 13:

Summaries of marginal (overall) effects by organ system for (A) total percentage change at average marginal effect (AME), (B) the mean percentage change at AME and (C) the median percentage change at AME for the additive IPW multivariable panel model

This exercise can be repeated for the three cannabis metrics providing complete coverage of the 89 anomalies shown in bivariate analysis to be cannabis related. As shown in eTables 35–41 and eFigs 28–33, this procedure selects out 69 CAs and finds that last month cannabis use × daily cannabis use interpolated is by far the most powerful predictive covariate in this set.

Interactive Panel Models

Zero Temporal Lags

An IPW interactive model, including a three-way interaction between tobacco, herb THC concentration and LMC_Resin_Daily, was examined. Two hundred sixty-seven significant positive terms were returned from this procedure (eTable 44), of which 155 included cannabis-related terms (eTable 45, ordered by anomaly). These are ordered by P-value (eTable 45), by mEV (eTable 46) and by marginal effect size (eTable 47). Finally, the most significant terms for the 76 CAs implicated are listed in the descending order of mEV from 2.87 × 1016 in eTable 48. Data are summarized by substance exposure term in eTable 49. P-value negative exponents (eFig. 34) and mEV exponents (eFig. 35) are aggregated and illustrated in eFig. 36, and summary data for marginal effects are shown in eFig. 37.

Tabular data summarizing these data by organ system are shown in eTable 50, and they are illustrated graphically in eFigs 38 and 39.

Two Temporal Lags

When all the independent variables are lagged to two years, the results shown in eTable 51 are obtained for significant positive terms in an IPW model featuring an interaction between tobacco and LMC_Resin_Daily. The results are summarized in eTable 52, which show that LMC_Resin_Daily is the most salient term by number of anomalies implicated and the sum of negative P-value exponents and mEV exponents and that herb THC concentration is the most salient factor by marginal effect size estimates. eTable 53 extracts cannabis-related terms noting duplicate entries for congenital heart disease under various cannabis-related terms and PDA. eTable 54 lists these anomalies in order, eTable 55 lists them by ascending P-values and eTable 56 lists them by descending mEV. In eTable 56, duplicate anomaly entries have been removed, leaving only the most signifcant terms. eTable 57 lists the applicable marginal effects. eFigure 40 lists these P-value negative exponents graphically, and eFig. 41 similarly lists the exponents of mEV. Summary data for numbers of anomalies and P- and E-values are shown in eFig. 42, and cannabis-related terms are noted to clearly outperform tobacco, alcohol and other substances. Similar findings are clearly shown upon consideration of marginal effects in eFig. 43.

Data are summarized by organ system in eTable 58 and presented graphically in eFigs 44 and 45. When considering E- and P-values, the cardiovascular system and CNS are the most affected (eFig. 44). When marginal effects are considered, the gastrointestinal tract is also particularly affected (eFig. 45).

Four Temporal Lags

In an IPW panel regression model featuring an interaction between LMC_Resin_Daily and LMC_Herb_Daily, the 169 positive and significant terms noted in eTable 59 are returned. Data are summarized in eTable 60. Fifty-five cannabis-related terms are extracted in eTable 61. The most significant terms for 50 unique anomalies are extracted and appear in eTable 62. P- and E-values are illustrated by anomaly in eFigs 47 and 48 and summary E- and P-values are illustrated in eFig. 48. This chart is again dominated by terms including daily cannabis exposure. Marginal effects are shown in eFig. 49. The mean and median marginal effect charts are again dominated by LMC_Herb_Daily.

These data are summarized in eTable 63 and are depicted graphically in the following eFigures. From eFig. 50, it is clear that cardiovascular and CNS anomalies dominate the picture. In considering marginal effects, it is clear that genital, face and limb anomalies dominate the total, mean and median marginal effects graphs.

Six Temporal Lags

When six temporal lags are considered in an IPW panel model again, featuring an interaction between LMC_Herb_Daily and LMC_Resin_Daily, 119 significant positive terms were returned (eTable 64). These are summarized in Table 7. Thirty cannabis-related terms may be extracted from these results (eTable 65). When the most significant of these are retained and listed in order of descending mEV, the results shown in Table 8 are revealed, which relate to 29 distinct CAs. The results in this table are remarkable for the very high mEVs noted, which decline from 3.61 × 1033 to 18.45. Genetic, cardiovascular, face and limb anomalies are notable. The table may be ordered in ascending order by P-values (eTable 66) or by marginal effects (eTable 67). Figure 14 illustrates the applicable P-value negative exponents, and Fig. 15 the mEV exponents. Summaries of the number of anomalies affected by substance, and E- and P-values are shown graphically in Fig. 16. Marginal effects are shown in Fig. 17. Here, the most efficacious terms are all related to compound daily indices of cannabis exposure, with the most salient term in each case being the interaction between LMC_Herb_Daily and LMC_Resin_Daily.

Table 7:

Summary table of multivariate interactive IPW panel regression results by substance at 6 years lag

Covariate	Number of positive terms	Sum mEV exponents	Sum P-value exponents	Mean % increment	Median % increment	Sum % increment
Alcohol	11	0	18	−6.43	−5.05	−70.76
Herb_THC	12	59	21	197.78	107.17	2373.35
Amphetamines	7	0	8	−15.53	−6.75	−108.74
Cocaine	12	0	15	−6.05	−4.91	−72.63
LM.Cann_Herb_THC_Daily.Int	10	69	22	469.61	200.11	4696.14
LM.Cann_Resin_THC_Daily.Int	10	21	16	10.31	4.58	103.05
LM.Cann_Herb_THC_Daily.Int x LM.Cann_Resin_THC_Daily.Int	8	115	14	1103.35	539.30	8826.77
Median_Income	27	0	50	0.00	0.00	0.04
Tobacco	22	0	67	4.3718	2.2772	96.1797

Table 8:

Multivariate interactive IPW panel regression results by substance at 6 years lag for most significant anomaly term

Anomaly	Term	Mean Anomaly Rate	ß-Estimate	Standard error	Sigma	T-statistic	Adj. R2	P-value	E-value estimate	E-value (lower bound)	P-value exponent	Lower E-value exponent	% Increment	P—Bonferroni	P—false discovery rate	P—Holm
Ear, face and neck	lag(LpmResinDailyInt, 6):lag(LpmHerbDailyInt, 6)	2.4050	124.8613	31.4481	0.7528	3.9704	0.2145	0.0005	7.09E + 65	3.61E + 33	3	33	198.54	0.0139	0.0028	0.0120
Hip dysplasia	lag(LpmResinDailyInt, 6):lag(LpmHerbDailyInt, 6)	6.3491	107.5064	34.7070	0.8308	3.0975	−0.3248	0.0045	2.76E + 51	1.40E + 19	2	19	75.21	0.1310	0.0131	0.0908
Hirschsprungs	lag(LpmHerbDailyInt, 6)	0.9900	21.0139	4.0517	0.2960	5.1864	0.5242	0.0000	2.30E + 28	6.01E + 17	4	17	457.40	0.0005	0.0005	0.0005
VATER/VACTERL	lag(LpmResinDailyInt, 6):lag(LpmHerbDailyInt, 6)	0.4798	66.1544	22.7293	0.5441	2.9105	−0.3275	0.0071	2.25E + 48	1.15E + 16	2	16	995.22	0.2072	0.0139	0.1260
PDA	lag(LpmHerbDailyInt, 6)	2.9086	59.4896	12.9845	0.9484	4.5816	0.0630	0.0001	1.23E + 25	3.21E + 14	4	14	155.69	0.0027	0.0009	0.0025
Anophthalmos	lag(LpmResinDailyInt, 6):lag(LpmHerbDailyInt, 6)	0.2059	22.2169	7.9565	0.1905	2.7923	−0.1420	0.0095	2.52E + 46	1.28E + 14	2	14	2319.08	0.2754	0.0155	0.1274
Conjoined twins	lag(LpmResinDailyInt, 6):lag(LpmHerbDailyInt, 6)	0.1128	20.6507	7.4064	0.1773	2.7882	−0.2137	0.0096	2.16E + 46	1.10E + 14	2	14	4233.14	0.2781	0.0155	0.1274
Congenital heart	lag(Herb, 6)	82.6649	10.2501	2.2265	0.1814	4.6037	0.4037	0.0001	4.35E + 22	1.40E + 13	4	13	5.48	0.0026	0.0009	0.0025
Single ventricle	lag(Herb, 6)	0.7525	24.7634	5.8151	0.4736	4.2585	0.1922	0.0002	9.19E + 20	2.97E + 11	3	11	601.77	0.0065	0.0016	0.0058
Skeletal dysplasias	lag(LpmHerbDailyInt, 6)	1.8518	30.4938	8.0752	0.5898	3.7762	−0.0234	0.0008	5.40E + 20	1.41E + 10	3	10	244.54	0.0231	0.0039	0.0191
Syndactyly	lag(LpmResinDailyInt, 6):lag(LpmHerbDailyInt, 6)	4.2302	49.3778	19.3916	0.4642	2.5464	−0.3173	0.0169	2.19E + 42	1.11E + 10	1	10	112.88	0.4902	0.0213	0.1403
Tot anom pul V ret	lag(Herb, 6)	0.6242	21.1817	5.7989	0.4723	3.6527	0.2061	0.0011	1.06E + 18	3.41E + 08	2	8	725.46	0.0319	0.0046	0.0253
Coarctation aorta	lag(Herb, 6)	3.3779	23.2510	7.0103	0.5710	3.3167	0.3213	0.0026	2.48E + 16	8.00E + 06	2	6	134.06	0.0756	0.0095	0.0574
Limb	lag(LpmResinDailyInt, 6):lag(LpmHerbDailyInt, 6)	37.7465	42.1716	18.0460	0.4320	2.3369	−0.3162	0.0271	7.63E + 38	3.88E + 06	1	6	12.65	0.7862	0.0302	0.1403
Anencephalus and ∼	lag(LpmHerbDailyInt, 6)	3.5488	14.5166	4.6598	0.3404	3.1153	0.6029	0.0043	1.43E + 17	3.74E + 06	2	6	127.60	0.1254	0.0131	0.0908
Hydrocephalus	lag(LpmHerbDailyInt, 6)	5.3737	15.9980	5.4132	0.3954	2.9554	0.3363	0.0064	1.95E + 16	5.10E + 05	2	5	84.27	0.1858	0.0139	0.1218
Cong. glaucoma	lag(LpmHerbDailyInt, 6)	0.2853	21.4190	7.3379	0.5360	2.9189	−0.0901	0.0070	1.24E + 16	3.24E + 05	2	5	1587.21	0.2031	0.0139	0.1260
Down syndrome	lag(LpmHerbDailyInt, 6)	20.7996	15.6869	5.3918	0.3938	2.9094	0.0850	0.0072	1.10E + 16	2.88E + 05	2	5	21.77	0.2078	0.0139	0.1260
VSD	lag(Herb, 6)	37.2628	8.3593	2.8717	0.2339	2.9109	0.2060	0.0071	2.66E + 14	8.60E + 04	2	4	12.15	0.2070	0.0139	0.1260
Edward syndrome	lag(Herb, 6)	4.9591	22.3653	7.9587	0.6483	2.8102	0.1626	0.0091	8.63E + 13	2.79E + 04	2	4	91.31	0.2639	0.0155	0.1274
ASD	lag(Herb, 6)	21.4739	20.3234	7.4426	0.6062	2.7307	−0.1801	0.0110	3.55E + 13	1.15E + 04	1	4	21.09	0.3189	0.0168	0.1274
Patau syndrome	lag(LpmHerbDailyInt, 6)	1.7507	20.3273	7.7899	0.5690	2.6094	0.2459	0.0146	2.63E + 14	6.86E + 03	1	3	258.66	0.4237	0.0202	0.1403
Eye	lag(Herb, 6)	3.6805	20.1616	7.6754	0.6252	2.6268	0.1091	0.0140	1.11E + 13	3.59E + 03	1	3	123.03	0.4069	0.0202	0.1403
Foetal alcohol	lag(LpmHerbDailyInt, 6)	0.2577	7.1157	2.7893	0.2037	2.5511	0.1258	0.0167	1.27E + 14	3.31E + 03	1	3	1757.20	0.4849	0.0213	0.1403
Aortic atresia ∼	lag(LpmResinDailyInt, 6):lag(LpmHerbDailyInt, 6)	0.5426	32.6821	15.3377	0.3672	2.1308	−0.1736	0.0424	3.02E + 35	1.54E + 03	1	3	880.06	1.0000	0.0424	0.1403
Respiratory	lag(Herb, 6)	3.2843	18.7230	7.6259	0.6211	2.4552	0.4014	0.0208	1.64E + 12	527.8115	1	2	137.88	0.6034	0.0251	0.1403
Digestive system	lag(Herb, 6)	17.0099	6.9501	2.9249	0.2382	2.3762	0.2960	0.0248	6.76E + 11	217.9977	1	2	26.62	0.7206	0.0288	0.1403
All anomalies	lag(Herb, 6)	249.8954	4.5283	2.0545	0.1673	2.2041	0.5142	0.0362	9.90E + 10	31.4602	1	1	1.81	1.0000	0.0389	0.1403
Choanal atresia	lag(Herb, 6)	0.9191	18.4092	8.5330	0.6950	2.1574	−0.0434	0.0400	5.87E + 10	18.4553	1	1	492.69	1.0000	0.0415	0.1403

Figure 14:

Scatterplot of negative log (base 10) P-value by congenital anomaly type from interactive multivariable IPW panel model temporally lagged by six years

Figure 15:

Scatterplot of log (base 10) mEV by congenital anomaly type from interactive multivariable IPW panel model temporally lagged by six years

Figure 16:

Summaries of E- and P-values by substance type for (A) number of anomalies with elevated mEV, (B) the sum of the mEV exponents and (C) the sum of the negative exponents of the significant P-values for the interactive multivariable IPW panel model temporally lagged by six years

Figure 17:

Summaries of marginal (overall) effects by substance type for (A) total percentage change at average marginal effect (AME), (B) the mean percentage change at AME and (C) the median percentage change at AME for the interactive multivariable IPW panel model temporally lagged by six years

Table 9 summarizes the organ systems implicated. Results are illustrated graphically in Fig. 18. Cardiovascular, general, limb and respiratory systems are all prominently represented. General, face and CNS effects are all prominently represented on the marginal effects summary shown in Fig. 19.

Table 9:

Summary table of multivariate interactive IPW panel regression results by organ system at 6 years lag

System	Number positive terms	System count	% Anomalies by system	Sum mEV exponents	Sum P-value exponents	Sum % increments	Mean % increments	Median % increments
Limb	3	6	50	35	4	200.74	66.91	75.21
Respiratory	1	2	50	2	1	137.88	137.88	137.88
General	5	11	45.45	44	9	7231.90	1446.38	995.22
Cardiovascular	9	23	39.13	68	21	2557.51	284.17	134.06
Central nervous	4	11	36.36	28	7	2653.98	663.50	125.32
Face	3	9	33.33	39	6	2278.44	759.48	492.69
Chromosomal	2	7	28.57	7	3	349.97	174.98	174.98
Gastrointestinal	2	8	25	19	5	484.03	242.01	242.01

Figure 18:

Summaries of E- and P-values by organ system. (A) Number of anomalies with elevated mEV. (B) Percentage of anomalies with elevated mEV. (C) The sum of the mEV exponents. (D) The sum of the negative exponents of the significant P-values for the interactive multivariable IPW panel model temporally lagged by six years

Figure 19:

Summaries of marginal (overall) effects by organ system for (A) total percentage change at average marginal effect (AME), (B) the mean percentage change at AME and (C) the median percentage change at AME for the interactive multivariable IPW panel model temporally lagged by six years

Discussion

Main Results

By using available European national datasets on congenital birth anomalies against various metrics of cannabis or other drug exposure, sophisticated analysis was able to associate all 90 tracked congenital birth anomalies with various metrics of cannabis exposure or combinations thereof. This is true both for birth defects considered in aggregate across the spectrum (Fig. 1 and accompanying analysis) and for CAs considered individually and by organ system. It was also shown convincingly both by random forest regression and by multivariable panel regression that compound indices of daily cannabis use were generally the most powerfully predictive covariates, confirming earlier concerns that it is the convergence of metrics of high cannabinoid exposure that most places infants and populations at genotoxic risk [15, 16].

Since all multivariable regression models were IPW, a pseudo-randomized analytical paradigm was created from which causal inferences could properly be drawn. Formal processes of causal inference were reinforced by the widespread use of E-values from both categorical two-by-two tables and linear and panel models.

Considered in aggregate whole, the CA rate was 71.8% higher in countries with high or increasing daily cannabis use, a result that was shown to be significant at the P = 4.74 × 10−17 level and was associated with an mEV of 1.52, which exceeds the epidemiological threshold for causality [81].

At bivariate analysis, 89 of 90 CAs were shown to be linked with the metrics of cannabinoid exposure, the sole exception being indeterminate sex. However, this anomaly was found to be cannabis related at multivariable panel regression in additive models and in interactive models lagged to zero and four years (Table 4 and eTables 34–35, 43–48, 59–62). For the series of substances tobacco, alcohol, amphetamines and cocaine, 3, 12, 23 and 68 anomalies had elevated mEV, respectively. For the metrics last month cannabis use, cannabis herb THC, cannabis resin THC, daily cannabis use interpolated, LMC_Herb_Daily and LMC_Resin_Daily, 23, 45, 34, 41, 41 and 42 mEVs were elevated.

Categorical Metrics

At categorical analysis, 84/90 CAs demonstrated elevated highest:lowest quintile ratios for LMC_Resin_Daily and 100% of CAs in the uronephrology, respiratory, limb, general, chromosomal and body wall classes were implicated. Eighty-two and 83 CAs had elevated mEVs for parameters including LMC_Resin_Daily and LMC_Herb_Daily, respectively (eTable 30 and eFig. 27).

The cannabinoid parameter with the highest variable importance on random forest and panel regression was LMC_Herb_Daily. Results for categorization on this covariate are shown in Table 2. It is of interest that the table is headed by VACTERL syndrome with highest to lowest quintile RR of 54.56 (17.55, 169.57), AFE 98.17% (94.30%, 99.41%) and PAR 97.76% (93.08%, 99.28%), P = 2.43 × 10−36 and mEV = 34.61.

Table 2:

Summary table of bivariate categorical relationships comparing the highest and lowest quintiles of last month cannabis × herb THC concentration × daily cannabis use interpolated

Anomaly	Total numbers of anomalies in quintile 5	Total numbers of normals in quintile 5	Total numbers of anomalies in quintile 1	Total numbers of normals in quintile 1	RR	RR (lower bound)	RR (upper bound)	AFE	AFE (lower bound)	AFE (upper bound)	PAR	PAR (lower bound)	PAR (upper bound)	Chi-square	P-value	E-value esimate	E-value (95% lower bound)
VATER/VACTERL	726	12 566 202	3	2 833 401	54.5626	17.5566	169.5705	0.9817	0.9430	0.9941	0.9776	0.9308	0.9928	157.1010	2.43E-36	108.62	34.61
Hip dysplasia	11 592	12 555 336	359	2 833 045	7.2802	6.5544	8.0864	0.8626	0.8474	0.8763	0.8367	0.8192	0.8525	1887.9451	0.00E + 00	14.04	12.59
Respiratory	6122	12 560 806	230	2 833 174	6.0013	5.2611	6.8456	0.8334	0.8099	0.8539	0.8032	0.7766	0.8266	924.2930	2.57E-203	11.48	10.00
Matern infect malform	1279	12 565 649	54	3 405 791	6.4191	4.8890	8.4280	0.8442	0.7955	0.8813	0.8100	0.7533	0.8537	237.0546	8.63E-54	12.32	9.25
Teratogenic synds	2284	12 564 644	111	3 405 734	5.5766	4.6093	6.7468	0.8207	0.7830	0.8518	0.7826	0.7393	0.8188	397.6438	8.97E-89	10.63	8.69
Edward syndrome	9689	12 557 239	676	3 405 169	3.8844	3.5931	4.1994	0.7426	0.7217	0.7619	0.6941	0.6710	0.7156	1354.3539	8.72E-297	7.23	6.65
Cystic lung	1413	12 565 516	86	3 405 759	4.4529	3.5818	5.5358	0.7754	0.7208	0.8194	0.7309	0.6697	0.7809	217.0697	1.97E-49	8.37	6.62
Foetal alcohol	784	12 566 144	45	3 405 800	4.7217	3.4964	6.3764	0.7882	0.7140	0.8432	0.7454	0.6618	0.8084	124.8482	2.75E-29	8.91	6.45
Bile duct A	600	12 566 328	34	3 405 811	4.7826	3.3854	6.7566	0.7909	0.7046	0.8520	0.7485	0.6512	0.8186	96.2678	5.02E-23	9.04	6.23
Mitral valve anomalies	2120	12 564 808	119	2 833 285	4.0167	3.3395	4.8312	0.7510	0.7006	0.7930	0.7111	0.6559	0.7575	255.3200	8.99E-58	7.50	6.13
Severe microcephaly	5083	12 561 845	381	3 405 464	3.6157	3.2582	4.0124	0.7234	0.6931	0.7508	0.6730	0.6397	0.7032	670.9093	3.17E-148	6.69	5.97
Double outlet RV	1768	12 565 160	105	2 833 298	3.7964	3.1180	4.6224	0.7366	0.6793	0.7837	0.6953	0.6331	0.7470	204.1802	1.28E-46	7.05	5.69
Turner syndrome	3364	12 563 564	260	3 405 585	3.5065	3.0909	3.9780	0.7148	0.6765	0.7486	0.6635	0.6217	0.7007	432.5232	2.29E-96	6.47	5.63
Cong. glaucoma	773	12 566 155	53	3 405 792	3.9528	2.9925	5.2211	0.7470	0.6658	0.8085	0.6991	0.6096	0.7681	109.4087	6.60E-26	7.37	5.43
Lateral anomalies	3182	12 563 747	213	2 833 191	3.3682	2.9320	3.8694	0.7031	0.6589	0.7416	0.6590	0.6117	0.7006	332.4899	1.38E-74	6.19	5.31
Genetic syndromes	9209	12 557 719	815	3 405 030	3.0623	2.8507	3.2897	0.6735	0.6492	0.6960	0.6187	0.5928	0.6430	1040.5508	1.38E-228	5.58	5.15
Hirschsprungs	1877	12 565 051	152	3 405 693	3.3467	2.8368	3.9482	0.7012	0.6475	0.7467	0.6487	0.5906	0.6985	231.4082	1.47E-52	6.15	5.12
Nervous system	43 028	12 523 900	3327	2 830 077	2.9159	2.8150	3.0205	0.6571	0.6448	0.6689	0.6099	0.5969	0.6225	3899.4204	0.00E + 00	5.28	5.08
Patau syndrome	3497	12 563 431	303	3 405 542	3.1279	2.7815	3.5174	0.6803	0.6405	0.7157	0.6260	0.5834	0.6643	403.7380	4.23E-90	5.71	5.01
Klinefelter	989	12 565 939	81	3 405 764	3.3091	2.6384	4.1503	0.6978	0.6210	0.7591	0.6450	0.5623	0.7120	120.6415	2.29E-28	6.07	4.72
Chromosomal	61 670	12 505 258	6213	3 399 632	2.6901	2.6209	2.7611	0.6283	0.6185	0.6378	0.5708	0.5605	0.5808	6018.9417	0.00E + 00	4.82	4.68
Anencephalus and ∼	5369	12 561 559	517	3 405 328	2.8145	2.5716	3.0803	0.6447	0.6111	0.6754	0.5881	0.5527	0.6206	551.8625	2.48E-122	5.07	4.58
Situs inversus	1252	12 565 676	109	3 405 736	3.1130	2.5596	3.7860	0.6788	0.6093	0.7359	0.6244	0.5503	0.6863	143.8198	1.95E-33	5.68	4.56
Craniosynostosis	5180	12 561 748	503	3 405 342	2.7910	2.5469	3.0585	0.6417	0.6074	0.6730	0.5849	0.5488	0.6181	527.1098	6.01E-117	5.03	4.53
Holoprosencephaly ∼	2849	12 564 079	274	3 405 571	2.8180	2.4894	3.1899	0.6451	0.5983	0.6865	0.5885	0.5393	0.6325	293.2214	4.94E-66	5.08	4.42
Small intestine S/A	1223	12 565 705	113	3 405 732	2.9332	2.4191	3.5566	0.6591	0.5866	0.7188	0.6033	0.5268	0.6675	131.8108	8.23E-31	5.31	4.27
An/microphthalmos	1558	12 565 370	159	3 405 686	2.6556	2.2558	3.1263	0.6234	0.5567	0.6801	0.5657	0.4964	0.6255	148.9350	1.48E-34	4.75	3.94
Down syndrome	34 645	12 532 283	4052	3 401 793	2.3172	2.2431	2.3938	0.5684	0.5542	0.5823	0.5089	0.4944	0.5230	2722.9178	0.00E + 00	4.06	3.91
Limb	61 519	12 505 409	6073	2 827 331	2.2839	2.2246	2.3449	0.5622	0.5505	0.5735	0.5117	0.4998	0.5232	4007.1357	0.00E + 00	4.00	3.88
Choanal atresia	1427	12 565 502	150	3 405 695	2.5783	2.1791	3.0506	0.6121	0.5411	0.6722	0.5539	0.4806	0.6169	131.1477	1.15E-30	4.60	3.78
Multicystic renal dys	6976	12 559 952	814	3 405 031	2.3226	2.1600	2.4975	0.5695	0.5370	0.5996	0.5099	0.4770	0.5408	549.2833	9.01E-122	4.08	3.74
Omphalocele	4382	12 562 547	503	3 405 342	2.3610	2.1529	2.5892	0.5765	0.5355	0.6138	0.5171	0.4754	0.5555	354.1107	2.70E-79	4.15	3.73
Skeletal dysplasias	2792	12 564 137	323	3 405 522	2.3427	2.0878	2.6286	0.5731	0.5210	0.6196	0.5137	0.4608	0.5614	222.8260	1.09E-50	4.12	3.59
Congenital heart	118 263	12 448 665	12 766	2 820 638	2.0887	2.0510	2.1271	0.5212	0.5124	0.5299	0.4704	0.4617	0.4791	6594.4726	0.00E + 00	3.60	3.52
Annular pancreas	274	12 566 654	24	3 405 821	3.0941	2.0386	4.6961	0.6768	0.5095	0.7871	0.6223	0.4457	0.7426	31.2763	1.12E-08	5.64	3.49
Cong. cataract	2204	12 564 724	259	3 405 586	2.3063	2.0277	2.6231	0.5664	0.5068	0.6188	0.5068	0.4466	0.5605	171.4996	1.74E-39	4.04	3.47
Hydrocephalus	7481	12 559 447	961	3 404 884	2.1098	1.9727	2.2563	0.5260	0.4931	0.5568	0.4661	0.4334	0.4970	497.3813	1.76E-110	3.64	3.36
Conjoined twins	202	12 566 726	17	3 405 828	3.2203	1.9631	5.2827	0.6895	0.4906	0.8107	0.6360	0.4253	0.7694	24.0045	4.81E-07	5.89	3.34
Amniotic band	960	12 565 968	112	3 405 733	2.3230	1.9101	2.8251	0.5695	0.4765	0.6460	0.5100	0.4162	0.5888	75.5777	1.76E-18	4.08	3.23
Digestive system	24 256	12 542 673	2761	2 830 643	1.9808	1.9043	2.0602	0.4951	0.4749	0.5146	0.4445	0.4246	0.4638	1205.8838	1.61E-264	3.37	3.22
Vascular disruptions	8534	12 558 394	1018	2 832 385	1.8901	1.7712	2.0170	0.4709	0.4354	0.5042	0.4207	0.3861	0.4534	381.4750	2.97E-85	3.19	2.94
Club foot	16 165	12 550 763	2412	3 403 433	1.8163	1.7403	1.8957	0.4494	0.4254	0.4725	0.3911	0.3680	0.4133	770.9346	5.63E-170	3.03	2.88
Neural tube defects	12 692	12 554 236	1891	3 403 954	1.8190	1.7332	1.9090	0.4502	0.4230	0.4762	0.3919	0.3658	0.4169	607.4520	2.00E-134	3.04	2.86
AVSD	5451	12 561 478	796	3 405 049	1.8559	1.7229	1.9992	0.4612	0.4196	0.4998	0.4024	0.3624	0.4400	274.2795	6.63E-62	3.12	2.84
Genital	29 065	12 537 863	3695	2 829 709	1.7735	1.7139	1.8352	0.4361	0.4165	0.4551	0.3870	0.3681	0.4053	1108.3316	2.55E-243	2.94	2.82
Hypoplastic left ht	3339	12 563 589	486	3 405 359	1.8620	1.6930	2.0479	0.4629	0.4093	0.5117	0.4041	0.3525	0.4516	169.3353	5.17E-39	3.13	2.78
Urinary	54 851	12 512 077	8623	3 397 222	1.7239	1.6853	1.7635	0.4199	0.4066	0.4329	0.3629	0.3503	0.3752	2274.3562	0.00E + 00	2.84	2.76
PV atresia	1449	12 565 479	203	3 405 642	1.9345	1.6702	2.2406	0.4831	0.4013	0.5537	0.4237	0.3445	0.4934	80.3873	1.54E-19	3.28	2.73
Abdominal wall defx	6953	12 559 976	1072	3 404 773	1.7578	1.6483	1.8746	0.4311	0.3933	0.4665	0.3735	0.3376	0.4075	303.5942	2.71E-68	2.91	2.68
Hydronephrosis	17 099	12 549 829	2725	3 403 120	1.7006	1.6332	1.7707	0.4120	0.3877	0.4353	0.3553	0.3325	0.3774	679.2268	4.92E-150	2.79	2.65
Hypospadias	24 229	12 542 699	3233	2 830 171	1.6897	1.6289	1.7528	0.4082	0.3861	0.4295	0.3601	0.3391	0.3805	804.4393	2.92E-177	2.77	2.64
Polydactyly	17 729	12 549 199	2904	3 402 941	1.6546	1.5909	1.7207	0.3956	0.3714	0.4189	0.3399	0.3173	0.3618	646.9932	5.03E-143	2.70	2.56
Diaphragmatic hernia	3749	12 563 179	590	3 405 255	1.7221	1.5789	1.8783	0.4193	0.3667	0.4676	0.3623	0.3126	0.4084	154.3957	9.49E-36	2.84	2.54
Severe CHD	28 086	12 538 842	4742	3 401 103	1.6052	1.5566	1.6553	0.3770	0.3576	0.3959	0.3226	0.3045	0.3401	927.5716	4.98E-204	2.59	2.49
TV S/A	851	12 566 077	123	3 405 722	1.8751	1.5521	2.2653	0.4667	0.3557	0.5586	0.4078	0.3014	0.4979	43.8914	1.74E-11	3.16	2.48
Bilat renal agenesis	1477	12 565 451	226	3 405 619	1.7712	1.5398	2.0373	0.4354	0.3506	0.5092	0.3776	0.2973	0.4488	65.8245	2.46E-16	2.94	2.45
Valproate syndrome	78	12 566 850	6	3 405 839	3.5232	1.5357	8.0827	0.7162	0.3489	0.8763	0.6650	0.2758	0.8451	10.0679	0.0008	6.50	2.44
Aortic atresia ∼	775	12 566 153	92	2 833 312	1.8993	1.5301	2.3575	0.4735	0.3465	0.5758	0.4232	0.3003	0.5246	35.0195	1.63E-09	3.21	2.43
Oesophageal S/A	3351	12 563 577	553	3 405 292	1.6423	1.5010	1.7968	0.3911	0.3338	0.4435	0.3357	0.2824	0.3851	119.2583	4.60E-28	2.67	2.37
Tetralogy of Fallot	4324	12 562 604	723	3 405 122	1.6209	1.4981	1.7536	0.3830	0.3325	0.4298	0.3282	0.2813	0.3720	147.3533	3.28E-34	2.62	2.36
All anomalies	358 462	12 208 466	64 516	3 341 329	1.5058	1.4934	1.5184	0.3359	0.3304	0.3414	0.2847	0.2796	0.2897	9542.4380	0.00E + 00	2.38	2.35
Blader extr/epispad	812	12 566 116	125	3 405 720	1.7605	1.4583	2.1253	0.4320	0.3143	0.5295	0.3744	0.2635	0.4686	35.5906	1.22E-09	2.92	2.28
Transpos grt vess	4214	12 562 714	733	3 405 112	1.5581	1.4405	1.6852	0.3582	0.3058	0.4066	0.3051	0.2571	0.3500	124.8471	2.75E-29	2.49	2.24
Hypoplastic right ht	891	12 566 037	142	3 405 703	1.7005	1.4245	2.0300	0.4119	0.2980	0.5074	0.3553	0.2489	0.4466	35.3481	1.38E-09	2.79	2.20
Limb reductions	7801	12 559 127	1431	3 404 414	1.4774	1.3965	1.5631	0.3231	0.2839	0.3602	0.2731	0.2376	0.3069	186.6616	8.51E-43	2.32	2.14
Spina bifida	5766	12 561 163	1075	3 404 770	1.4537	1.3620	1.5514	0.3121	0.2658	0.3554	0.2630	0.2215	0.3024	128.3342	4.74E-30	2.27	2.06
PDA	3491	12 563 437	639	3 405 206	1.4806	1.3609	1.6109	0.3246	0.2652	0.3792	0.2744	0.2208	0.3243	84.2905	2.14E-20	2.32	2.06
Anorectal S/A	3764	12 563 164	698	3 405 147	1.4615	1.3481	1.5844	0.3158	0.2582	0.3688	0.2664	0.2146	0.3147	85.8211	9.85E-21	2.28	2.03
PV stenosis	5594	12 561 334	1067	3 404 778	1.4209	1.3308	1.5170	0.2962	0.2486	0.3408	0.2488	0.2063	0.2889	111.7550	2.02E-26	2.19	1.99
Tot anom pul V ret	662	12 566 266	114	3 405 731	1.5738	1.2901	1.9198	0.3646	0.2249	0.4791	0.3110	0.1837	0.4185	20.3466	3.23E-06	2.52	1.90
Eye	6732	12 560 196	1334	3 404 511	1.3677	1.2897	1.4504	0.2688	0.2246	0.3105	0.2244	0.1854	0.2615	110.1067	4.64E-26	2.08	1.90
Single ventricle	1051	12 565 877	191	3 405 654	1.4913	1.2782	1.7399	0.3294	0.2177	0.4252	0.2788	0.1783	0.3670	26.1624	1.57E-07	2.35	1.87
Arterial truncus	841	12 566 088	150	3 405 695	1.5195	1.2772	1.8078	0.3419	0.2170	0.4468	0.2901	0.1774	0.3874	22.6105	9.92E-07	2.41	1.87
Encephalocele	1558	12 565 370	297	3 405 548	1.4217	1.2558	1.6095	0.2966	0.2037	0.3787	0.2491	0.1666	0.3234	31.2049	1.16E-08	2.20	1.82
VSD	55 036	12 511 892	11 869	3 393 976	1.2567	1.2321	1.2818	0.2043	0.1883	0.2199	0.1680	0.1544	0.1815	514.0649	4.14E-114	1.82	1.77
Ebstein’s anomaly	514	12 566 414	95	3 405 750	1.4663	1.1781	1.8251	0.3180	0.1512	0.4521	0.2684	0.1200	0.3918	11.8921	0.0003	2.29	1.64
Orofacial clefts	15 026	12 551 902	2772	2 830 632	1.2222	1.1737	1.2727	0.1818	0.1480	0.2143	0.1535	0.1240	0.1819	94.6207	1.15E-22	1.74	1.63
Posterior ureth valv	1806	12 565 122	389	3 405 456	1.2582	1.1277	1.4039	0.2052	0.1132	0.2877	0.1689	0.0905	0.2405	16.9660	1.90E-05	1.83	1.51
Ear, face and neck	3203	12 563 725	711	3 405 134	1.2209	1.1256	1.3242	0.1809	0.1116	0.2449	0.1481	0.0895	0.2029	23.2624	7.07E-07	1.74	1.50
Coarctation aorta	3788	12 563 141	874	3 404 971	1.1746	1.0913	1.2642	0.1487	0.0837	0.2090	0.1208	0.0666	0.1718	18.4386	8.77E-06	1.63	1.41
Cleft lip ± palate	8444	12 558 484	2027	3 403 818	1.1290	1.0756	1.1851	0.1143	0.0703	0.1562	0.0921	0.0560	0.1269	24.1059	4.56E-07	1.51	1.36
Cleft palate	6486	12 560 443	1286	2 832 118	1.1371	1.0711	1.2072	0.1206	0.0664	0.1717	0.1006	0.0546	0.1444	17.7599	1.25E-05	1.53	1.35
Gastroschisis	1953	12 564 975	453	3 405 392	1.1684	1.0549	1.2942	0.1441	0.0520	0.2273	0.1170	0.0406	0.1873	8.9283	0.0014	1.61	1.30
Anophthalmos	211	12 566 717	51	3 405 794	1.1213	0.8258	1.5224	0.1081	−0.2109	0.3431	0.0871	−0.1679	0.2864	0.5387	0.2315	1.49	1.00
Anotia	368	12 566 560	87	3 405 758	1.1464	0.9075	1.4481	0.1277	−0.1019	0.3094	0.1033	−0.0833	0.2577	1.3150	0.1257	1.56	1.00
Aortic valve S/A	1432	12 565 496	452	3 405 393	0.8586	0.7725	0.9544	−0.1647	−0.2946	−0.0478	−0.1252	−0.2193	−0.0383	7.9990	0.0023	1.60	—
ASD	17 175	12 549 753	5504	3 400 341	0.8457	0.8204	0.8717	−0.1825	−0.2189	−0.1471	−0.1382	−0.1646	−0.1123	117.5198	1.10E-27	1.65	—
Duodenal S/A	1390	12 565 539	451	3 405 394	0.8353	0.7511	0.9289	−0.1972	−0.3313	−0.0766	−0.1489	−0.2448	−0.0604	11.0618	0.0004	1.68	—
Indeterminate sex	421	12 566 507	161	3 405 684	0.7087	0.5910	0.8498	−0.4111	−0.6921	−0.1767	−0.2974	−0.4795	−0.1376	13.9470	9.40E-05	2.17	—
Syndactyly	5199	12 561 729	1713	3 404 132	0.8225	0.7788	0.8687	−0.2157	−0.2840	−0.1512	−0.1623	−0.2110	−0.1155	49.3507	1.07E-12	1.73	—

All anomalies are shown in this table as having RR of 1.51 (1.49, 1.52), AFE 33.59% (33.04%, 34.14%) and PAR 28.47% (27.96%, 28.97%), P < 2.2 × 10−307 and mEV = 2.35.

The CNS anomalies severe microcephaly, nervous system, anencephalus and similar hydrocephalus, neural tube defects, spina bifida and encephalocele were associated with AFEs of 72.34% (69.31%, 75.08%), 65.71% (64.48%, 66.89%), 64.47% (641.11%, 67.54%), 52.60% (49.31%, 55.68%), 52.60% (49.31%, 55.68%), 45.02% (42.30%, 47.62%), 31.21% (26.58%, 35.54%) and 29.66% (20.37%, 37.87%) and mEVs of 5.97, 5.08, 4.58, 3.36, 2.86, 2.06 and 1.82, respectively.

Facial anomalies have been shown to be developmentally related to abnormalities of CNS development since the head and the brain form due to the facial organizer and the forebrain organizer, which have interrelated and interconnected form and functions developmentally [57, 82]. Congenital glaucoma, holoprosencephaly, anophthalmos/microphthalmos, choanal atresia, congenital cataract, eye anomalies, orofacial clefts, ear, face and neck anomalies, cleft lip ± cleft palate and cleft palate were associated with AFEs of 74.70% (66.58%), 64.51% (59.83%, 68.65%), 62.34% (55.67%, 68.01%), 61.21% (54.11%, 67.22%), 56.64% (50.68%, 61.88%), 26.88% (22.46%,m 31.05%), 18.18% (14.80%, 21.43%), 18.09% (11.16%, 24.49%), 11.43% (7.03%,15.62%) and 12.06% (6.64%, 17.76%) and mEVs of 5.43, 4.42, 3.94, 3.78, 3.47, 1.90, 1.63, 1.50, 1.36 and 1.35, respectively.

The cardiovascular anomalies mitral valve anomalies, double outlet right ventricle, congenital heart disorders, pulmonary valve atresia, aortic atresia and similar, vascular disruptions, tricuspid valve stenosis or atresia, hypoplastic left heart, atrioventricular septal defect (AVSD), hypoplastic right heart, tetralogy of Fallot, severe congenital heart disease, total anomalous pulmonary venous return, transposition of the great vessels, arterial truncus, single ventricle, PDA, Ebstein’s anomaly, pulmonary valve stenosis, ventricular septal defect (VSD) and coarctation of the aorta were associated with AFEs of 75.10% (70.6%, 79.30%), 73.66% (67.93%, 78.37%), 52.12% (51.24%, 52.99%), 48.31% (40.13%, 55.37%), 47.35% (34.65%, 57.58%), 47.09% (43.54%, 50.42%), 46.67% (35.57%, 55.86%), 46.29% (40.93%, 51.17%), 46.12% (41.96%, 49.98%), 41.19% (29.80%, 503.74%), 38.30% (33.25%, 42.98%), 37.70% (35.76%, 39.59%), 36.46% (22.49%, 47.91%), 35.82% (30.58%, 40.66%), 34.19% (21.70%, 44.68%), 32.94% (21.77%, 42.52%), 32.46% (26.52%, 37.92%), 31.80% (15.12%, 45.12%), 29.62% (24.86%, 34.08%), 20.43% (18.83%, 21.99%) and 14.87% (8.37%, 20.90%) and mEVs of 6.13, 5.69, 3.52, 2.73, 2.43, 2.94, 2.48, 2.78, 2.84, 2.20, 2.36, 2.49, 1.90, 2.24, 1.87, 1.87, 2.06, 1.64, 1.99, 1.77, and 1.41, respectively.

The limb anomalies hip dysplasia, limb anomalies, skeletal dysplasia, polydactyly and limb reductions were associated with AFEs of 86.26% (84.74%, 87.63%), 56.22% (55.05%, 57.35%), 57.31% (52.10%, 61.96%), 39.56% (37.14%, 41.89%) and 32.31% (28.39%, 36.02%) and mEVs of 12.59, 3.83, 3.59, 2.56 and 2.14, respectively.

The genetic syndromes Trisomy 18 (Edwards syndrome), Turner syndrome (monosomy X), genetic syndromes, Trisomy 13 (Patau syndrome), Klinefelter (male disomy X), chromosomal anomalies and Trisomy 21 (Downs syndrome) were noted to have AFEs of 74.26% (72.17%, 76.19%), 71.48% (67.65%, 74.86%), 67.35% (64.92%, 69.60%), 68.03% (64.05%, 71.57%), 69.78% (62.10%, 75.91%), 6.83% (61.85%, 63.78%) and 56.84% (55.42%, 58.23%) and mEVs of 6.65, 5.63, 5.15, 5.01, 4.72, 4.68 and 3.91, respectively.

Considering LMC_Resin_Daily and comparing the highest and lowest quintiles for all anomalies, genetic disorders, central nervous, cardiovascular, facial, limb and VACTERL disorders, the following notable observations were reported (Table 3). For the sequence of disorders all anomalies, teratogenic syndromes, Trisomy 18, Trisomy 13, genetic syndromes, Turners syndrome, Trisomy 21 and chromosomal disorders, the applicable P-values were <2.2 × 10−307, 3.69 × 10−79, 6.31 × 10−252, 1.50 × 10−93, 7.42 × 10−227, 8.52 × 10−58, 9.84 × 10−318 and <2.23 × 10−307, with associated mEVs of 2.49, 5.55, 4.46, 4.29, 4.01, 3.18, 2.61 and 3.18, respectively. For the CNS disorders severe microcephaly and anencephalus and similar, the P-values were 9.68 × 10−101 and 1.57 × 10−101 and the mEVs were 3.57 and 3.44, respectively. For the cardiovascular disorders, congenital heart disease, severe congenital heart disease, VSD, AVSD, vascular disruptions, aortic atresia and similar, and tetralogy of Fallot and Ebstein’s anomaly, the applicable P-values were <2.2 × 10−307, 1.79 × 10−232, <2.2 × 10−307, 1.65 × 10−101, 6.50 × 10−65, 2.45 × 10−45, 1.54 × 10−34 and 3.23 × 10−17, respectively, and the relevant mEVs were 2.08, 2.46, 2.74, 3.38, 2.34, 2.94, 2.25 and 4.09, respectively. For the facial anomalies holoprosencephaly and orofacial clefts, the applicable P-values were 6.11 × 10−44 and 6.07 × 10−51, respectively, and the applicable mEVs were 2.90 and 1.86, respectively. For the limb anomalies limb anomalies and limb reduction anomalies, the P-values were <2.3 × 10−307 and 1.94 × 10−54, respectively, and the relevant mEVs were 3.49 and 2.25, respectively. For VACTERL, the relevant P-value was 1.45 × 10−32 and mEV was 6.00.

Table 3:

Summary table of bivariate categorical relationships comparing the highest and lowest quintiles oflast month cannabis × resin THC concentration × daily cannabis use interpolated

Anomaly	Total numbers of anomalies in quintile 5	Total numbers of normals in quintile 5	Total numbers of anomalies in Quintile 1	Total numbers of normals in quintile 1	RR	RR (Lower bound)	RR (upper bound)	AFE	AFE (lower bound)	AFE (upper bound)	PAR	PAR (lower bound)	PAR (upper bound)	Chi-squared	P-value	E-value esimate	E-value (95% lower bound)
Foetal alcohol	686	10 335 677	5	3 967 367	52.6610	21.8492	126.9237	0.9810	0.9542	0.9921	0.9739	0.9375	0.9891	251.5786	5.88E-57	104.82	43.19
Amniotic band	862	10 335 501	59	3 967 313	5.6078	4.3077	7.3002	0.8217	0.7679	0.8630	0.7690	0.7044	0.8196	209.0833	1.09E-47	10.69	8.08
Multicystic renal dys	6362	10 330 000	684	3 966 688	3.5700	3.2993	3.8630	0.7199	0.6969	0.7411	0.6500	0.6242	0.6741	1143.2133	6.70E-251	6.60	6.05
VATER/VACTERL	726	10 335 637	61	3 690 483	4.2494	3.2723	5.5182	0.7647	0.6944	0.8188	0.7054	0.6251	0.7685	139.8292	1.45E-32	7.97	6.00
Hypoplastic right ht	961	10 335 402	97	3 967 276	3.8027	3.0861	4.6856	0.7370	0.6760	0.7866	0.6695	0.6004	0.7266	182.0095	8.82E-42	7.07	5.62
Teratogenic synds	2000	10 334 363	219	3 967 154	3.5053	3.0489	4.0300	0.7147	0.6720	0.7519	0.6442	0.5965	0.6862	353.4852	3.69E-79	6.47	5.55
Hip dysplasia	11 179	10 325 184	1325	3 689 220	3.0124	2.8457	3.1889	0.6680	0.6486	0.6864	0.5972	0.5762	0.6172	1593.9238	0.00E + 00	5.47	5.14
Omphalocele	4163	10 332 199	562	3 966 810	2.8432	2.6035	3.1050	0.6483	0.6159	0.6779	0.5712	0.5366	0.6032	591.8580	4.94E-131	5.13	4.65
Valproate syndrome	72	10 336 290	4	3 967 368	6.9089	2.5243	18.9091	0.8553	0.6039	0.9471	0.8102	0.5075	0.9269	19.1507	6.04E-06	13.30	4.49
Edward syndrome	8773	10 327 590	1265	3 966 107	2.6619	2.5096	2.8235	0.6243	0.6015	0.6458	0.5457	0.5216	0.5685	1147.9334	6.31E-252	4.77	4.46
Matern infect malform	1134	10 335 229	147	3 967 225	2.9609	2.4935	3.5160	0.6623	0.5990	0.7156	0.5863	0.5183	0.6446	169.0136	6.08E-39	5.37	4.42
Patau syndrome	3186	10 333 177	457	3 966 915	2.6759	2.4260	2.9515	0.6263	0.5878	0.6612	0.5477	0.5072	0.5849	419.5932	1.50E-93	4.79	4.29
Ebstein’s anomaly	449	10 335 913	56	3 967 316	3.0775	2.3311	4.0628	0.6751	0.5710	0.7539	0.6002	0.4882	0.6877	69.8286	3.23E-17	5.61	4.09
Lateral anomalies	2900	10 333 462	405	3 690 140	2.5566	2.3042	2.8367	0.6089	0.5660	0.6475	0.5342	0.4898	0.5749	336.8852	1.52E-75	4.55	4.04
Genetic syndromes	9047	10 327 316	1435	3 965 937	2.4198	2.2888	2.5584	0.5868	0.5631	0.6091	0.5064	0.4821	0.5296	1032.5849	7.42E-227	4.27	4.01
Cystic lung	1210	10 335 153	176	3 967 196	2.6388	2.2529	3.0908	0.6210	0.5561	0.6765	0.5422	0.4744	0.6012	156.3959	3.47E-36	4.72	3.93
Situs inversus	1098	10 335 264	159	3 967 214	2.6506	2.2445	3.1302	0.6227	0.5545	0.6805	0.5440	0.4727	0.6056	142.7689	3.30E-33	4.74	3.92
Hydronephrosis	15 960	10 320 403	2758	3 964 614	2.2211	2.1332	2.3127	0.5498	0.5312	0.5676	0.4688	0.4502	0.4868	1580.8321	0.00E + 00	3.87	3.69
Posterior urethral valve	1700	10 334 662	277	3 967 096	2.3556	2.0747	2.6746	0.5755	0.5180	0.6261	0.4949	0.4366	0.5471	185.8450	1.28E-42	4.14	3.57
Severe microcephaly	4513	10 331 849	774	3 966 599	2.2380	2.0737	2.4153	0.5532	0.5178	0.5860	0.4722	0.4367	0.5054	452.6242	9.68E-101	3.90	3.57
Limb	55 575	10 280 787	9536	3 681 008	2.0808	2.0362	2.1265	0.5194	0.5089	0.5297	0.4433	0.4329	0.4536	4590.7979	0.00E + 00	3.58	3.49
Anencephalus and ∼	4844	10 331 519	860	3 966 512	2.1619	2.0107	2.3245	0.5374	0.5027	0.5698	0.4564	0.4219	0.4889	456.2601	1.57E-101	3.75	3.44
AVSD	4990	10 331 373	900	3 966 472	2.1281	1.9823	2.2846	0.5301	0.4955	0.5623	0.4491	0.4150	0.4812	456.1556	1.65E-101	3.68	3.38
Polydactyly	15 794	10 320 569	2964	3 964 408	2.0453	1.9666	2.1271	0.5111	0.4915	0.5299	0.4303	0.4112	0.4488	1334.9209	1.46E-292	3.51	3.35
Abdominal wall defx	6525	10 329 837	1234	3 966 138	2.0296	1.9098	2.1568	0.5073	0.4764	0.5364	0.4266	0.3965	0.4552	542.2794	3.01E-120	3.48	3.23
Urinary	48 201	10 288 162	9559	3 957 814	1.9354	1.8935	1.9783	0.4833	0.4719	0.4945	0.4033	0.3923	0.4141	3621.1633	0.00E + 00	3.28	3.19
Turner syndrome	2965	10 333 397	551	3 966 822	2.0654	1.8859	2.2620	0.5158	0.4698	0.5579	0.4350	0.3900	0.4767	255.4275	8.52E-58	3.55	3.18
Chromosomal	53 757	10 282 606	10 726	3 956 646	1.9237	1.8843	1.9639	0.4802	0.4693	0.4908	0.4003	0.3899	0.4106	3983.8113	0.00E + 00	3.26	3.18
Blader extr/epispad	809	10 335 554	144	3 967 229	2.1564	1.8061	2.5746	0.5363	0.4463	0.6116	0.4552	0.3668	0.5313	75.8075	1.56E-18	3.74	3.01
Club foot	14 653	10 321 710	2999	3 964 373	1.8754	1.8032	1.9505	0.4668	0.4454	0.4873	0.3875	0.3672	0.4071	1018.4416	8.81E-224	3.16	3.01
Hypoplastic left ht	2946	10 333 417	580	3 966 793	1.9496	1.7835	2.1311	0.4871	0.4393	0.5308	0.4069	0.3612	0.4495	224.1839	5.53E-51	3.31	2.97
Aortic atresia ∼	695	10 335 668	115	3 690 429	2.1578	1.7714	2.6284	0.5366	0.4355	0.6195	0.4604	0.3608	0.5444	61.3019	2.45E-15	3.74	2.94
Conjoined twins	190	10 336 173	28	3 967 345	2.6045	1.7516	3.8729	0.6161	0.4291	0.7418	0.5369	0.3456	0.6723	24.1230	4.52E-07	4.65	2.90
Holoprosencephaly ∼	2601	10 333 762	519	3 966 853	1.9236	1.7506	2.1136	0.4801	0.4288	0.5269	0.4003	0.3513	0.4456	191.9022	6.11E-44	3.26	2.90
Respiratory	5185	10 331 178	1028	3 689 517	1.8009	1.6843	1.9255	0.4447	0.4063	0.4806	0.3711	0.3350	0.4053	305.6756	9.56E-69	3.00	2.76
VSD	46 424	10 289 939	10 401	3 956 971	1.7132	1.6772	1.7499	0.4163	0.4038	0.4286	0.3401	0.3285	0.3514	2532.8001	0.00E + 00	2.82	2.74
Neural tube defects	11 338	10 325 024	2487	3 964 885	1.7498	1.6755	1.8274	0.4285	0.4032	0.4528	0.3514	0.3279	0.3741	655.9895	5.57E-145	2.90	2.74
Bilat renal agenesis	1487	10 334 875	306	3 967 066	1.8652	1.6493	2.1094	0.4639	0.3937	0.5259	0.3847	0.3186	0.4444	101.8620	2.98E-24	3.14	2.68
PV stenosis	5118	10 331 245	1128	3 966 244	1.7415	1.6328	1.8575	0.4258	0.3876	0.4616	0.3489	0.3136	0.3824	291.9487	9.35E-66	2.88	2.65
Choanal atresia	1240	10 335 123	256	3 967 116	1.8592	1.6251	2.1269	0.4621	0.3847	0.5298	0.3830	0.3103	0.4481	84.2572	2.17E-20	3.12	2.63
Down syndrome	29 415	10 306 948	6805	3 960 567	1.6591	1.6160	1.7034	0.3973	0.3812	0.4129	0.3226	0.3080	0.3370	1450.7519	9.84E-318	2.70	2.61
Nervous system	36 794	10 299 568	8033	3 682 511	1.6354	1.5964	1.6753	0.3885	0.3736	0.4031	0.3189	0.3053	0.3322	1632.9280	0.00E + 00	2.65	2.57
Double outlet RV	1617	10 334 745	325	3 690 219	1.7764	1.5769	2.0012	0.4371	0.3658	0.5003	0.3639	0.2976	0.4240	91.8471	4.68E-22	2.95	2.53
Encephalocele	1438	10 334 925	310	3 967 062	1.7805	1.5748	2.0129	0.4383	0.3650	0.5032	0.3606	0.2927	0.4220	87.2573	4.76E-21	2.96	2.53
Skeletal dysplasias	2380	10 333 983	531	3 966 841	1.7204	1.5659	1.8900	0.4187	0.3614	0.4709	0.3423	0.2898	0.3910	130.9758	1.25E-30	2.83	2.51
Bile duct A	509	10 335 854	101	3 967 272	1.9343	1.5625	2.3947	0.4830	0.3600	0.5824	0.4030	0.2866	0.5005	38.0366	3.47E-10	3.28	2.50
TV S/A	767	10 335 596	159	3 967 214	1.8515	1.5609	2.1964	0.4599	0.3593	0.5447	0.3809	0.2869	0.4626	51.5808	3.44E-13	3.11	2.50
All anomalies	309 546	10 026 817	75 745	3 891 627	1.5686	1.5563	1.5810	0.3625	0.3574	0.3675	0.2912	0.2867	0.2957	12 889.1510	0.00E + 00	2.51	2.49
Severe CHD	25 055	10 311 308	6065	3 961 307	1.5856	1.5418	1.6307	0.3693	0.3514	0.3868	0.2974	0.2813	0.3130	1058.4314	1.79E-232	2.55	2.46
Anorectal S/A	3430	10 332 933	792	3 966 581	1.6623	1.5387	1.7958	0.3984	0.3501	0.4431	0.3237	0.2799	0.3648	169.8277	4.03E-39	2.71	2.45
Craniosynostosis	4319	10 332 043	1018	3 966 355	1.6284	1.5210	1.7435	0.3859	0.3425	0.4264	0.3123	0.2732	0.3493	199.8699	1.11E-45	2.64	2.41
Cong. cataract	1905	10 334 457	435	3 966 938	1.6809	1.5146	1.8654	0.4051	0.3398	0.4639	0.3298	0.2705	0.3843	97.6924	2.44E-23	2.75	2.40
Eye	5774	10 330 589	1380	3 965 992	1.6060	1.5144	1.7031	0.3773	0.3397	0.4128	0.3045	0.2708	0.3367	254.7824	1.18E-57	2.59	2.40
Vascular disruptions	7593	10 328 770	1732	3 688 812	1.5653	1.4857	1.6491	0.3611	0.3269	0.3936	0.2941	0.2634	0.3234	288.0838	6.50E-65	2.51	2.34
Cong. glaucoma	650	10 335 713	140	3 967 232	1.7821	1.4846	2.1391	0.4388	0.3264	0.5325	0.3611	0.2575	0.4502	39.5357	1.61E-10	2.96	2.33
PV atresia	1435	10 334 928	331	3 967 042	1.6640	1.4766	1.8753	0.3990	0.3228	0.4667	0.3243	0.2553	0.3868	71.2770	1.55E-17	2.72	2.32
Anotia	298	10 336 065	59	3 967 313	1.9387	1.4663	2.5632	0.4842	0.3180	0.6099	0.4042	0.2477	0.5281	22.3830	1.12E-06	3.29	2.29
Hirschsprungs	1562	10 334 800	366	3 967 006	1.6381	1.4619	1.8355	0.3895	0.3159	0.4552	0.3156	0.2495	0.3759	73.7104	4.52E-18	2.66	2.28
Tetralogy of Fallot	3815	10 332 548	943	3 966 430	1.5528	1.4460	1.6675	0.3560	0.3084	0.4003	0.2854	0.2434	0.3251	148.8534	1.54E-34	2.48	2.25
Limb reductions	6606	10 329 757	1663	3 965 710	1.5247	1.4449	1.6089	0.3441	0.3079	0.3785	0.2749	0.2431	0.3054	240.0232	1.94E-54	2.42	2.25
Genital	24 721	10 311 642	6108	3 684 437	1.4451	1.4052	1.4861	0.3080	0.2884	0.3271	0.2470	0.2299	0.2637	672.8740	1.18E-148	2.25	2.16
Hypospadias	20 687	10 315 676	5119	3 685 426	1.4429	1.3995	1.4877	0.3069	0.2854	0.3278	0.2461	0.2274	0.2643	558.8957	7.31E-124	2.24	2.15
Spina bifida	5056	10 331 306	1317	3 966 056	1.4735	1.3868	1.5656	0.3214	0.2789	0.3613	0.2549	0.2182	0.2899	159.0619	9.07E-37	2.31	2.12
Diaphragmatic hernia	3337	10 333 026	866	3 966 507	1.4790	1.3725	1.5938	0.3239	0.2714	0.3726	0.2571	0.2117	0.2999	106.7022	2.59E-25	2.32	2.09
Congenital heart	96 438	10 239 925	24 784	3 665 761	1.3893	1.3701	1.4088	0.2802	0.2701	0.2902	0.2229	0.2143	0.2315	2169.7116	0.00E + 00	2.12	2.08
Annular pancreas	286	10 336 077	63	3 967 309	1.7425	1.3265	2.2889	0.4261	0.2461	0.5631	0.3492	0.1862	0.4795	16.3330	2.66E-05	2.88	1.98
Cleft lip ± palate	7984	10 328 379	2236	3 965 136	1.3705	1.3077	1.4363	0.2703	0.2353	0.3038	0.2112	0.1818	0.2396	175.0977	2.85E-40	2.08	1.94
Transpos grt vess	3798	10 332 564	1042	3 966 330	1.3990	1.3063	1.4983	0.2852	0.2345	0.3326	0.2238	0.1809	0.2644	93.0853	2.50E-22	2.15	1.94
Oesophageal S/A	2927	10 333 436	814	3 966 559	1.3802	1.2771	1.4916	0.2755	0.2169	0.3296	0.2155	0.1664	0.2618	66.7118	1.57E-16	2.10	1.87
Hydrocephalus	6444	10 329 919	1841	3 965 531	1.3435	1.2757	1.4149	0.2557	0.2161	0.2932	0.1989	0.1659	0.2305	125.8285	1.68E-29	2.02	1.87
Orofacial clefts	13 773	10 322 590	3735	3 686 809	1.3166	1.2699	1.3651	0.2405	0.2125	0.2674	0.1892	0.1658	0.2119	223.9970	6.07E-51	1.96	1.86
Gastroschisis	1815	10 334 548	497	3 966 875	1.4017	1.2693	1.5479	0.2866	0.2122	0.3540	0.2250	0.1622	0.2831	44.9227	1.02E-11	2.15	1.85
Arterial truncus	750	10 335 613	195	3 967 177	1.4763	1.2611	1.7282	0.3226	0.2070	0.4213	0.2560	0.1570	0.3435	23.7801	5.40E-07	2.31	1.83
Digestive system	20 789	10 315 573	5754	3 684 791	1.2900	1.2529	1.3282	0.2248	0.2019	0.2471	0.1761	0.1570	0.1947	294.3454	2.81E-66	1.90	1.82
Single ventricle	913	10 335 450	245	3 967 127	1.4303	1.2422	1.6470	0.3009	0.1950	0.3928	0.2372	0.1475	0.3175	25.0121	2.85E-07	2.21	1.79
Syndactyly	4503	10 331 859	1346	3 966 027	1.2841	1.2082	1.3647	0.2212	0.1724	0.2672	0.1703	0.1305	0.2083	65.1528	3.47E-16	1.89	1.71
Cleft palate	5688	10 330 674	1609	3 688 936	1.2622	1.1943	1.3340	0.2077	0.1627	0.2504	0.1619	0.1250	0.1973	68.3469	6.86E-17	1.84	1.68
An/microphthalmos	1347	10 335 016	401	3 966 972	1.2893	1.1533	1.4414	0.2244	0.1329	0.3062	0.1729	0.0987	0.2410	20.0633	3.75E-06	1.90	1.57
Aortic valve S/A	1234	10 335 128	388	3 966 985	1.2207	1.0891	1.3682	0.1808	0.0818	0.2691	0.1376	0.0594	0.2093	11.7826	2.99E-04	1.74	1.40
Mitral valve anomalies	1674	10 334 689	520	3 690 024	1.1494	1.0417	1.2683	0.1300	0.0400	0.2115	0.0992	0.0290	0.1643	7.7067	0.0028	1.56	1.25
PDA	2871	10 333 491	1012	3 966 360	1.0889	1.0136	1.1698	0.0816	0.0134	0.1451	0.0604	0.0092	0.1088	5.4322	0.0099	1.40	1.13
Coarctation aorta	3241	10 333 122	1161	3 966 212	1.0715	1.0020	1.1458	0.0667	0.0020	0.1272	0.0491	0.0010	0.0949	4.0768	0.0217	1.35	1.05
Duodenal S/A	1275	10 335 088	473	3 966 900	1.0346	0.9310	1.1498	0.0335	−0.0741	0.1303	0.0244	−0.0536	0.0967	0.3999	0.2636	1.22	1.00
Klinefelter	798	10 335 564	279	3 967 094	1.0978	0.9579	1.2582	0.0891	−0.0439	0.2052	0.0660	−0.0332	0.1558	1.8022	0.0897	1.43	1.00
Small intestine S/A	1053	10 335 310	394	3 966 978	1.0258	0.9137	1.1517	0.0252	−0.0945	0.1317	0.0183	−0.0680	0.0976	0.1862	0.3330	1.19	1.00
Anophthalmos	191	10 336 172	78	3 967 295	0.9399	0.7223	1.2231	−0.0640	−0.3845	0.1824	−0.0454	−0.2604	0.1329	0.2130	0.3222	1.32	—
ASD	13 356	10 323 007	6183	3 961 190	0.8291	0.8045	0.8545	−0.2061	−0.2430	−0.1703	−0.1409	−0.1646	−0.1176	149.0682	1.39E-34	1.70	—
Ear, face and neck	2475	10 333 888	1232	3 966 140	0.7711	0.7202	0.8256	−0.2969	−0.3886	−0.2112	−0.1982	−0.2541	−0.1448	55.9168	3.78E-14	1.92	—
Indeterminate sex	303	10 336 060	241	3 967 132	0.4826	0.4075	0.5715	−1.0722	−1.4542	−0.7497	−0.5972	−0.7550	−0.4536	74.4764	3.07E-18	3.56	—
Tot anom pul V ret	526	10 335 837	216	3 967 156	0.9347	0.7978	1.0951	−0.0699	−0.2535	0.0868	−0.0495	−0.1742	0.0619	0.6988	0.2016	1.34	—

Panel Regression

At additive IPW panel regression, 74 CAs were shown to be cannabis related. The introduction of interactive terms in the regression had the effect of increasing this number somewhat to 76 CAs. The number of implicated CAs dropped to 31, 50 and 29 CAs in interactive panel models lagged by two, four and six years.

In additive panel models, the most affected organ systems were uronephrology, face, central nervous body wall and cardiovascular systems, but 50% or more of CAs in all organ systems studied had elevated mEVs. Chromosomal, face, central nervous and cardiovascular systems had highest cumulative mEVs.

This pattern persisted with the introduction of interactive terms and temporal lagging. At six years of temporal lag in an interactive IPW panel, model body wall, genital, CNS and cardiovascular system had the highest percentage of CAs implicated. Cardiovascular, uronephrology, central nervous and general classes had the highest cumulative mEVs.

Interpretation

Several features stand out as major implications from the substantial body of evidence presented, namely the strength of reported effects (documented in the tables); the breadth of reported affects comprehensively across the spectrum of CAs and anomaly classes; the consistency with in vitro mechanistic studies (discussed below); the consistency with animal studies [83-85]; the consistency with results reported from elsewhere [4–9, 18, 55, 86–88]; concordance with recent cancer analyses [17, 18, 20, 89] and the implications of presumptive cannabinoid genotoxicity for cellular ageing community-wide [19].

A critical concern is that increasing European cannabis exposure is associated with a broad spectrum of CAs, in fact extending to the totality of anomalies reported by European data. This finding is compounded when one notes that cannabis exposure has also been linked with breast, thyroid, liver and pancreatic cancer and acute lymphoid leukaemia [17] and acute myeloid leukaemia of childhood [18], which latter constitute heritable mutagenic and teratogenic malignant syndromes.

However, cancer and CAs are relatively rare disorders compared to the widespread availability in the community of a known genotoxin, including allowing its entry into the food chain. Moreover, as our results clearly indicate that it is the convergence of increased cannabis prevalence, intensity and concentration of exposure driving the expected total dose–response relationships, which is the leading risk factor for heritable genotoxic outcomes, it is likely that these changes will leverage multiplicatively off each other in populations where policies favouring cannabis liberalization are in play.

Since cannabis is a known major epigenomic toxin [31–37, 90] and since DNA methylation has been shown to be a direct cause of mammalian ageing [19], it becomes difficult to avoid the conclusion that increasing cannabis use will actually increase the epigenomic age of widely exposed populations.

Furthermore, since many epigenomic effects are known to be heritable to several subsequent generations, this implies the multigenerational transmission of these changes in addition to described diagnosable genotoxic [18] and/or neurotoxic [91-94] outcomes. It also includes the disconcerting possibilities that babies will be born ‘old’ with advanced epigenomic ages and may continue to age in an accelerated manner as has been demonstrated in clinical populations [95] or develop abnormally due to large-scale epigenomic dysregulation.

Commentary on Specific Defects by Systems

Total Anomalies

All anomalies were significantly elevated on bivariate continuous (P = 0.0006, mEV = 8.41, eTable 12) and categorical [RR = 1.57 (1.56, 1.58), P < 2.2 × 10−307 and mEV = 2.49; Table 3] analyses. All anomalies were significantly associated with cannabis exposure on additive IPW panel regression (P = 0.0022, mEV = 451.56.41; Table 4) and in interactive IPW panel models lagged to zero (P = 0.0067, mEV = 7.60 × 1012; eTable 46), four (P = 4.20 × 10−6, mEV = 1.11 × 104; eTable 62) and six (P = 0.0362, mEV = 31.46; Table 8) years.

This important finding has been replicated in Colorado and Canada [5, 7]. No comparable metric exists for this datapoint in USA datasets generally.

Cardiovascular System

Cardiovascular disorders are widely acknowledged to be the commonest CAs. Cardiovascular disorders associated with metrics of cannabis use on bivariate analysis were: aortic atresia and similar, aortic valve stenosis or atresia, arterial truncus, ASD, AVSD, coarctation of aorta, congenital heart disease, double outlet right ventricle, Ebstein’s anomaly, hypoplastic left heart syndrome, hypoplastic right heart syndrome, mitral valve anomalies, PDA, pulmonary valve atresia, pulmonary valve stenosis, severe congenital heart disease, single ventricle, tetralogy of Fallot, total anomalous pulmonary venous return, transposition of the great vessels, tricuspid valve stenosis or atresia, vascular disruptions and VSD.

It is thus of interest that total cardiovascular disorders were noted to also be increased in association with cannabis use in the northern Canadian Provinces and in the US state of Colorado [5, 7]. Similarly, it is important to note that ASD was noted to occur with elevated rates in the US states of Hawaii and Colorado, across the USA and in Australia, (and presumably also in Canada as it is the commonest cardiovascular anomaly) [4–7, 55]; VSD was noted to be increased in Australia, USA Colorado, Hawaii and in other series [4–7, 55, 86]; PDA was noted to be increased in Colorado, Australia and USA [5, 8, 55]; tetralogy of Fallot was increased in Hawaii, Australia and USA [4, 8, 55, 96]; AVSD was noted to be elevated also in USA (manuscript submitted) and pulmonary valve atresia and stenosis was associated in USA [55].

Chromosomal Disorders

The chromosomal disorders identified in this study as being linked to indices of cannabis exposure were: Chromosomal, Downs syndrome (Trisomy 21), Edward syndrome (Trisomy 18), Genetic syndromes, Klinefelter (Male XXY), Patau syndrome (Trisomy 13) and Turner syndrome (Female XO). It is of considerable importance that elevated rates of Downs syndrome have previously been identified in association with cannabis exposure in Colorado [5], Australia [8], Canada [7], Hawaii [4] and USA [11, 55]; of Edwards syndrome in USA (manuscript submitted); of Patau syndrome in USA [18, 55]; of Turner syndrome in USA [18, 55]; of Turner syndrome in Australia [8] and of chromosomal anomalies in USA, Canada and Australia [7, 8, 18, 55].

In this context, it becomes important to observe that cannabis use has also been linked with testicular cancer in several studies [97-102] and with acute lymphoid leukaemia [17, 18, 20]. These disorders invariably or generally involve major rearrangements, translocations or deletions of chromosomes 12 and 19, respectively [100, 102–111].

If one aggregates all of these chromosomal disorders together, noting that chromosomes 12, 13, 18, 9, 21 and X are of 133, 114, 80, 57, 48 and 153 MB in length, it becomes clear that this provides clinical evidence of direct cannabinoid genotoxicity to 585 MB of the 3000 MB of the human DNA or 19.5%, which is clearly a not inconsiderable fraction of the human genome.

Central Nervous System

The disorders of the CNS that were noted to be elevated on bivariate analyses were: anophthalmos/microphthalmos, anencephalus and similar, anophthalmos, craniosynostosis, encephalocele, eye, hydrocephalus, nervous system, neural tube defects, severe microcephaly and spina bifida. This is consistent with published data. The neural tube defects anencephaly, spina bifida and encephalocele were previously noted to be elevated in Canada after cannabis exposure. Hydrocephaly and microcephaly were previously identified in the Hawaiian series [4]. Hydrocephalus and microcephalus were positively identified with prenatal cannabis exposure in USA [55]. Microcephalus was also associated with cannabis use in Colorado [5].

Face

The facial anomalies that were positively associated on bivariate testing with metrics of cannabis exposure were: anotia, choanal atresia, cleft lip ± palate, cleft palate, congenital cataract, congenital glaucoma, ear, face and neck, holoprosencephaly and similar and orofacial clefts.

Cleft lip with and without cleft palate and cleft palate were also associated with cannabis exposure in USA [55]. Holoprosencephaly was also likely cannabis related in USA [55]. Choanal atresia was positively associated with cannabis use in Hawaii and USA [4, 55].

Gastrointestinal Disorders

The gastrointestinal disorders that were identified in bivariate analyses as being cannabis related were: annular pancreas, anorectal stenosis or atresia, bile duct atresia, digestive system disorders, duodenal stenosis or atresia, Hirschsprungs disease, oesophageal stenosis or atresia and small intestine stenosis or atresia. Of this list, gastrointestinal disorders, biliary atresia and small bowel stenosis or atresia were strongly identified in the USA datasets [55]. Pyloric stenosis and atresia and large bowel and anorectal stenosis and atresia were positively related to antenatal cannabis exposure in Hawaii [4] but were not tracked in Europe. Large bowel disorders and Hirschsprung’s disease were strongly cannabis related in USA [55]. Gastrointestinal disorders, including small and large bowel stenoses and atresias, were also noted to be positively associated with cannabis exposure in Australia [8].

Limb Disorders

The limb disorders that were positively associated with cannabis exposure on bivariate analysis were: club foot, hip dysplasia, limb, limb reductions, polydactyly and syndactyly. Interestingly, syndactyly, polydactyly and reduction deformity of the upper limbs were first identified in the Hawaiian series [4]. Reduction of the lower limbs was a positive finding in the USA series [55]. Musculoskeletal disorders were also noted to be elevated in Colorado in association with increased cannabis exposure.

Interestingly, recent reports from France and Germany [112-115], where cannabis use is rising, relate to ‘outbreaks’ of babies borne without limbs but no such reports have been issued from nearby Switzerland where cannabis is not allowed to enter the food chain.

An increased incidence of limb defects is reminiscent of the notorious teratogenic agent thalidomide to which we directly owe the modern system of pharmacological regulation and safeguards [116]. Concerningly 21 of the 31 CAs described following thalidomide exposure are documented in the present series in association with cannabis exposure [116-122]. Similarly, cannabis shares most (12/13) of the same mechanisms of molecular action as thalidomide .

Body Wall

One of the CAs for which the strongest evidence exists for cannabis association is gastroschisis [4, 123–128]. Moreover, gastroschisis was shown to be 3-fold elevated after multivariable adjustment for all likely confounders in a careful Canadian study [125]. It was therefore of considerable interest to investigate the relationships in the European dataset.

The body wall anomalies with which cannabis was found to be significantly associated were: abdominal wall defect, diaphragmatic hernia, gastroschisis and omphalocele. Gastroschisis has been previously linked with cannabis use also in Hawaii [4], Colorado [5], Australia [8] and Canada [7, 125, 129–131]. Omphalocele was also associated with cannabis use both in Australia [8] and in preclinical series in animals [85]. Diaphragmatic hernia and gastroschisis have previously been linked with cannabis use in USA [87].

Uronephrological Disorders

Uronephrolgical disorders that were found to be elevated on bivariate analyses included: bilateral renal agenesis, bladder extrophy/epispadias, hydronephrosis, hypospadias, multicystic renal disorders, posterior urethral valve and urinary disorders. Obstructive genitourinary disorders was positively associated with prenatal cannabis use in Hawaii [4] and USA [55] and congenital posterior urethral valve was also positively associated in USA [55]. Hence, the European series significantly extends work on this body system.

Genital Disorders

The genital disorder that was found to be elevated after increased population cannabis exposure in bivariate analyses was genital disorder. Indeterminate sex was also noted to be elevated in several multivariable IPW panel analyses. Genitourinary anomalies were also noted to be elevated in Colorado [5], and epispadias was found to be elevated in USA data [55] and hypospadias in Australia [8].

Respiratory

The respiratory defects that were identified on bivariate analysis to be linked with cannabis exposure were cystic lung and respiratory anomalies. Respiratory anomalies were previously identified as a cannabis-associated defect both in Colorado [5] and in USA [55]. Cystic anomalies of the lung have been described as occurring in association with cannabis exposure in later life [132-134].

General

The CAs considered under the ‘general’ category and found to be elevated at bivariate analyses were: all anomalies, amniotic band, conjoined twins, foetal alcohol syndrome, lateral anomalies, maternal infections resulting in malformations, situs inversus, skeletal dysplasias, teratogenic syndromes, valproate syndrome and VATER/VACTERL.

The presence of foetal alcohol syndrome on this list is noteworthy for several reasons. There is evidence of co-abuse of substances of addiction. Many studies document the gateway effect of both cannabis and alcohol leading on to further drug use. Interestingly the foetal alcohol syndrome has been shown to be mediated epigenetically partly via the cannabinoid receptor type 1 (CB1R) [135-138].

Of particular interest is the rare syndrome VATER/VACTERL. Whilst for a long period the reason that these multiple syndromes were observed together was unclear, this has recently been attributed to the inhibition of the key human morphogen sonic hedgehog [57]. Sonic hedgehog plays a pivotal and key and irreplaceable role in human morphogenesis at many points in most organ systems. The positive association of cannabis exposure with this syndrome is important in that it demonstrates at once both that the general view that cannabis is not associated with defined CAs is incorrect and also that cannabis can affect the development of many body systems

Cannabinoids have also been shown to inhibit other key body morphogens, including retinoic acid signalling [139-141], wnt signalling [142-147], the hippo pathway [31], notch signalling [148-152], fibroblast growth factor [153, 154], including transactivation of the FGF1R by CB1R [155], and bone morphogenetic proteins [156-158].

Epidemiological Causal Assignment

Inverse probability weighting effectively removes the non-comparability of groups of interest and creates a situation where groups can be properly compared. It is now well established that inverse probability weighting has the effect of ‘pseudo-randomizing’ an observational dataset, which transforms its findings from merely observational-level associations to replicate a randomized trial in a real-world setting. The reanalysis of several observational trials using this procedure has been done and was shown to reliably replicate the results of subsequent randomized clinical trials [159]. Therefore, the use of this procedure strengthens the present findings.

The other major weakness of observational studies is that uncontrolled extraneous factor(s) might account for the apparently causal nature of an association. The E-value quantifies the degree of association required of some unmeasured confounding covariate with both the exposure of concern and the outcome of interest to obviate the described effect. The E-value usually associated with causality is 1.25 [81], and an E-value of 9, such as is represented by the tobacco–lung cancer relationship, represents a large effect [160]. Many of the mEV reported herein are much larger than this, and range up to 2.46 × 1039. The E-values in Table 8 range at six years lag from 18.45 to 3.61 × 1033 so that they are clearly in a range far beyond what could reasonably be ascribed to extraneous confounding.

Hence, it can be properly stated that the reported results in this study fulfil the epidemiological criteria of causality. Naturally, this comment is not intended to imply that further experimental, laboratory, genomic and epigenomic studies are not required for, indeed, our results underscore the importance of such further and critically important mechanistic investigations.

Mechanisms and Pathways

The mechanisms by which cannabis exposure might induce genotoxicity or epigenotoxicity are numerous and diverse. Plethoric effects on many reproductive organs, sperm stem cells, spermatids, oocytes, cell division, chromosomes, DNA bases, genes, histones and the epigenome have all been described [21-38]. Partly because there are eight different receptors for cannabinoids described [161], the subject is complex and has been reviewed in considerable detail elsewhere [18, 31, 35, 36, 53, 90, 162–167]. It is intended here to make only a few remarks by way of introduction and overview.

Epigenomic Mechanisms

Much elegant work has been undertaken in recent years documenting that quite widespread changes of DNA methylation are induced by cannabis exposure [37] and that these can be passed to subsequent generations of mice [33-37] and can affect the epigenetic regulation of key medium spiny neurons in the nucleus accumbens of the brain, which is a key appetitive driver centre in subsequent generations, which, in turn, affects the proclivity to develop major disorders, such as (opioid and cocaine) drug dependency syndromes [33-37]. It has also been shown that the epigenomic changes in rat and human sperm are equivalent at the pathway level and that such changes in rats were transmissible to subsequent generations [31]. Indeed, it was recently shown that the cessation of cannabis exposure for 17 weeks allows many of the epigenomic changes seen in rats and humans to reverse [32].

Relative Histone Deficiency

It has been shown both in classical focussed investigations of histone synthesis [28–30, 41, 168, 169] and of proteomic screens [40] that histone synthesis in cannabis-exposed tissues is greatly reduced, sometimes by as much as 50%. This has far-reaching implications for epigenomic dysregulation but remains relatively unexplored by modern methods. This area of induced epimutagensis is clearly a very fertile area for further research in much the same way that classical mutagenesis studies have proved to be invaluable as a resource for investigation of gene and protein structure and function in classical genomics.

Epigenomic Implications of Mitochondriopathy

The mitochondriopathic effects of cannabinoids have been richly investigated in considerable detail but remain largely overlooked in the current broader investigative environment.

In their inner and outer membrane complexes and intermembrane spaces, mitochondria carry the full complement of cannabinoid signal reception and transduction machinery as is found in the plasmalemma. They are therefore competent to receive and transduced signals from lipophilic cannabinoids, which will clearly move freely across lipid bilayer membranes.

Whilst many of the proteins found within mitochondria are encoded by the 16 KB of mitochondrial DNA, some are not. Hence, there must be a coordination in the expression of the mitochondrial and nuclear genes to allow normal mitochondrial function. Similar remarks apply to the supply of nuclear energy and epigenomic substrates. This mitonuclear communication occurs through various small-molecule shuttles, including nicotinamide mononucleotide and glutamate/aspartate [52]. This indirect system is known as mitonuclear balance [52].

Mitochondrial inhibition from cannabinoids has been shown to occur through many mechanisms, including uncoupling oxidative phosphorylation by way of uncoupling protein 2 and increasing transmembrane proton leak, and a reduction of many of the key cytochromes of the electron transport chain, including the F1 ATPase itself, which finally harnesses the proton motive force to drive ATP synthesis [40].

Since mitochondria generate the bulk of the small molecules that are used as epigenomic substrates and co-factors (such as methyl, acetyl, phosphate, phosphoribosyl and many other groups) and the bulk of cellular ATP for the largely energy-dependent genomic and epigenomic reactions, it follows that inhibition of intermediary mitochondrial metabolism necessarily has a profound impact on epigenomic expression by limiting the supply of both substrates and energy. Moreover, mitochondrial inhibition will also greatly perturb the indirect pathways of mitonuclear balance.

Indeed, membranous structural continuity has recently been demonstrated between several subcellular organelles and the nucleus [170].

These areas would appear to provide fertile areas for further detailed mechanistic studies.

Cannabinoid DNA Methylome—Ageing DNA Methylome

Derangement of the epigenome is one of the major potential concerns to emerge from this study. Given cannabis exposure has previously been linked with accelerated organismal ageing clinically [95] and since both ageing and cannabinoid pathophysiology are mediated importantly through the epigenome, one fertile field of further investigation would be the intersection and interaction of these two DNA methylomes, including DNA methylation ages [171-173]. Such studies would powerfully inform the multisystem derangements described herein.

Cannabis-induced changes to histone phosphorylation and acetylation states have also previously been documented [41]. The serious changes induced in sperm motility by altered tubulin glycosylation have also recently been documented and were alluded to above as part of the important ‘tubulin code’ [42]. Hence, important work on comparative cannabinoid perturbations of proteomes, phosphoproteomes and post-translational proteomic changes, including altered glycan physiology (as a major post-translational modification) [174, 175] and their altered functional metabolomic and epigenomic implications, are similarly large gaps in the literature, which remain to be addressed.

Non-Coding RNA Expression

The areas of cannabinoid impacts on signalling RNAs of various classes, such as short interfering RNAs, long non-coding RNAs and promoter, enhancer and superenhancer expression and activity are relatively unexplored to our knowledge and also constitute very important areas for further research, especially in view of the organ-specific CAs and tumourigenesis described and alluded to above. Cannabinoid-induced changes to transfer RNAs and gonadal piwi interacting RNAs have been described [176].

Further Research

The present findings strongly indicate the need for further research in both epidemiological and mechanistic fields. Geospatial statistical studies could be extended by the use of smaller geographical units and the use of formal spatiotemporal geostatistical methods. Mechanistic studies are also indicated into the myriad of cannabinoid effects on reproductive toxicity, its effects on mitosis and meiosis, chromosomal and genomic toxicity and the many complex and interacting epigenetic mechanisms, which are likely to be interactively in play.

Generalizability

We feel that these data are widely generalizable for several reasons, including the internal consistency of results within this study, their external consistency with other studies of large datasets published internationally [4–9, 17, 18, 55, 86, 87], the demonstration that the criteria of epidemiological causality are fulfilled (by inverse probability weighting and E-value studies), their robust mechanistic support from the cellular and molecular literature in laboratory investigations [24–26, 28–30, 168, 169, 177–180] and strong support from the preclinical prenatal cannabinoid exposure literature involving several model animal species [83, 85].

The European and USA datasets together comprise the majority of the global datasets on these issues. The fact that similar results have been found in both datasets provides strong external validation a posteriori to both sets of studies of the veracity of their findings.

Strengths and Limitations

This study has a number of strengths and limitations. Strengths include the use of a full panel of substance-related and economic covariates, the use of the formal quantitative techniques of causal inference, the use of multipanelled graphs to visualize whole datasets at one glance and the use of space-time panel regression. We also used formal random forest regression techniques to evaluate variable importance formally. The study has been done on a background of a current understanding of both the cannabis-related epidemiological and mechanistic literature. Its shortcomings relate mainly to the fact that in common with many large epidemiological studies individual participant level data on cannabinoid and substance exposure is not available. Linear interpolation was used to complete some drug data fields.

Conclusion

Using various metrics of cannabis exposure, this study provides compelling evidence of cannabis exposure in the aetiology of a broad range of CAs observed in Europe. Particularly notable on this list were all anomalies and anomalies of the cardiovascular, central nervous, gastrointestinal, chromosomal and genetic, genital, uronephrological, limb and body wall systems and the multisystem disorder VACTERL. These results are consistent with and substantially extend recently documented results from several other jurisdictions [4–9, 17, 18, 55, 86, 87]. They are also consistent with in vitro and preclinical studies from the 1960s to the present time [24–26, 28–30, 168, 169, 177–180]. Data specifically implicate high-intensity daily use and amply confirm concerns raised in relation to the triple epidemiological convergence of rising cannabis use prevalence, rising cannabis use intensity and rising THC concentrations in cannabis herb and resins as a particularly potent driver of cannabis-related disorders [15, 16]. Results are also in accord with recent reports of the epidemiological implication of cannabis with cancers of many types [17, 18, 20, 89] and, since epigenotoxicity has been definitively linked with cellular ageing [19], carry far-reaching implications for programmes and policies, which would lead to widespread genotoxic/epigenotoxic damage across the community under the erroneous assumption that the known epigenotoxic/genotoxic effects of cannabinoids are of minimal clinical significance. The present results amply document the fallacy of this assumption and its severe and ominous implications for population health policy. Further laboratory research, particularly relating to the genotoxic, epigenotoxic, metabolomic-mitochondriopathic-epigenomic and chromosomal toxicity of diverse cannabinoids together with high-resolution spatiotemporal epidemiological studies are strongly indicated.

Click here for additional data file.

Abbreviations

Acronym	Expanded meaning
AME	Average marginal effect
AFE	Attributable fraction in the exposed
ASD	Atrial septal defect
AVSD	Atrioventricular septal defect
CA	Congenital anomaly
CAR	Congenital anomaly rate
CB1R	Cannabinoid receptor type 1
EMCDDA	European Monitoring Centre for Drugs and Drug Addiction
E-Value	Expected value
IPW	Inverse probability weighted
LMC_Herb_Daily	Last month cannabis use × herbal THC content × daily use interpolated
LMC_Resin_Daily	Last month cannabis use × cannabis resin THC content × daily use interpolated
mEV	Minimum E-value
PAR	Population attributable risk
PDA	Patent ductus arteriosus
RR	Relative rate
THC	Δ9-tetrahydrocannabinol
VATER/VACTERL	Vertebral, anorectal, cardiac, tracheo-esophageal fistula ± oesophageal atresia, renal anomalies and limb abnormalities syndrome
VSD	Ventricular septal defect
WHO	World Health Organization