European Epidemiological Patterns of Cannabis- and Substance-Related Congenital Neurological Anomalies: Geospatiotemporal and Causal Inferential Study

Introduction. Of the many congenital anomalies (CAs) recently linked with community cannabis exposure, arguably the most concerning are neurological CAs (NCAs). We therefore conducted a detailed study of this in fourteen European nations. Methods. Congenital anomaly data were from Eurocat. Drug exposure data were from European Monitoring Centre for Drugs and Drug Addiction. Income from World bank. Results. The Netherlands, Spain, France and Bulgaria reported increasing rates of many NCAs. The NCA rate (NCAR) was higher in nations with increasing daily cannabis use when compared to those without (p = 0.0204, minimum E-value (mEV) = 1.35). At bivariate analysis, the mEVs of the following NCAs were significantly cannabis related: severe microcephaly 2.14 × 1013, craniosynostosis 5.27 × 1011, nervous system 4.87 × 1011, eye 2.73 × 107, microphthalmos 4.07 × 106, anencephalus 710.37, hydrocephalus 245.64, spina bifida 14.86 and neural tube defects 13.15. At inverse probability, weighted panel regression terms including cannabis were significantly related to the following series of anomalies: nervous system, anencephalus, severe microcephalus, microphthalmos, neural tube defect and spina bifida from p = 5.09 × 10−8, <2.2 × 10−16, <2.2 × 10−16, 4.84 × 10−11, <2.2 × 10−16 and 9.69 × 10−7. At geospatial regression, this same series of anomalies had terms including cannabis significant from p = 0.0027, 1.53 × 10−7, 3.65 × 10−6, 2.13 × 10−8, 0.0002 and 9.76 × 10−12. 88.0% of 50 E-value estimates and 72.0% of mEVs > 9. This analysis therefore demonstrates both close association of cannabis exposure with multiple NCAs across space-time and also fulfills the quantitative criteria of causal inferential analysis. Conclusions. Nine NCARs on bivariate and six NCARs on multivariable regression were cannabis related and fulfilled quantitative epidemiological criteria for causality and are consistent with other series. Particular concerns relate to exponential dose–response effects demonstrated in the laboratory and epidemiological studies. Great caution with community cannabinoid penetration is warranted. Data indicate that cannabis is a significant environmental teratogen and thus imply that cannabinoids should be regulated similarly to the manner in which all other important genotoxins are carefully controlled by communities for their self-sustaining longevity and the protection of generations yet to come.


Introduction
With increasing recent reports describing many congenital anomalies (CAs) attributable to prenatal cannabis exposure, one of the most concerning groups of anomalies is that affecting the brain and neurological development generally [1,2]. With classical studies in the preclinical literature describing hydrocephalus, encephalocele, spina bifida and anencephalus in rats, rabbits and hamsters gestationally exposed to cannabis [3][4][5], and recent confirmatory epidemiological findings relating to encephalocele and hydrocephaly [6] and CDC reports of anencephalus [7], reports of neural tube elevations in Canada [8], increased hydrocephalus in Australia [9] and several increased rates of encephalocele, spina bifida and microphthalmia in the USA [10,11], we were concerned about the profile that would be encountered in Europe where some nations such as Spain, France, Portugal, the Netherlands have reported increasing rates of cannabis use [12,13].
Many of these anomalies affected or were known to relate to brain development. Since congenital causes of being intellectually challenged can be severe, can remain throughout life and may often impose severe degrees of disability, this is an intrinsically severe and important group of congenital anomalies to understand. Moreover, since prenatal cannabis exposure had previously been linked with anencephalus [6,7], autism [14][15][16] and smaller head sizes [17][18][19], we enquired if this apparent continuum of increasing brain damage would be continued in the European dataset, particularly as that database includes an anomaly called "severe microcephalus" which is not tracked elsewhere.
Both cannabinoid genotoxicity [20][21][22][23][24][25][26][27][28][29][30] and cannabinoid-induced disruptions of mitochondrial activity [31][32][33][34][35][36], which support and supply key substrates and energy to genomic and epigenomic reactions [37], have been extensively documented to demonstrate exponential cannabis dose-response relationships. Moreover, these exponential effects [22][23][24][25][26]30,35] have been confirmed in several epidemiological studies which document a substantial discontinuous step in congenital anomaly rates from the fourth to the highest quintile of cannabis exposure [1,10,11]. Since many parts of Europe are well documented to have seen abrupt rises in community exposure to cannabinoids resulting from increased prevalence of cannabis use, cannabis use daily intensity, and increased ∆9-tetrahydrocannabinol (THC) potency in marketed cannabis products [12,13], it seemed likely that this would transition the community from the low risk genotoxic zone into the level at high risk for major and serious genotoxic outcomes. Indeed, there is some powerful evidence that this may indeed be occurring via contamination of the food chain at least in those parts of northeastern France where cannabis crops are grown and both calves and human babies are being born limbless [38][39][40].
We decided to include some eye anomalies in this group as the eye is formed as the confluence of two sets of tissues: one coming from the face, which is responsible for the formation of the anterior segments of the eye including the lens, and an outgrowth of the forebrain and optic nerve which is responsible for the posterior parts of the eye including the retina and neural elements [41]. For this reason, we treat anterior eye anomalies separately under a companion paper addressing orofacial disorders. However, the decision was made to include whole-of-eye anomalies in this section. The data report both anophthalmos and anophthalmos/microphthalmos. These two disorders behave quite differently analytically. It seems to us that the former condition referred to anophthalmos properly so-called and the latter actually refers to microphthalmos. These disorders have therefore been treated this way in the following analysis.
For all of these reasons, we were keen to study central nervous system anomalies in the European dataset both in terms of the bivariate associations and also in terms of their multivariable modelling within causal inferential and space-time paradigms. The following report sets out this detailed account. The specific hypotheses investigated were as follows: (1) Did previously described bivariate relationships between cannabis exposure and neurological CAs (NCAs) persist after multivariate adjustment; (2) Did these relationships persist when considered formally in a space-time context; and (3) Did these relationships fulfill the formal criteria of quantitative causal inferential analysis.

Data
Data on all available congenital anomaly rates were downloaded for each of 14 nations by each individual year from the European Network of Population-Based Registries for the Epidemiological Surveillance of Congenital Anomalies (EUROCAT) website [42] and analyzed. It is important to note that EUROCAT total congenital anomaly rate includes anomaly rates amongst live births, stillbirths and cases where early termination for anomaly was practiced all combined together so that it represents a total overall picture across all classes of births. The nations selected were chosen on the basis of the availability all or most of their congenital anomaly data across the years 2010-2019. National alcohol (liters of pure alcohol consumed per capita annually) and tobacco (percent daily tobacco use prevalence) use data were downloaded from the World Health Organization [43]. Drug use data for amphetamines cannabis and cocaine were taken from the European Monitoring Centre for Drugs and Drug Addiction (EMCDDA) [44]. Last month cannabis use data were also supplemented by data on the tetrahydrocannabinol (THC) content of cannabis herb and resin provided in recent reports which have been published [13]. Data on daily cannabis use are also available from EMCDDA and was consistent with data collated in recent reports [13]. Median household income data (in USD) were taken from the World Bank [45].

National Assignment
Nations were categorized as being either high and/or rising daily cannabis use or low and/or falling daily cannabis use based on a recent European epidemiological study of cannabis use in Europe (see Supplementary Figure S4 in reference [13]). Thus, Italy, the Netherlands, Norway, France, Germany, Belgium, Croatia, Portugal and Spain were categorized as nations experiencing increasing daily use, while Hungary, Poland, Finland, Bulgaria, and Sweden were nations which were experiencing low or falling levels of daily cannabis use.

Derived Data
Because several metrics of cannabis use, exposure and consumption were available, it was possible to calculate various derived compound metrics. For this reason, last month cannabis use prevalence data were multiplied by the THC content of cannabis herb and resin to derive their products. These metrics were then multiplied by imputed daily cannabis use rates to derive further compound metrics for both cannabis herb and resin.

Data Imputation
Linear interpolation was used to complete missing data. This was particularly used for daily cannabis use. In all, 59 data points on daily cannabis use from EMCDDA were available for these 14 nations in this period. Linear interpolation was used to expand this dataset to 129 datapoints (further details provided in Results section). Swedish data on cannabis resin THC concentration were not available. However, it was noted that the resin/herb THC concentration was virtually constant in nearby Norway at 17.7 so this ratio was multiplied by the Swedish cannabis herb THC concentration data to derive estimates of Swedish cannabis resin THC concentration. Similarly, Polish data for the cannabis resin THC concentration were unavailable. The resin to herb THC concentration ratio of nearby Germany was used to estimate the resin THC content in Poland from the known herb THC concentrations observed in Poland. Since geospatial analytical techniques do not tolerate missing data, the techniques of last observation carried forward or backwards for Croatia in 2018 and 2019 and the Netherlands in 2010 were used to complete missing data. Multiple imputation methods could not be applied to this dataset as it is not possible to apply multiple imputed datasets in panel or spatial multivariable regression techniques.

Statistics
R Studio version 1.4.1717, based on R version 4.1.1 from the Comprehensive R Archive Network and the R Foundation for Statistical Computing [46], was used for data processing. The analysis was conducted in December 2021. Dplyr from the tidyverse was used for data manipulation [47]. Log transformation of datasets was employed as appropriate to improve compliance with normality assumptions based on the results of the Shapiro-Wilks test. ggplot2 from the tidyverse was used to draw graphs. Maps were drawn using ggplot2 together with sf (simple features) [48] and custom color palettes and palettes taken from the viridis and viridisLite packages were used to generate fill schemas [49]. The R package colorplaner was used to draw and color bivariate maps [50]. Illustrations have not been published previously and are original. Linear regression was performed from Base R. The package nlme was used for mixed effects regression [51]. The classical technique for model reduction, namely serial deletion of the least significant term, was used in all multivariable models to yield a final reduced model which has been model presented. Multiple linear models were processed in a single pass using combined techniques from R packages broom and purrr [47,52,53]. The overall effect of covariates in multivariable models may be quantified and this is denoted as the "marginal effect". In this study, the overall marginal effect in multivariate models was calculated using the R package margins [54].

Covariate Selection
The presence of so many different metrics for cannabis exposure and consumption created a major problem for analysis as it was not clear which was the most appropriate metric to employ in any particular analytical scenario. Indiscriminate use of excessive covariates in a multivariable model would make models impossible to analyze by unnecessarily consuming degrees of freedom and thereby restrict ability to assess interactions. This issue was formally addressed by employing random forest regression using the R package ranger [55] with variable importance being formally studied using the R package vip (variable importance plot) [56]. The most highly predictive set of covariates from this process was entered into the regression modelling equations. The Results section presents the tables from this analysis.

Panel and Geospatial Analysis
R package plm [57] was used to conduct panel analysis across both space and time simultaneously using the "twoways" effect. The sparse spatial weights matrix was calculated using R package spdep (spatial dependency) [58] by using the edge and corner "queen" relationships. The spatial panel random effects maximum likelihood (spreml) function from the package spml, which allows detailed modelling and correction of model error structures [59,60] was used for geospatial modelling. Such models may produce four model coefficients of interest and these in turn are useful in determining the most appropriate error structure for the model. These coefficients are rho the spatial coefficient, phi the random error effect, psi the serial correlation effect, and theta the spatial autocorrelation coefficient. In each case, the most appropriate error structure was chosen for each spatial model generally, taking care to preserve the model error specification across closely related models. The appropriate error structure was determined by the backwards methods from the full general model to the most specific model, as has been previously described [61]. Temporal lagging by one or two years, as indicated, was applied to both panel and geospatial models.

Causal Inference
The formal tools of causal inference were used for this analysis. Inverse probability weighting (ipw) is the technique of choice to convert a purely observational study into a pseudo-randomised study and from these models it is entirely appropriate to draw causal inferences [62]. All multivariable panel models presented in the present study were inverse probability weighted. The R package ipw was used for inverse probability weighting. Similarly E-values (expected values) is a powerful form of sensitivity analysis and may be used to quantify the correlation required of some hypothetical unmeasured confounder covariate with both the exposure of concern and the outcome of interest in order to explain away some apparently causal relationship [63][64][65]. It therefore provides a quantitative measure of the robustness of the model to extraneous covariates which have not been accounted for within the parameters which have been selected and measured. E-Values are associated with a corresponding confidence interval and the 95% lower bound of this confidence interval is reported herein as the minimum E-Value confidence interval. E-Value estimates greater than 1.25 have been reported to indicate causality [66] and E-values greater than nine have been described as high [67]. The R package EValue [68] was used for the calculation of E-values. Both E-values and inverse probability weighting are foundational and pivotal techniques used in formal causal inferential methods in order to allow causal relationships to be assessed from pseudorandomized real world observational studies.

Data Availability
Raw datasets including 3800 lines of computation code in R has been made freely available through the Mendeley data repository at the following URLs: 10.17632/vd6mt5r5jm.1 and 10.17632/gr5ntsbp7p.1.

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

Input Data
Supplementary Table S1 sets out the baseline data for this study. It shows the 14 nations which contributed data via the Eurocat database and the 11 congenital anomalies which were studied in this group. As can be seen, there were 1327 datapoints in this dataset with most anomalies having 110 data points each. Drug use data, including several metrics for cannabis exposure including compound indices, are also shown along with median household income.
Eye anomalies were included with disorders of the CNS as the eye is formed embryologically as a confluence of tissues contributed both from the face and the developing forebrain. The facial tissues contribute the anterior segments of the eye and the retina and neural tissues come from the brain. For this reason, congenital anomalies of the front section of the eye are considered in our companion manuscript looking at facial disorders and disorders of the eye as a whole are considered in the present paper.
Daily drug use data were notable incomplete and is shown in Supplementary Table S2. 59 points are listed there. In order to allow these data to be used in analyses, a further 70 datapoints were added by linear interpolation as shown in Supplementary Table S3. The regression lines for tobacco, alcohol and amphetamine in Supplementary Figure S1 are non-descript and not significant. For some anomalies, the regression relationships of cocaine are positive and strongly significant. For four of the anomalies listed, the compound cannabis metric used is clearly positive. It is noted in Supplementary Figure S1 that the slopes of the regression lines for the cannabis metric and cocaine are opposite for the two anomalies Anophthalmia/Microphthalmia and anophthalmos. For this reason, these two designations are believed to be reporting on different conditions, the first one is believed to be a surrogate for microphthalmia and the second relates to Anophthalmia properly so-called. Scatterplots for six anomalies are shown in Figure 2 and for five NCAs in Supplementary Figure S2. The regression lines for tobacco, alcohol and amphetamine in Supplementary Figure  S1 are non-descript and not significant. For some anomalies, the regression relationships of cocaine are positive and strongly significant. For four of the anomalies listed, the compound cannabis metric used is clearly positive. It is noted in Supplementary Figure S1 that the slopes of the regression lines for the cannabis metric and cocaine are opposite for the two anomalies Anophthalmia/Microphthalmia and anophthalmos. For this reason, these two designations are believed to be reporting on different conditions, the first one is believed to be a surrogate for microphthalmia and the second relates to Anophthalmia properly so-called. Scatterplots for six anomalies are shown in Figure 2 and for five NCAs in Supplementary Figure S2.    Figure 3 is a graphical map showing the rate of central nervous system conditions across Europe over the decade 2010-2019. The area of France appears to have constantly relatively high rates. The rate in Germany started high but has declined from there. The rate has increased across the decade in Spain. Figure 4 shows the pattern of severe microcephaly across the decade. For this anomaly, the situation in Spain and France has declined whilst that in Germany and Belgium has improved. Rates in the Netherlands have fluctuated.

Bivariate Relationships
The space-time patterns of anencephaly are as shown in Figure 5. Rates in Spain, France, Germany and Belgium have increased. Bulgaria has often had a high rate. Polish rates have been relatively lower and have declined across the decade. Rates in the low countries have fluctuated across a relatively high range.       The pattern of microphthalmia is shown in Supplementary Figure S3. For the reasons described above, it is felt that these data reflect microphthalmia more than anophthalmia which appear to be different entities. Rates in Poland have been persistently low and have declined. Those elsewhere have fluctuated.
Rates of neural tube defects (Supplementary Figure S4) have been persistently low in Poland but have deteriorated in France, Germany and Spain. They have often been relatively high in Bulgaria.
Eye anomalies (Supplementary Figure S5) were persistently high in Finland when data were available and rates have increased in France. German and Italian rates have declined across this period.
Rates of one of the compound cannabis metrics, last month cannabis use: cannabis resin THC are shown in Supplementary Figure S6. This index is noted to have increased in Spain, the Netherlands, France and Portugal across this period. Figure 6 is a bivariate plot of overall central nervous system anomaly rates compared to the compound cannabis exposure metric last month cannabis use: cannabis resin THC concentration. One reads the graph by noting that the areas which are high for both covariates are shaded purple or pink. Thus, the maps clearly illustrate the emergence of simultaneously high rates of both parameters in France, Spain, Bulgaria and the Netherlands across the decade.
The pattern of severe microcephaly plotted against the same compound cannabis exposure parameter is shown in Figure 7. Again, dually high rates of both parameters are noted to have emerged in France, Spain and the Netherlands across this decade.
When microphthalmos is considered, the areas of France and Spain are noted to have turned shades of purple across this period, which again indicates dually concordant elevated rates (Supplementary Figure S7).
The convergent pattern of neural tube defects across the continent is shown in Supplementary Figure S8 where simultaneously elevated rates are noted to have occurred in France, the Netherlands and Bulgaria. The rate in Germany is noted to have risen but without a rise in this particular metric of cannabis exposure.

Comparing High and Low Cannabis Use Countries
As described in the Methods section, this group of nations may be separated into those where daily cannabis use is increasing and those where it is decreasing or stable. When nations are dichotomized in this way, the appearance for all the anomalies considered together is as shown in Supplementary Figure S9. At linear regression, it is found that the rate in the nations where cannabis use is higher is significantly greater than those where it is not (β-est. = 0.155, t = 2.332, p = 0.0265; model Adj.R.Squ. = 0.0029, F = 4.934, p = 0.0265). When this issue is considered by each anomaly, the appearance shown in Supplementary Figure S10 is seen. In this Figure, it is clear that the rates in the two sets of nations overlap when considered overall, a finding confirmed by mixed effects regression. However, in the severe microcephaly group there is a significant difference between the groups (β-est. = 0.372, t = 2.349, p = 0.0205; model Adj.R.Squ. = 0.0360, F = 5.517, p = 0.0205).
Supplementary Table S4 shows details from 132 linear regression models for each anomaly considered against each of the substance exposure metrics. From these models, those with positive regression coefficients and significant p-values were extracted as shown in Supplementary Table S5. Sixty models were thus selected. They are listed in descending order of their minimum E-values. Importantly, it is noted that daily cannabis use interpolated and severe microcephaly head up this list of bivariate associations. Forty-seven of these 60 (78.3%) substance exposure terms include cannabis compared to only seven (11.7%) for cocaine, two (3.3%) for tobacco and only one (1.7%) for alcohol.

Multivariable Analysis
The next step is to move from the bivariate context to multivariable regression. However, it is unclear which selection of the 13 exposure covariates is most suitable for each CA as the set of independent variables. Random Forrest regression was used in tandem with variable importance plots to decide this issue. The six variable importance tables for the six CAs of primary interest are shown as Supplementary Tables S6-S11. The reasons for choosing this set of CAs are described in the Discussion section.

Multivariable Relationships-Panel Regression
Supplementary Table S12 shows a series of inverse probability weighted panel regression models for nervous system disorders as an additive and an interactive model and a model lagged by two years. The use of inverse probability weighting is important as it allows us to move from a merely observational study into a pseudo-randomized paradigm from which causal inference can properly be made. One notes in this Table that cannabis exposure metrics are positive in each model with high levels of statistical significance (from 5.09 × 10 −8 in the interactive model).
This pattern is repeated across Supplementary Tables S13-S17 which consider, respectively, anencephalus, severe microcephalus, microphthalmos, neural tube defect and spina bifida. These important findings at multivariable regression confirm that in all six cases metrics of cannabis exposure survive model reduction and persist at high levels of statistical significance in final models. This indicates that the strong effects observed at bivariate regression persist after careful adjustment for relevant covariates in additive, interactive and lagged contexts.

Multivariable Relationships-Geotemporospatial Regression
The next issue of relevance was to consider these data in their native space-time context, where important analytical confounding factors such as serial correlation, spatial correlation, random error structures and spatial autocorrelation could be properly and formally accounted for. For this reason, the geospatial links between countries shown in Supplementary Figure S11 were defined, edited and finalized as shown, which then became the basis for the sparse spatial weights matrix entered into the spatial regression equations. Table 1 presents the final geospatial regression additive, interactive and temporally lagged models along with details of the individual regression parameters employed for nervous system disorders. In the additive and interactive models, terms including cannabis exposure have positive regression coefficients and are statistically significant. This signal is not seen in the temporally lagged models. Table 2 presents final spatiotemporal regression results for anencephalus. Terms relating to cannabis exposure are not seen in the additive model but are notable with high levels of statistical significance in the interactive and the models lagged to one and two years.  Table 3 shows final geospatial models for severe microcephalus. In this case, cannabisrelated terms are seen in each of the models presented and their magnitude dominates the regression result. The effect is not removed by temporal lagging.
Final additive, interactive and lagged geospatial models for microphthalmos are presented in Table 4. Terms with positive regression coefficients are high levels of statistical significance seen in all models. When neural tube defects are considered in Table 5, cannabis metrics with positive coefficients appear in the additive and both lagged models.
Review of final models for spina bifida confirms this pattern with terms for the cannabis metric being noted in all the geospatial models studied (Table 6).

Causal Inference-E-Values
From each of the regression models, it is possible to extract applicable E-values. The E-value quantitates the expected magnitude of the correlation of some unknown confounding variable with both the exposure of concern and the outcome of interest to explain away an apparently causal effect. It is therefore an important part of quantitative causal inferential techniques. E-values for panel and geospatial models are presented in Supplementary Tables S18 and S19, respectively. These Tables are then combined and the 50 E-values are shown together in Supplementary Table S20 listed in descending order of minimum E-value. This table is notable for several features. Severe microcephalus and microphthalmos head up the table, as does daily cannabis exposure interpolated. Furthermore, many of the E-values reported are very high, especially when it is noted that values above nine are said to be high [67].
When these 50 E-value estimates are listed out in order as shown in Table 7, 44/50 (88.0%) are in the high range [67] and 50/50 (100%) exceed the threshold of causality [66].  Supplementary Table S21 lists these E-values in order of the anomaly to which they apply. These data are then summarized by anomaly in Table 8 where they are listed in descending order of their minimum E-values. Interestingly, all the median minimum E-values are in the high range from 8.79 to 1.04 × 10 29 . Income 0 (0, 0) 4.71 × 10 −5 Abbreviations-Please see Table 1.
Similarly, these values can be ordered by the regression term involved, grouped into the primary covariate of interest (herb or resin THC content or daily cannabis use, Supplementary Table S22), and summarized (Supplementary Table S23). Groups can be compared using the Wilcoxson test where it is noted that all the inter-group comparisons are statistically significant (Supplementary Table S24). Thus, the order of severity in these multivariable studies is daily cannabis use > herb THC concentration > resin THC concentration.

Main Results
Nine of the eleven NCAs identified were related to metrics of cannabis exposure on bivariate testing and these relationships held through multivariable modelling for all six of the NCAs selected for advanced analysis. These overall results were consistent with other similar recent epidemiological reports in other jurisdictions [6][7][8][9][10][11].
Mapping studies showed that rates of most disorders were relatively severe in France. Rates of most disorders deteriorated in Spain (Figures 3 and 4 and Supplementary Figures S4 and S5). On bivariate maps, graphs with last month cannabis use × resin THC concentration as the independent variable, rates of nervous system disorders deteriorated in the Netherlands, Spain, France and Bulgaria; rates of severe microcephaly deteriorated in the Netherlands, Spain and France; rates of microphthalmos deteriorated in the Netherlands, Bulgaria and France; and rates of neural tube defects deteriorated in the Netherlands, Bulgaria and France (Figures 6 and 7 and Supplementary Figures S7 and S8).
Neurological congenital anomaly rates were higher in nations with increasing daily cannabis use compared to those without increasing daily use (p = 0.0204, minimum E-value (mEV) = 1.35; Supplementary Figure S9). Daily cannabis use interpolated had the most significant slopes of all the covariates (Figures 1 and 2, Supplementary Figures S1 and S2, Table 1). NCAs most closely related to cannabis use were severe microcephalus, nervous system disorders, eye disorders and microphthalmos ( Table 1).
Considering E-value estimates and lower bounds, 88.0% and 72.0% of 50 values exceeded nine and were thus in the high zone [67], and 100.0% and 96.0% of lower E-values exceeded the threshold for causality at 1.25 [66] (Table 7).
The order of importance of the primary cannabis metrics was daily cannabis use > herb THC concentration > resin THC concentration (Supplementary Tables S23 and S24).

Choice of Anomalies
The six anomalies analyzed in detail were selected to enable more comprehensive assessment of nervous system disorders overall. Anencephalus has been previously identified by researchers form the Centres for disease control [7] and it was clearly of interest to investigate if this finding would be confirmed in an independent dataset. Severe microcephaly was chosen as it arguably falls on the developmental pathway between anencephalus, reduced head size [17,18] and autistic spectrum disorder, which has previously been noted from prenatal cannabis exposure [14][15][16]69,70]. Further, as this anomaly classification does not exist in the US nomenclature it was of great interest to study it in the European context. Microphthalmia was of interest as it was first identified in the Hawaiian series but not found to be significantly cannabis-associated. It was also identified in a recent US series where it was found to be weakly cannabis associated. The signal in the USA data for both microtia and orofacial clefts was very strong so it was naturally of interest to see how the signal for microphthalmia would perform in the European data. Neural tube anomalies including spina bifida were of great interest as anencephalus is included in this group and had been identified with cannabis exposure, and neural tube anomalies generally had been positively linked with community cannabis exposure in Canada [8].

Qualitative Causal Inference
In 1965, the great epidemiologist A.B. hill described nine criteria which have become standard qualitative criteria against which to judge potentially causal relationships [71]. These criteria were strength of association, consistency amongst studies, specificity, temporality, coherence with known data, biological plausibility, biological dose-response curve, analogy with similar situations elsewhere and experimental confirmation. It is clear that all of these criteria are well fulfilled by the present results. It has been known for several decades that severe central nervous system anomalies were observed in the offspring of animals fed significant doses of cannabinoids during gestation. These defects included anencephalus and hydrocephalus [3][4][5]. Some of the genomic and epigenomic laboratory and clinical evidence is described below.

Quantitative Causal Inference
One of the typical concerns in generalizing the results of observational studies is that the underlying subgroups are not comparable. This issue is well addressed through the technique of choice in causal inferential studies: inverse probability weighting. This technique was deployed in all panel regression studies in this series which effectively overcomes this major interpretational concern and transforms these results from those from an observational study into a pseudorandomized study from which it is entirely proper to draw causal inferences.
The second major concern of such studies is that some extraneous covariate might exist which may account for the effects observed and negate the findings. Use of the E-values sets tight constraints on the associational behavior of such a hypothetical covariate by defining the degree of association required with both the outcome of interest and the exposure of concern of any such extraneous uncontrolled confounder. E-values above nine are generally said to be high [67]. E-values above 1.25 are usually required for causal effects [66]. As noted, 88% of the E-value estimates in the present study were in the high range so that uncontrolled confounding can be effectively discounted in interpreting the present results. Indeed, if the results reported in Supplementary Table S21 are examined closely, it is noted that the median minimum E-value of all six anomalies studied in detail was above 8.7, again giving great assurance that the results reported are indeed robust.

Morphogen Gradients
Forebrain specification is under the ventralizing control of sonic hedgehog (shh) and fibroblast growth factor 8 (FGF8) released from midline axial structures and opposed by GLI3 dorsally. shh specifies the early forebrain and divides the visual fields into left and right halves [41]. shh greatly stimulates forebrain growth. If the shh signal is blocked, forebrain growth and development is greatly inhibited and the mid-face region does not grow normally. In severe cases, this can lead to holoprosencephaly and sometimes cyclopia (single eye). These data link congenital anomalies affecting the facial region to brain development [41].
Spinal cord formation happens under the tight control of morphogen gradients acting in each of the three dimensional planes [41]. Dorsally high gradients of bone morphogenetic proteins (BMP) and Wnts exist which decline ventrally. Ventrally high gradients of sonic hedgehog exist which decline dorsally. The exterior of the spinal cord releases high gradients of retinoic acid which decline centrally towards the spinal canal. Caudally, there are high gradients of FGF8 and GDF11. Cephalad, there are high gradients of retinoic acid. These morphogen gradients control various gene cassettes (Pax genes dorsally and Nkx genes ventrally) to specific cells exactly in their correct positions including the well known dorsal-ventral split between sensory and motor neurons. Longitudinally, patterns of HOXA and HOXB genes control spinal development.
Since forebrain development is closely linked with shh signaling [41], shh inhibition by ∆9THC or cannabidiol [29] necessarily implies significant and severe impairment of brain development.
Cannabinoids can also interfere with morphogenetic signaling via epigenomic pathways as described below.
The pattern which emerges from cannabinoid inhibition of these guiding morphogens will depend on the timing and dose exposure of the developing fetus to the xenoteratogen. As noted, the potential implications of disruption of these key controllers of embryological morphogenesis are very serious.

Epigenomic Controls of Brain formation and Neurological Development
The recent paper by Schrott and colleagues in 2021 identified several pathways as importantly disrupted in cannabis dependence in rats and humans including cerebral disorder, neurodevelopmental disorder, agenesis and organismal growth [86]. After 11 weeks of cannabis withdrawal, which is one human sperm cycle in duration, they noted that pathways in hippocampal formation, quantity of pyramidal neurons, organismal death, cognitive impairment, encephalopathy and learning disabilities were increased. Importantly, cannabis dependence and withdrawal were associated with strong epigenomic links with many key neurotransmitter receptors of both the adult and developing brain including (Table 9; data from [86]) the following: GRIA (Glutamate Ionotropic AMPA Receptor) the predominant excitatory receptor of the brain also involved in longterm depression and potentiation in the hippocampus, mental retardation, depression and neuropathic pain syndromes [87]; GRIK (Glutamate Ionotropic Kainate Receptor) which has role in mental retardation, epilepsy and cold sensation [88]; GRIN (Glutamate Ionotropic NMDA Receptor) which plays critical roles in synaptic plasticity, long-term potentiation and depression and in mental retardation, chronic pain, motor movement and convulsive syndromes and schizophrenia, and Alzheimer's and Parkinson's diseases [89]; GRID (Glutamate Ionotropic Delta Type Receptor) is expressed in Purkinje cells of the cerebellum and throughout the brain and when mutated is involved in ataxic movement disorders [90]; GRM (Glutamate Metobotropic Receptor) is implicated in many diseases including schizophrenia, chronic pain syndromes, bipolar disorder, epilepsy and ADHD [91]. Other positive hits include the following: GABR (Gamma-aminobutyric acid receptor) forms the main group of inhibitory receptors in the brain which fine tune and control the action of excitatory receptors. It is also involved in synaptogenesis [92]; GABRA (Gamma-aminobutyric acid receptor Subunit A) is involved in numerous brain functions and importantly in epilepsy [93]; GABRB (Gamma-aminobutyric acid receptor Subunit B) is involved in many brain functions including epilepsy, autism, taste and somatosensory sensation, as a histamine receptor and in synaptogenesis [94]; HTR (5-Hydroxytryptamine Receptors) is involved in many disorders such as depression, schizophrenia, obsessive compulsive disorder, mood, perception, cognition and interaction with many psychoactive substances [95]; DRD (Dopamine Receptor Genes) is related to schizophrenia, addictions, endocrine actions, tremor [96]; DAT1/SLC6A3 (Dopamine Transporter/Solute Carrier Family 6, member 3) is involved in autism, ADHD, depression, bipolar disorder, schizophrenia, epilepsy, addictions, Parkinsonism and movement disorders [97]; OPRM (Opioid Receptor Mu) is the principal target of opioids, endorphins and enkephalins and is also involved in many other addictive drugs actions [98] (it is also involved in heterooligomerization and synaptogenesis); OPRD (Opioid Receptor Delta) is the receptor for enkephalins and is involved in transducing pain perception and immune signaling of certain interleukins (4 and 13) (it is found in the nucleus accumbens, olfactory bulb and neocortex [99]); OXTR (Oxytocin receptor) is involved in parturition and facial recognition, and via projections to the ventral tegmental area, pro-social and pro-sexual stimuli [100]; NRXN (Neurexin Receptor) cell surface ligand of neurexin family ligands includes cerebellins 1 and 2 and forms a key part of synaptic scaffold machinery (it exists in over 3000 spliced isoforms [101] and may be involved in schizophrenia); NLGN (Neuroligin) is a neuronal cell surface protein involved in synaptic formation with neurexins (it is implicated in schizophrenia, autism, mental retardation and pervasive developmental delay and nervous system development [102]); the neurexin-neuroligin ligand-receptor complex is known to play a key role in scaffolding of synapses, synapse formation and synaptic maintenance [103][104][105][106][107]. This system has also been shown to be inhibited by cannabis both directly [103,108] and epigenomically [86].
Details of the epigenomic DNA methylation changes at the opioid receptors are shown in Supplementary Table S25 (data from [86]). Changes are noted to be particularly marked in cannabis dependence.

Downs Syndrome Cell Adhesion Molecule
Supplementary Table S26 (data from [86]) lists 14 epigenomic change annotations for in Downs Syndrome Cell Adhesion Molecule (DSCAM) and in its closely related homologue DSCAM-Like 1 (DSCAML1), proteins which have previously been implicated in central and peripheral nervous system development and involved in neurite guidance and selfavoidance in the brain and retina, and axonal crossing of the midline in the spinal cord [109]. The protein has numerous isoforms and splice variants which allows unique self-coding for cortical and other neurons and thus forms the substrate for self-avoidance guidance cues.

Discs-Large Associated Protein 2
Discs-large associated protein 2 (DLGAP2) is a synaptic protein previously linked with schizophrenia and autism which has been found to be epigenomically altered in cannabis dependence with changes heritable through the paternal line of rats into the nucleus accumbens of offspring [70]. It was also found to be differentially hypomethylated at 17 sites in human sperm [86]. In the present dataset, it was identified at high levels of statistical significance at eight points, as detailed in Supplementary Table S27 (data from [86]).

Retinoic Acid in Forebrain Development-Direct and Epigenomic Effects
A powerful recent study comparing forebrain development in mouse, macaque and human embryos identified strong gradients of retinoic acid (RA) declining from the frontal pole to the premotor cortex as being primarily responsible for the surge in prefrontal cortex (PFC) neuronal numbers, spinogenesis, synaptogenesis, layer 4 lamination with its biomarker RORB and long-range connections to the mediodorsal thalamic nuclei [110]. The RA gradient was maintained by local synthesis by ALDH1A1, transduced by RXR-A, RXR-G and RARB-B receptors and RA was then catabolized towards the premotor cortex by CYP26B1.
As noted above, direct evidence for interaction by cannabinoids with RA has been previously described [81][82][83].
Moreover, there were annotations in the epigenomic study by Schrott and colleagues [86] describing two hits for ALDH1A1 one each in dependence and withdrawal with Bonferroni-adjusted p-values of 0.0106 and 0.0291, and nine other hits for members of the ALDH1 family in cannabis dependence and withdrawal. The RA-associated genes cadherin 8 (CDH8) and protocadherin 17 (PCDH17) had one and 156 hits, respectively (Supplementary Table S28, data from [86]). There were five hits for RARB (Bonferroni-adjusted p-values from 0.0006 to 0.0248) and one hit each for RXRA, RXRG, RXRA and RORB in cannabis dependence or withdrawal with Bonferroni-adjusted p-values from 0.0009 to 0.0191 (Supplementary Table S29, data from [86]).
Moreover, key amongst the secondary messengers induced by RA in the fetal middle trimester was cerebellin 2 (cbln2), which mediated the massive growth in dendritic spines and synaptogenesis of the human forebrain [111]. RA was found to be a transcription factor acting on the cbln2 enhancer. Cerebellins ligate neurexins 1, 2 and 3 and GRID 1 and 2 which were all highly expressed during midfetal forebrain development. Cannabis has also previously been shown to interfere with neurexin binding [105,107].

Cerebellins
Cerebellin 2 (cbln2) has previously been linked with diverse neurological syndromes including obsessive compulsive disorder, Tourette's syndrome and schizophrenia and the binding partners of cbln2, neurexins and GRID2 have been linked with autism and schizophrenia [111].
In the epigenomic screen of Schrott, there was one hit for a cbln2 intron in cannabis withdrawal on page 183 with a Bonferroni-adjusted p-value of 0.01723 (page 183) [86]. Interestingly, there were 71 hits for cbln4 in from Bonferroni-adjusted p-value = 0.0086 (cannabis withdrawal, page 29) and from p = 7.73 × 10 −20 (page 239 in cannabis dependence). Cbln4 has a different distribution in the brain to cbln2 and mainly occupies subcortical and hindbrain sites [112].

Slit-Robo
The robo2-slit1 receptor-ligand system has previously been shown to be key to the exuberant development of the massive human neocortex characterized by very high numbers of neurons and synaptic connections [113,114]. These results were validated by the recent demonstration of the importance of the slit-robo Rho GTPase Activating Protein 2C (SRGAP2C) in a comparative neurodevelopmental study where overexpression of human SRGAP2C in mice led to increased frontal lobe connectivity and improved learning of a task by the modified animals [115]. The slit-robo system can be directly antagonized by cannabinoids [116,117]. Examining both SRGAP2C and its antagonist SRGAP2B in the Schrott data provided the detailed results showing that cannabis also antagonizes this system epigenomically (Supplementary Table S30, data from [86]). Data on annotations relating to slit and robo are reported in Supplementary Table S31 and Table 10, respectively; data from [86].
Of greater concern is that discrete jumps in the incidence of many congenital anomalies have been described for many anomalies at the highest level of cannabis exposure [1,10,11]. This important finding indicates directly that the concerning findings repeatedly demonstrated in the laboratory are in fact confirmed in the profiles of public health and neonatal seen epidemiologically.

Generalizability
This study is remarkable for its notable consistency with previous reports in animals [3][4][5] and with other epidemiological reports from Australia, Canada and the USA [8][9][10]119,120]. The results also greatly strengthen earlier bivariate reports from Europe and confirm results in the space-time context and formal causal inferential frameworks [119,121]. That is to say that there is an impressive uniformity of many major datasets in various countries around the world which generally indicate similar results. Together, this impressive body of evidence presents a strong and uniform set of results linking cannabinoid exposure with NCA teratological outcomes. A causal relationship between cannabis exposure, on the one hand, and NCA teratology, on the other, is further supported by the many mechanistic pathways which form a key and pivotal plank of formal causal algorithms of establishing causality [71] and further indicating the causal nature of this link.
Of particular importance in this discussion is the evidently graded series of NCAs from mildly impaired intellectual development to microcephaly, to severe microcephaly to anencephaly which demonstrates a clear and graded spectrum of adverse effects of cannabinoid exposure on infant and child neurological development.
In view of this impressive and remarkably uniform body of work, it seems clear that the relationship described is causal in nature. Because of this notable confluence of large epidemiological datasets on this issue, the strength of both the statistical and mechanistic arguments for causality and the concordance with data from preclinical models, we feel that this relationship is widely generalizable to other places wherever data of sufficient size and quality exist for it to be reliably assessed.

Strengths and Limitations
This study had a number of strengths including the use of one of the world's largest and most comprehensive congenital anomaly datasets; the use of advanced statistical modelling; the use of inverse probability weighting and E-values to engage the techniques for formal quantitative causal inference and transfer the analysis from an observational study to a pseudorandomized study; the liberal use of multi-paneled maps and graphs to display multiple covariates across time and space; the use of bivariate maps, which is unusual, and the provision of extensive supplementary material, online resources, data and computational code. Ranger regression was used for formal variable selection. Like many epidemiological studies, the present work did not have access to individual cannabis exposure data. Additionally, some of the data, particularly those relating to daily cannabis use, had to be interpolated due to severe missing data problems. This limitation should be borne in mind when interpreting results for relevant parameters.

Future Directions
The present strong and now widely replicated results reported herein and in similar comparable studies internationally clearly place a major spotlight on the whole issue of the impacts of cannabinoids on neural development. Hence, the impact of the present report is to increase the research importance of these issues and justify further detailed studies in these areas. Further fine-resolution spatiotemporal and causal inferential epidemiological studies need to be performed. Mechanistic and particularly epigenomic studies need to be explored in detail to further understand these effects and are all strongly indicated.

Conclusions
This study found that nine of the eleven NCAs in this European dataset were closely related to metrics of cannabis exposure on bivariate analysis and this relationship persisted after adjustment in all six of the NCAs selected for further detailed multivariable analysis. In this regard, these results are consistent with those from Canada, Hawaii, Colorado and the USA [1,10,11] and also with previous reports from CDC [7]. The major tools of formal quantitative causal inference were widely employed in this study. Inverse probability weighting of all panel models transferred the findings from those of a merely observational study into a pseudorandomized study from which it is quite appropriate to draw casual inferences. Collected E-values were predominantly in the high range; effectively excluding uncontrolled confounding was a possible explanation for these results. Therefore, these results fulfilled quantitative epidemiological criteria of causality. Beyond the very concerning findings of a number of NCAs closely related to metrics of cannabis exposure, findings raise two particular concerns. The first relates to interference with brain development during in utero life, which is typically a long-term disability from which recovery is difficult and the degree of disability severe. The second relates to the now well documented exponential dose-response relationships of cannabinoids both in the laboratory [22][23][24][25][26]30,35] and in a number of recent epidemiological studies [1,10,11]. Both these concerns strongly indicate that penetration of cannabinoids into the community should rationally be tightly restricted if we are to seriously steward our responsibilities as custodians of the human brain, genome and epigenome for the generations to come. Author Contributions: 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, 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. All authors have read and agreed to the published version of the manuscript.

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