key: cord-1015922-j2w1tsf8
authors: Reece, Albert Stuart; Hulse, Gary Kenneth
title: Cannabinoid and substance relationships of European congenital anomaly patterns: a space-time panel regression and causal inferential study
date: 2022-02-03
journal: Environ Epigenet
DOI: 10.1093/eep/dvab015
sha: 201677b5754dc04b9358a6f5bdd90f9b62cf1693
doc_id: 1015922
cord_uid: j2w1tsf8

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 × 10(39)) 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.

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.

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] [22] [23] [24] [25] [26] [27] [28] [29] [30] [31] [32] [33] [34] [35] [36] [37] [38] . Several studies have also shown cannabis to have a major effect perturbing DNA methylation [31] [32] [33] [34] [35] [36] [37] , with these changes shown to be inheritable to subsequent generations [31] [32] [33] [34] [35] [36] [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] [44] [45] [46] [47] [48] [49] [50] [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 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 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.

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

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.

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

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.

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 

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.

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.

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.

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.

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.

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.

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 .

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 ( 

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.

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

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.

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 × 10 12 ; eTable 46), four (P = 4.20 × 10 −6 , mEV = 1.11 × 10 4 ; 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 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] .

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] [98] [99] [100] [101] [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] [103] [104] [105] [106] [107] [108] [109] [110] [111] .

If one aggregates all of these chromosomal disorders together, noting that chromosomes 12, 13, 18, 9, 21 and 

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

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

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

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

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

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.

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

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

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

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.

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] [22] [23] [24] [25] [26] [27] [28] [29] [30] [31] [32] [33] [34] [35] [36] [37] [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] [163] [164] [165] [166] [167] . It is intended here to make only a few remarks by way of introduction and overview.

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

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.

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.

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.

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

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.

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.

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.

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.

The risks of marijuana use during pregnancy

Marijuana use during stages of pregnancy in the United States

Prevalence of substance use disorders by time since first substance use among young people in the US

Risk of selected birth defects with prenatal illicit drug use

Cannabis teratology explains current patterns of coloradan congenital defects: the contribution of increased cannabinoid exposure to rising teratological trends

Cannabis consumption patterns explain the east-west gradient in Canadian neural tube defect incidence: an ecological study

Canadian cannabis consumption and patterns of congenital anomalies: an ecological geospatial analysis

Broad Spectrum epidemiological contribution of cannabis and other substances to the teratological profile of northern New South Wales: geospatial and causal inference analysis

Cannabis and pregnancy don't mix

Cannabis in pregnancy -rejoinder, exposition and cautionary tales

Epidemiological overview of multidimensional chromosomal and genome toxicity of cannabis exposure in congenital anomalies and cancer development

Rates of prenatal cannabis use among pregnant women before and during the COVID-19 pandemic

Trends in selfreported and biochemically tested marijuana use among pregnant females in California from

Recommendations from cannabis dispensaries about first-trimester cannabis use

Quadruple convergence -rising cannabis prevalence, intensity, concentration and use disorder treatment

Public health monitoring of cannabis use in Europe: prevalence of use, cannabis potency, and treatment rates

Cannabinoid exposure as a major driver of pediatric acute lymphoid Leukaemia rates across the USA: combined geospatial, multiple imputation and causal inference study

Epidemiological overview of multidimensional chromosomal and genome toxicity of cannabis exposure in congenital anomalies and cancer development

Reprogramming to recover youthful epigenetic information and restore vision

A geospatiotemporal and causal inference epidemiological exploration of substance and cannabinoid exposure as drivers of rising US pediatric cancer rates

Cannabinoid receptor 1 influences chromatin remodeling in mouse spermatids by affecting content of transition protein 2 mRNA and histone displacement

Human sperm express cannabinoid receptor Cb1, the activation of which inhibits motility, acrosome reaction, and mitochondrial function

The cannabinoid system and male reproductive functions

Morphological and cytochemical effects of marijuana cigarette smoke on epithelioid cells of lung explants from mice

Effects of marijuana and tobacco smoke on human lung physiology

Chromosome breakage in users of marihuana

Low doses of widely consumed cannabinoids (cannabidiol and cannabidivarin) cause DNA damage and chromosomal aberrations in human-derived cells

Genetic and Perinatal Effects of Abused Substances

Effects of cannabinoids on spermatogensis in mice

Influence of cannabinoids on somatic cells in vivo

Cannabinoid exposure and altered DNA methylation in rat and human sperm

Refraining from use diminishes cannabis-associated epigenetic changes in human sperm

Maternal cannabis use alters ventral striatal dopamine D2 gene regulation in the offspring

Parental THC exposure leads to compulsive heroin-seeking and altered striatal synaptic plasticity in the subsequent generation

Epigenetic effects of cannabis exposure

High times for cannabis: Epigenetic imprint and its legacy on brain and behavior

Genome-wide DNA methylation profiling reveals epigenetic changes in the rat nucleus accumbens associated with cross-generational Effects of adolescent THC exposure

Effects of cannabis and natural cannabinoids on chromosomes and ova

Maternal cannabis use alters ventral striatal dopamine D2 gene regulation in the offspring

Genes and pathways co-associated with the exposure to multiple drugs of abuse, including alcohol, amphetamine/ methamphetamine,cocaine, marijuana, morphine, and/or nicotine: a review of proteomics analyses

Influence of psychoactive and nonpsychoactive cannabinoids on chromatin structure and function in human cells

Tubulin glycylation controls axonemal dynein activity, flagellar beat, and male fertility

Regulatory effects of cannabidiol on mitochondrial functions: a review

Cannabidiol targets mitochondria to regulate intracellular Ca2+ levels

Dose-dependent cannabidiol-induced elevation of intracellular calcium and apoptosis in human articular chondrocytes

A cannabinoid link between mitochondria and memory

The influence of delta9-tetrahydrocannabinol, cannabinol and cannabidiol on tissue oxygen consumption

Hypothalamic POMC neurons promote cannabinoid-induced feeding

S)Pot on mitochondria: cannabinoids disrupt cellular respiration to limit neuronal activity

Inhaled marijuana smoke disrupts mitochondrial energetics in pulmonary epithelial cells in vivo

Delta 9-tetrahydrocannabinol disrupts mitochondrial function and cell energetics

NAD(+) metabolism and the control of energy homeostasis: a balancing act between mitochondria and the nucleus

Chromothripsis and epigenomics complete causality criteria for cannabis-and addiction-connected carcinogenicity, congenital toxicity and heritable genotoxicity

Cannabinoid exposure as a major driver of pediatric acute lymphoid Leukaemia rates across the USA: combined geospatial, multiple imputation and causal inference study

Cannabis in pregnancy -rejoinder, exposition and cautionary tales

VACTERL association: information for families Great Ormond St

Cannabinoids Exacerbate Alcohol Teratogenesis by a CB1-Hedgehog Interaction

European Commission

The World Bank: crude data: birth rate 2021

European Monitoring Centre for Drugs and Drug Addiction (EMCDDA): Statistical Bulletin 2021 -prevalence of drug use Lisbon, Portugal: European Monitoring Centre for Drugs and Drug Addiction (EMCDDA)

The World Bank: Crude Data: Adjusted net national income per capita (current US$): The World Bank

Welcome to the Tidyverse

Margins: marginal effects for model objects

Web site and r package for computing E-values

Sensitivity analysis in observational research: introducing the E-value

Commentary: developing bestpractice guidelines for the reporting of E-values

Tools for the analysis of epidemiological data

Ranger: a fast implementation of random forests for high dimensional data in C++ and R

Variable importance plots-an introduction to the vip package

gplots: Various R programming tools for plotting data. R package version 311

Broom: convert statistical objects into tidy tibbles: R package version

mixed: tidying methods for mixed models

Teoria statistica delle classi e calcolo delle probabilità Editor

Simple Sequentially A. Rejective Multiple Test Procedure

Controlling the false discovery rate: a practical and powerful approach to multiple testing

CRAN: Central ' 'R' 'Archive Network

IPW: an R package for inverse probabilty weighting

Technical considerations in the use of the E-value

Human Embryology and Developmental Biology 6th edn

Teratogenicity of marihuana extract as influenced by plant origin and seasonal variation

Effect of marihuana extract on fetal hamsters and rabbits

Cannabis and Health. 1 1st edn

Noninherited risk factors and congenital cardiovascular defects: current knowledge: a scientific statement from the American Heart Association Council on Cardiovascular Disease in the Young: endorsed by the American Academy of Pediatrics

Using Bayesian models to assess the effects of under-reporting of cannabis use on the association with birth defects, national birth defects prevention study

Maternal periconceptional illicit drug use and the risk of congenital malformations

Causal inference multiple imputation investigation of the impact of cannabinoids and other substances on ethnic differentials in US testicular cancer incidence

Impacts of cannabinoid epigenetics on human development: reflections on Murphy et. al. 'cannabinoid exposure and altered DNA methylation in rat and human sperm

Epidemiological associations of various substances and multiple cannabinoids with autism in USA

Effect of cannabis legalization on US autism incidence and medium term projections

Gastroschisis and autism-dual canaries in the Californian coalmine

Co-occurrence across time and space of drug-and cannabinoid-exposure and adverse mental health outcomes in the National Survey of Drug Use and Health: combined geotemporospatial and causal inference analysis

Cannabis exposure as an interactive cardiovascular risk factor and accelerant of organismal ageing: a longitudinal study

Contemporary epidemiology of rising atrial septal defect trends across USA 1991-2016: a combined ecological geospatiotemporal and causal inferential study

Cannabis use and incidence of testicular cancer: a 42-year follow-up of Swedish men between

Association of marijuana use and the incidence of testicular germ cell tumors

Cannabis exposure and risk of testicular cancer: a systematic review and meta-analysis

Population-based casecontrol study of recreational drug use and testis cancer risk confirms an association between marijuana use and nonseminoma risk

Incident testicular cancer in relation to using marijuana and smoking tobacco: A systematic review and meta-analysis of epidemiologic studies

Marijuana use and testicular germ cell tumors

Mechanisms of leukemia evolution: lessons from a congenital syndrome

Worldwide comparison of survival from childhood leukaemia for 1995-2009, by subtype, age, and sex (CONCORD-2): a population-based study of individual data for 89 828 children from 198 registries in 53 countries

DNA methylation protects hematopoietic stem cell multipotency from myeloerythroid restriction

Molecular cytogenetics of human germ cell tumours: i(12p) and related chromosomal anomalies

International trends in the incidence of testicular cancer: lessons from 35 years and 41 countries

Testicular cancer-discoveries and updates

A genome-wide association study of testicular germ cell tumor

Testicular germ cell tumours: predisposition genes and the male germ cell niche

Germ cell tumors from a developmental perspective: cells of origin, pathogenesis, and molecular biology (emerging patterns)

Babies born with deformed hands spark investigation in Germany. CNN News

France to investigate cause of upper limb defects in babies. The Guardian [Internet], 3 November

Scientists are baffled by spatter of babies born without hands or arms in France, as investigation fails to discover a cause. Daily Mail

Baby arm defects prompt nationwide investigation in France

The rise, fall and subsequent triumph of thalidomide: lessons learned in drug development

Drug-induced (Thalidomide) malformations

Thalidomide and congenital malformations

The thalidomide saga

Thalidomide induces limb defects by preventing angiogenic outgrowth during early limb formation

Thalidomide-induced limb defects: resolving a 50-year-old puzzle

Thalidomide-induced teratogenesis: history and mechanisms

A case-control study of maternal periconceptual and pregnancy recreational drug use and fetal malformation using hair analysis

Recreational drug use: a major risk factor for gastroschisis?

Maternal risk factors for gastroschisis in Canada

A population-based study of gastroschisis: demographic, pregnancy, and lifestyle risk factors

Maternal periconceptional illicit drug use and the risk of congenital malformations

Association of vasoconstrictive exposures with risks of gastroschisis and small intestinal atresia

Canadian Pediatric Surgery N. Outcomes of early versus late intestinal operations in patients with gastroschisis and intestinal atresia: results from a prospective national database

Spatial variability of gastroschisis in Canada

A perinatal health surveillance report

Cannabis as a cause of giant cystic lung disease

Giant cystic lung disease with mediastinal compression in a short-term heavy cannabis smoker

Response to aldington: cannabis and lung cancer

Maternal alcohol exposure during mid-pregnancy dilates fetal cerebral arteries via endocannabinoid receptors

CB1R-mediated activation of caspase-3 causes epigenetic and neurobehavioral abnormalities in postnatal ethanol-exposed mice

Ethanol exposure induces neonatal neurodegeneration by enhancing CB1R Exon1 histone H4K8 acetylation and up-regulating CB1R function causing neurobehavioral abnormalities in adult mice

Postnatal ethanol exposure alters levels of 2-arachidonylglycerol-metabolizing enzymes and pharmacological inhibition of monoacylglycerol lipase does not cause neurodegeneration in neonatal mice

Lipid abundance in zebrafish embryos is regulated by complementary actions of the endocannabinoid system and retinoic acid pathway

Genetic susceptibility to posttraumatic stress disorder: analyses of the oxytocin receptor, retinoic acid receptor-related orphan receptor A and cannabinoid receptor 1 genes

Retinoic acids and hepatic stellate cells in liver disease

Effects of cannabidiol interactions with Wnt/beta-catenin pathway and PPARgamma on oxidative stress and neuroinflammation in Alzheimer's disease

Identification of synergistic interaction between cannabis-derived compounds for cytotoxic activity in colorectal cancer cell lines and colon polyps that induces apoptosis-related cell death and distinct gene expression

Identifying novel members of the Wntless interactome through genetic and candidate gene approaches

miR-23b-3p and miR-130a-5p affect cell growth, migration and invasion by targeting CB1R via the Wnt/beta-catenin signaling pathway in gastric carcinoma

Non-canonical Wnt signaling through Ryk regulates the generation of somatostatinand parvalbumin-expressing cortical interneurons

Analyzing the role of cannabinoids as modulators of Wnt/beta-catenin signaling pathway for their use in the management of neuropathic pain

Opposing actions of endocannabinoids on cholangiocarcinoma growth is via the differential activation of Notch signaling

The endocannabinoid, anandamide, augments Notch-1 signaling in cultured cortical neurons exposed to amyloid-beta and in the cortex of aged rats

Activation of type 1 cannabinoid receptor (CB1R) promotes neurogenesis in murine subventricular zone cell cultures

Ubiquitination-dependent CARM1 degradation facilitates Notch1-mediated podocyte apoptosis in diabetic nephropathy

Potentiation of the antitumor activity of adriamycin against osteosarcoma by cannabinoid WIN-55,212-2

The CB1 cannabinoid receptor mediates excitotoxicity-induced neural progenitor proliferation and neurogenesis

The FGF receptor uses the endocannabinoid signaling system to couple to an axonal growth response

Cannabinoid 1 receptordependent transactivation of fibroblast growth factor receptor 1 emanates from lipid rafts and amplifies extracellular signalregulated kinase 1/2 activation in embryonic cortical neurons

Manipulating molecular switches in brown adipocytes and their precursors: a therapeutic potential

Brown fat biology and thermogenesis

Orexin receptors: multi-functional therapeutic targets for sleeping disorders, eating disorders, drug addiction, cancers and other physiological disorders

Methods of public health research -strengthening causal inference from observational data

The Book of Why. The New Science of Cause and Effect

Microglial activation and cannabis exposure

Chronic toxicology of cannabis

Epigenetic regulation of immunological alterations following prenatal exposure to marijuana cannabinoids and its long term consequences in offspring

Adolescent exposure to Delta(9)-tetrahydrocannabinol alters the transcriptional trajectory and dendritic architecture of prefrontal pyramidal neurons

Prenatal cannabis exposure -The ' 'first hit' ' to the endocannabinoid system

DNA methylation and substance-use risk: a prospective, genome-wide study spanning gestation to adolescence

Cannabinoids and gene expression during brain development

Nonmutagenic action of cannabinoids in vitro

Genetic effects of marijuana

Alterations in inter-organelle crosstalk and Ca(2+) signaling through mitochondria during proteotoxic stresses

Epigenetic predictor of age

DNA methylation age of human tissues and cell types

Age-associated epigenetic drift: implications, and a case of epigenetic thrift?

Profile of immunoglobulin G N-glycome in COVID-19 patients: a case-control study

Pathways from epigenomics and glycobiology towards novel biomarkers of addiction and its radical cure

The epigenetics of the endocannabinoid system

Cannabis physiopathology epidemiology detection

Keep off the Grass

Effects of cannabinoids on macromolecular synthesis and replication of cultured lymphocytes

Inhibition of cellular mediated immunity in marihuana smokers

All authors had full access to all the data in the study and take responsibility for the integrity of the data and the accuracy of the data analysis. Prof. Mark Stevenson at the University of Melbourne is very gratefully thanked for specially preparing a new version of epiR to handle the very large integers involved in this analysis. Prof. Maya Mathur of Stanford Medical School is thanked for her advice, assistance and guidance with E-values.

All data generated or analysed during this study are included in this published article and its supplementary information files. Data along with the relevant R code have been made publicly available on the Mendeley Database Repository and can be accessed from this URL doi: 10.17632/vd6mt5r5jm.1.

Interestingly, recent reports from France and Germany [112] [113] [114] [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] [117] [118] [119] [120] [121] [122] . Similarly, cannabis shares most (12/13) of the same mechanisms of molecular action as thalidomide .

Supplementary data is available at EnvEpig online.

No funding was provided for this study. No funding organization played any role in the design and conduct of the study; collection, management, analysis and interpretation of the data; preparation, review, or approval of the manuscript; and decision to submit the manuscript for publication.

The authors declare that they have no competing interests. 

Not applicable.

A.S.R. assembled the data, designed and conducted the analyses and wrote the first manuscript draft. G.K.H. provided technical and logistic support, co-wrote the paper, assisted with gaining ethical approval and provided advice on manuscript preparation and general guidance to study conduct. A.S.R. had the idea for the article, performed the literature search, wrote the first draft and is the guarantor for the article.