{'fraction_non_alphanumeric': 0.030518545192653942, 'fraction_numerical': 0.019481454807346057, 'mean_word_length': 5.285051488061559, 'pattern_counts': {'":': 0, '<': 2, '<?xml version=': 0, '>': 1, 'https://': 0, 'lorem ipsum': 0, 'www.': 0, 'xml': 0}, 'pii_count': 0, 'substrings_counts': 0, 'word_list_counts': {'cursed_substrings.json': 5, 'profanity_word_list.json': 0, 'sexual_word_list.json': 0, 'zh_pornsignals.json': 0}}
Cannabinoid and substance relationships of European congenital anomaly patterns: a space-time panel regression and causal inferential study

Abstract 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][2][3], recent independent reports from Hawaii [4], Colorado [5], Canada [6,7], Australia [8] and USA [9][10][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.
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 cooccurring 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.

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 useintensity-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][67][68][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.

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

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.
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 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).
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.
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.

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,  Figure 6: Volcano plot of significance (negative log (P-value)) against log (mEV) for daily cannabis use interpolated 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  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.

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 × 10 9 .
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 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.
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.

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 × 10 16 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            Data are summarized by organ system in eTable 58 and presented graphically in eFigs 44 and 45. When considering E-and Pvalues, 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). 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 cannabisrelated 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 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 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.

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 (

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  Table 2. It is of interest that the table is headed by VACTERL     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 .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.

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. Log (Median Percent AME) 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

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][84][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][32][33][34][35][36][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][92][93][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.
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][5][6][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].
If one aggregates all of these chromosomal disorders together, noting that chromosomes 12, 13, 18, 9, 21 and

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,  [55]. Musculoskeletal disorders were also noted to be elevated in Colorado in association with increased cannabis exposure. Median Percent Increment at AME 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

Body Wall
One of the CAs for which the strongest evidence exists for cannabis association is gastroschisis [4,[123][124][125][126][127][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][130][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][133][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][136][137][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][140][141], wnt signalling [142][143][144][145][146][147], the hippo pathway [31], notch signalling [148][149][150][151][152], fibroblast growth factor [153,154], including transactivation of the FGF1R by CB1R [155], and bone morphogenetic proteins [156][157][158].

Epidemiological Causal Assignment
Inverse probability weighting effectively removes the noncomparability 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 'pseudorandomizing' 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 Log (Median Percent AME) 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 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 × 10 39 . The E-values in Table 8 range at six years lag from 18.45 to 3.61 × 10 33 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.

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][34][35][36][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][34][35][36][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 energydependent 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][172][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 posttranslational 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 highresolution spatiotemporal epidemiological studies are strongly indicated.