Epidemiological Patterns of Cannabis- and Substance- Related Congenital Uronephrological Anomalies in Europe: Geospatiotemporal and Causal Inferential Study

Introduction. Recent reports linking prenatal and community cannabis exposure to elevated uronephrological congenital anomaly (UCA) rates (UCAR’s) raise the question of its European epidemiology given recent increases in community cannabinoid penetration there. Methods. UCAR data from Eurocat. Drug use data from European Monitoring Centre for Drugs and Drug Addiction. Income from World bank. Results. UCAR increased across Spain, Netherlands, Poland and France. UCAR’s and cannabis resin THC increased simultaneously in France, Spain, Netherlands and Bulgaria. At bivariate analysis all UCA’s were related to cannabis herb and resin THC concentrations. All UCAR’s were bivariately related to cannabis metrics ordered by median minimum E-value (mEV) as hypospadias > multicystic renal disease > bilateral renal agenesis > UCA’s > hydronephrosis > posterior urethral valve > bladder exstrophy/epispadias. At inverse probability weighted multivariable analysis terms including cannabis were significant for the following series of anomalies: UCA’s, multicystic renal disease, bilateral renal agenesis, hydronephrosis, congenital posterior urethral valves from P = 1.91 × 10−5, 2.61 × 10−8, 4.60 × 10−15, 4.60 × 10−15 and 2.66 × 10−10. At geospatial analysis the same series of UCA’s were significantly related to cannabis from P = 7.84 × 10−15, 7.72 × 10−5, 0.0023, 6.95 × 10−5, and 8.82 × 10−5. 45/51 (88.2%) of E-value estimates and 31/51 (60.8%) of mEV’s >9. Conclusion. Analysis confirms a close relationship between cannabis metrics and all seven UCA’s and fulfill formal criteria for quantitative causal inference. Given the exponential cannabinoid genotoxicity dose–response relationship results provide a powerful stimulus to constrain community cannabinoid exposure including protection of the food chain to preserve the genome and epigenome of coming generations.


Background
Recent reports demonstrate increasing concern at the genotoxic activities of cannabinoids expressed as congenital anomalies, cancer and aging . Recent epidemiological teratological reports have causally linked cannabis exposure with cardiovascular, chromosomal, central nervous system, limb, gastrointestinal, orofacial and body wall congenital anomalies [6][7][8]12,15,18,[24][25][26][27][28][29][30][31]. A prominent part of cannabis teratology is uronephrological (urinary tract and renal) congenital anomalies (UCA's) which have featured in reports of elevated rates of renal agenesis in Colorado [13], obstructive genitourinary defect in Hawaii [19], congenital posterior urethral valve, obstructive genitourinary defect, renal agenesis/hypoplasia, hypospadias and epispadias in USA [7] and hypospadias in Australia [8]. It was therefore of interest to investigate contemporary trends of this group of anomalies in Europe particularly given the rising triple convergence of prevalence of daily use, intensity of daily use and rising cannabinoid potency in widely available products in many countries across the European continent in the recent decade [20,32].

Exponential Dose-Response Effect Curve
It is important to appreciate that the genotoxic action of cannabinoids is strongly exponential with effects rising dramatically in the higher dose range. This has been noted both for the direct DNA mutagenicity [43,[65][66][67][68][69][70][71][72] and for the mitochondrial metabolic reactions which underpin genomic and epigenomic reactions and provide substrates and energy to the genetic and epigenetic machinery [67,[73][74][75][76][77][78]. This implies in turn that the rising triple confluence of prevalence-intensity-potency mentioned above can relatively abruptly place communities into a zone of high genotoxic risk result in a relatively sudden appearance of major genotoxic outcomes [20,32]. This likely accounts for adverse experiences such as recent French and German outbreaks of limblessness at up to sixty times the background rate in areas where many hectares are sown to cannabis [79][80][81][82]. Interestingly not only are babies being born without limbs but so too are the calves suggesting cannabinoid exposure via the food chain [79][80][81]. However, none of these effects are happening in nearby Switzerland where cannabinoids are not permitted in foodstuffs.

Lessons from VACTERL Syndrome
Interestingly a companion paper to this one has shown that the VACTERL syndrome (Vertebral, anorectal, cardiac, tracheo-esophageal fistulae/esophageal atresia, renal and limb anomalies) [83] was strongly and causally linked with European cannabinoid exposure [15,84]. Since renal anomalies are part of the VACTERL syndrome this finding also imputes renal abnormalities in the overall pattern of uronephrological congenital anomalies (UCA's).

Study Questions
For these reasons it was considered important to assess UCA rates (UCAR's) in response to increased cannabinoid exposure across many European countries, particularly as UCA's appear to be a relatively sensitive mutagenic index of cannabinoid teratology. It was also considered important that this assessment be undertaken in a bivariate and multivariable context within a causal inferential framework, and in a spatiotemporal context which formally takes account of the native space-time environment in which this data was generated.

Data
The data on congenital anomaly rates was downloaded from the European Network of Population-Based Registries for the Epidemiological Surveillance of Congenital Anomalies (EUROCAT) website [85]. Annual data for all available anomalies was downloaded from both the summary and the detailed files for 14 nations for each individual year. The metric of interest was the total congenital anomaly rate which includes live births, stillbirths and cases where early termination for anomaly was practised. The total rate therefore represents the rate net of all the different birth exigencies to which pregnancies are subject. Nations were selected based upon most of their congenital anomaly data being available across most of the period 2010-2019. Data on national tobacco (percent daily tobacco use prevalence) and alcohol (litres of pure alcohol consumed per capita annually) consumption were obtained from the World Health Organization website [86]. The website of the European Monitoring Centre for Drugs and Drug Addiction (EMCDDA) [87] was used to source data on drug use relating to cannabis, amphetamines and cocaine. The major metric of cannabis exposure was last month cannabis use data. This was supplemented by data on the tetrahydrocannabinol (THC) content of cannabis herb and resin provided in recent published reports [32] and by data on daily cannabis use which was also available from EMCDDA and was collated in recent Public Health reports [32]. Data on median household income data (in $USD) was derived from the World Bank sources [88].

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

Derived Data
Because several metrics of cannabis exposure were available it was possible to combine and group these in various ways. Thus last month cannabis use prevalence data was multiplied by the THC content of cannabis resin and herb to derive compound metrics from their product. These metrics were again multiplied by imputed daily cannabis use prevalence rates to derive more complex compound metrics for both cannabis resin and herb.

Data Imputation
Linear interpolation across years was used to complete missing data. This was especially so for daily cannabis use. For these 14 nations 59 data points on daily cannabis use were available from EMCDDA across this period. Linear interpolation was able to expand this dataset out to 129 datapoints (further details are provided within Results Section). There was no available data on cannabis resin THC concentration for Sweden. It was noted however that the resin: herb THC concentration was virtually in nearby Norway at 17.7. This ratio was therefore applied to the Swedish cannabis herb THC concentration data to derive estimates of the Swedish cannabis resin THC concentration. In a similar manner data for the cannabis resin THC concentration in Poland were not available from EMCDDA. The resin to herb THC concentration ratio of nearby Germany across time was used to estimate the resin THC content in Poland from the known Polish herb THC concentrations. Since geospatial analytical techniques do not tolerate missing data the dataset was completed by interpolation by the last observation carried forward or backwards for Croatia in 2018 and 2019 and Netherlands in 2010. Multiple imputation methods could not be used for these analyses as multiple imputation techniques have not been applied to panel or spatial multivariable regression methods.

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 was used for data processing [89]. The analysis was performed in December 2021. Data manipulation was performed using dplyr from the R Core team tidyverse [90]. Data were log transformed where appropriate as guided by the Shapiro-Wilks test to improve compliance with normality assumptions. ggplot2 from tidyverse was used for constructing graphs. Both ggplot2 and sf (simple features) [91] were used for drawing maps using both custom colour palettes and palettes provided in the viridis and viridisLite packages [92]. Bivariate maps were drawn with package colorplaner by William Murphy [93]. Illustrations are all original and have not been previously published. Linear regression was conducted in R Base. R package nlme was used for mixed effects regression [94]. Model reduction was by the classical technique of serial deletion of the least significant term to yield a final reduced model which is the model presented in all multivariable models. Combined techniques from R packages purrr and broom were used to process multiple linear models in a single pass virtually simultaneously [90,95,96]. The overall effect of covariates particularly of interactive multivariable models may be quantified and is denoted as the marginal effect. The R package margins was used to calculate the overall marginal effect [97].

Covariate Selection
With numerous different metrics of cannabis exposure available it was unclear which was the most appropriate set of covariates to select for each individual model. Use of excessive cannabis related covariates could lead to issues of collinearity, consumption of excessive degrees of freedom and a concomitant limitation of interactions which could be considered, or model mis-specification or imbalancing more generally.
The presence of multiple different metrics for cannabis consumption and exposure created a problem for analysis as it was not clear which was the most appropriate metric to employ for any particular model. Indiscriminate use of excessive covariates in a multivariable model would unnecessarily consume degrees of freedom and thereby restrict ability to assess interactions. Random forest regression using the R package ranger [98] with variable importance being formally assessed via the R package vip (variable importance plot) [99] was used to formally address this issue of variable selection. The most predictive of the covariates identified by this process were then entered into the regression modelling equations. Applicable tables from these analysis are presented below in the Results Section.

Panel and Geospatial Analysis
The R package plm [100] was utilized to conduct panel analysis across both space and time simultaneously using the "twoways" effect. The spatial weights matrix for temporospatial modelling was calculated using the edge and corner "queen" relationships using R package spdep (spatial dependency) [101]. The spatial panel random effects maximum likelihood (spreml) function from the package spml which allows detailed modelling and correction of model error structures [102,103] was used for geospatial modelling. Such models can produce up to four descriptive model coefficients which are can be used to verify the most appropriate error structure for the model. These coefficients are psi the serial correlation effect, rho the spatial coefficient, phi the random error effect and theta the spatial autocorrelation coefficient. The most appropriate error structure was chosen for each series of spatial models taking care to preserve the model error specification across related models. The appropriate error structure was determined by the backwards method from the full general model to the most specific model by the method which has previously been published [104]. Temporal lagging by one or two years was employed as indicated in both panel and geospatial models.

Causal Inference
This analysis deployed the formal tools of causal inference. Inverse probability weighting (ipw) is the analytical technique of choice to convert a purely observational study into a pseudo-randomised study from which it is entirely appropriate to draw causal inferences [105]. Inverse probability weighting was employed in all multivariable panel models presented. Inverse probability weights were calculated via the R package ipw. Similarly E-values (expected values) 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 relationship which appears to be apparently causal in nature [106][107][108]. It therefore provides a quantitative measure of the robustness of the model to extraneous covariates which have not been accounted for within the measured parameters and is a popular and increasingly published form of sensitivity analysis. E-Values have a measure o uncertainty associated with them which can be quantified as a confidence interval and the 95% lower bound of this confidence interval is particularly relevant to situations of increased risk. For these reasons it has been extensively reported herein. E-Value estimates greater than 1.25 to indicate causality [109] and E-values larger than nine are described as being high [110]. The R package EValue was used for the calculation of both E-value estimates and their 95% lower bounds [111]. Both inverse probability weighting and E-values are foundational and pivotal techniques used in formal quantitative causal inferential methods in order to allow causal relationships to be assessed from real world observational studies and to control for the effect of observed and unobserved covariates respectively.

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 URL's: https://data.mende ley.com/datasets/c6psrbr34j/1 and https://data.mendeley.com/datasets/vd6mt5r5jm/1.

Ethics
The Human Research Ethics Committee of the University of Western Australia provided ethical approval for the study to be undertaken 24 September 2021 (No. RA/4/20/4724). All methods were carried out in accordance with relevant guidelines and regulations and were concordant with the Declaration of Helsinki.

Input Data
Supplementary Table S1 sets out the general overview of the covariates used in this study. As indicated 839 data points were downloaded from the EUROCAT database from the 14 nations shown and applying to the seven CA's of urinary disorders, bilateral renal agenesis, bladder exstrophy or epispadias, congenital hydronephrosis, multicystic renal dysplasia and posterior urethral valve or prune belly. Covariates for drug exposure including various primary and compound metrics of cannabis exposure are as indicted in the table along with median household income data.
Data on daily drug use by European nation and year was examined and was found to be incomplete as indicated in Supplementary Table S2 with 59 points available. In order to allow this information to be used in comparative studies a further 70 datapoints were added to this data by linear interpolation. This augmented dataset is shown in Supplementary Table S3. Figure 1 shows the relationship between the rates of the various UCA's and the substance tobacco, alcohol, daily cannabis use interpolated, amphetamine and cocaine exposure. It appears that tobacco exposure is not associated with any UCA incidence. There does seem to be a positive relationship between annual alcohol consumption and multicystic renal disease and posterior urethral valve as those regression lines have a positive slope. Amphetamine and cocaine exposure appear to be positively related to all defects mentioned. Daily cannabis interpolated appears to be related to hydronephrosis, multicystic renal disease and urinary anomalies globally. Figure 2 illustrates the relationship between various cannabis exposure metrics and rates of the various UCA's. Many of the regression lines appear to be positively up-sloping including those for last month cannabis use, the THC concentration of cannabis herb and resin, and most of the compound cannabis metrics. Figure 3 shows the rates of UCA's across Europe for the past decade. Rates seem to have risen in Spain, Portugal, France, Belgium, Bulgaria, Sweden and Croatia but have declined in Germany.

Bivariate Analysis
Rates of bilateral renal agenesis are shown in Figure 4. The appear to have increased in Spain, Norway, Germany, Netherlands, Belgium and Croatia but reduced in Poland.
Rates of posterior urethral valve are shown in Figure 5. Rates increased in Spain, Portugal, Poland, France, Bulgaria and Sweden but reduced in Germany. Rates in the low countries and France fluctuated.
Rates of multicystic renal disease are shown in Figure 6. Rates increased in France, Bulgaria, Norway, Poland, Netherlands, Belgium, Croatia and Sweden and fluctuated in Germany. Figure 7 shows the rate of the compound metric last month cannabis use: cannabis resin THC concentration. It has increased across the continent but the rises in Spain, Portugal, Netherland and France were most accentuated.              Figure 8 is a bivariate plot of the co-distribution of uronephrolgical congenital disorders and last month cannabis use: cannabis resin THC concentration across Europe over the decade. One reads the plot by noting the areas which turn purple or crimson which signify a high distribution of both covariates. The significance of the other colours is indicated in the colorplane key. The series of graphs makes clear the co-emergence of high rates of both uronephrolgical disease and the compound cannabis use metric in France, Netherlands and Bulgaria across the decade. Figure 9 is a similar bivariate series of maps showing the co-emergence of high rates of bilateral renal agenesis and the compound cannabis exposure metric in Netherlands. Figure 10 plays the same role for rates of posterior urethral valve. The co-emergence of high rates in France and Netherlands is made plain by this graphic. Figure 11 performs a similar function for multicystic renal disease. The co-emergence of high rates in France, Netherlands and Bulgaria is shown.

Bivariate Overview
Nations were grouped into those where daily cannabis use was increasing or high and those where it is not based on recent epidemiological reports [32]. When all CA's are groups together the appearances seen in Figure 12 are shown. Nations where daily use is increasing are shown to have higher rates of UCA's at linear regression (β-est. = 0.3121, t = 2.687, p = 0.0073; model Adj.R.Squ. = 0.074, F = 75.22, df = 12, 837, p = 0074). The E-value estimate for this effect is 1.71 with a 95% lower bound on the confidence interval (mEV) of 1.29 which exceeds the quoted threshold fror causality of 1.25 [109] When this comparison is done in an anomaly-specific manner the appearances shown in Figure 13 are seen. When these differences are studied formally by mixed effects regression the rate in the countries with increasing daily use is significantly higher than the remainder (β-est. = 0.2048, t = 3.943, p = 8.73 × 10 −5 ; model AIC = 1755.781, Log.Lik = −872.891). This effect is associated with an E-value estimate of 1.97 and mEV = 1.57.

Bivariate Linear Regressions
The regression lines shown graphically in Figures 1 and 2 may be studied by multiple bivariate linear regression models simultaneously using the purrr-broom workflow in R. when this is performed 85 regression models are derived which are listed in Supplementary Table S2 in descending order of mEV. Interestingly daily cannabis use appears at the top of this list of models.
From this list those which have positive regression coefficients and significant p-values may be extracted. 37 such models are retained in this way which are listed in Table 1.

Covariate Selection
Having identified so many bivariate relationships of interest the next question relates to how these different covariates might perform in a multivariable context. However, given that there are 13 substance exposure covariates it is not entirely obvious which might be the best combination of covariates to use in anomaly-specific regression equations.
This issue was addressed formally by the use of random Forest regression in tandem with variable importance assessment using the R packages ranger and vip. When this was done for the overall UCA's and for multicystic renal disease, bilateral renal agenesis, hydronephrosis and posterior urethral valve the variable importance tables shown in Supplementary Tables S3-S7 were derived.

IPW Panel Regression
These covariates were then used to define multivariable inverse probability weighted panel regression models for a series of multivariable models relating to the five congenital uronephrological disorders of interest. The inverse probability weighting is of considerable importance and allows the analysis to move from that of an observational study only to a pseudo-randomized milieu from which is it appropriate to draw causal inferences.
Supplementary Table S8 describes three such final inverse probability weighted panel regression models for uronephrological disease as a whole. Additive, interactive and temporally lagged models are shown. In each one cannabis metrics are significant in final models and survive model reduction, have positive regression coefficients and are highly statistically significant.
This pattern is continued across the other anomalies, namely multicystic renal dysplasia, bilateral renal agenesis, hydronephrosis and posterior urethral valve syndrome as shown in Supplementary Tables S9-S12.
This result is important because it means directly that in this set of congenital anomalies the strong trends identified at bivariate regression persist after multivariable adjustment in an inverse probability weighted and therefore casual paradigm.

Geospatial Regression
The next issue to consider was whether these relationships would continue when considered formally in their native space-time context in which major methodological concerns such as random errors, serial correlation, spatial correlation and spatial autocorrelation would be properly accounted for and managed.
Supplementary Figure S1 presents maps of the initial, edited and final geospatial links between countries which were derived and used as the basis for the sparse spatial weights matrix used in the spatial regressions. Table 2 presents additive, interactive and lagged geospatial models for uronephrological disorders broadly. The pattern identified earlier in panel regression is continued with some cannabis-related terms in each model persisting in the final model after model reduction and maintaining high levels of statistical significance from 7.84 × 10 −15 in the case of the first additive model. This pattern is continued when multicystic renal disease is considered as shown in Table 3. In each model the effect of cannabis-related covariates is overwhelmingly in the positive direction.
Similar results are also seen when bilateral renal agenesis is considered in interactive and lagged models as shown in Table 4.
When hydronephrosis is considered the pattern is continued in additive, interactive and models lagged to one year, but the signal fades by two years of temporal lagging ( Table 5).
The pattern is continued and confirmed when posterior urethral valve syndromes are considered as shown in Table 6. In each case terms including cannabis metrics persist after adjustment in final models, have positive regression coefficients and are statistically significant.

Causal Inference E-Values
E-values may be extracted from each of these regression models. Table 7 does this for panel models and Table 8 does this for spatial models.
Both lists are combined in Table 9 where 51 such E-value pairs are listed in descending order of mEV. The prominence of daily cannabis use interpolated at the head of the table is noted as is the high level of statistical significance of the E-values returned.                 The E-values may be extracted alone and listed and ordered as shown in Table 10. Here, it is noted that 45/51 (88.2%) E-value estimates exceed 9 and are thus in the high range [110], and all 51/51 (100%) exceed the threshold of causality at 1.25 [109]. Moreover, 31/51 (60.8%) mEV's exceed nine and therefore fall into the high range and 51/51 (100%) exceed the threshold for causality.   Table 9 may also be ordered by anomaly as shown in Table 11 and summarized in Table 12. Table 12 is also ordered by descending median minimum E-value. The table is headed up by hydronephrosis and posterior urethral valve. It is noted that the median mEV's are relatively high for the five anomalies studied. Table 9 may also be ordered by the regression term and grouped according to the major covariate between daily cannabis use, and the THC concentration of cannabis herb and resin. This is defined by the extra column which has now been added headed "Group" in Table 13. This table may then be summarized by the primary cannabis metric as shown in Table 14.
These results may then be compared using the Wilcoxson test a indicated in Table 15. One notes immediately that most of the apparent between group differences noted in Table 14 are statistically highly significantly different.

Main Results
The main results of this study is that all the UCAR's identified are closely related to various metrics of cannabis exposure. Strong relationships were demonstrated on bivariate analysis which persisted after multivariable adjustment in inverse probability weighted panel models and thereby fulfilled the formal criteria of quantitative causal inference. These results were also confirmed for all five UCA's studies by multivariable geospatial regression.

Detailed Main Results
A mapping analysis demonstrated that UCAR's increased across Spain, Netherlands, Poland and France. Congenital posterior urethral valve rates increased in Spain, France, Belgium, Poland and Netherlands. Multicystic renal disease increased in Spain, France, Netherlands and Norway. A bivariate map study showed that UCAR's and last month cannabis use x cannabis resin THC concentration increased simultaneously in France, Spain, Netherlands and Bulgaria. Bilateral renal agenesis increased along with last month cannabis use x cannabis resin THC concentration in Netherlands. Congenital posterior urethral valve increased along with last month cannabis use x cannabis resin THC concentration in Netherlands and France. Multicystic renal disease increased along with last month cannabis use x cannabis resin THC concentration in Netherlands and France. That is UCAR's increased in countries with increased daily cannabis exposure.
When nations with increasing daily cannabis use were compared with those without the former group has a higher rate of UCA's overall (p = 0.0073). When this comparison was made by UCA this relationship intensified (p = 8.73 × 10 −5 ).
At bivariate analysis tobacco exposure was not related to UCA's. Alcohol exposure was related to multicystic renal disease and posterior urethral valve. Cocaine and amphetamine were strongly related to most UCA's. All UCA's were strongly related to cannabis herb and resin THC concentrations (Figures 1 and 2). All anomalies could be related to cannabis metrics on bivariate analysis with strength of association judged by the median minimum E-value (mEV) as hypospadias (174.71) > multicystic renal disease (5.37) > bilateral renal agenesis (5.05) > UCA's (2.68) > hydronephrosis (2.21) > posterior urethral valve (1.73) > bladder exstrophy/epispadias (1.53) ( Table 1).

Choice of Anomalies
The selection of UCA's chosen for detailed analysis was based upon their implication in prior reports as mentioned in the Introduction, the frequency of the anomaly and the clinical significance. Hypospadias and epispadias were not chosen as they are relatively minor anomalies. Bladder exstrophy was not studied further as it is very rare.

Qualitative Causal Inference
It is worth considering how the present results perform against the Hill criteria of causality which are accepted criteria for assessing causal relationships [112]. These criteria were strength of association, consistency amongst studies, specificity, temporality, coherence with known data, biological plausibility, overresponse relationship, analogy with similar situations elsewhere and experimental confirmation. It is noted that the present study fulfills all of these criteria.

Quantitative Causal Inference
One of the major issues confronting the interpretation of observational studies is the issue of non-comparability of experimental groups. This issue can be addressed by inverse probability weighting which is the technique of choice in causal inference and has been widely adopted here in all panel models. This technique has the effect of transforming the analysis from that of an observational study only into a pseudorandomized quasiexperimental system from which it is entirely proper to draw causal inferences.
The other issue with which observational studies are frequently confronted is the potential for some uncontrolled extraneous confounding covariate to explain away the reported results. This issue is addressed by calculating the E-value which quantifies the degree of association required of such an hypothetical covariate with both the exposure of concern and the outcome of interest in order to obviate the reported apparently association. E-Values (expected values) in excess of nine are reported as being in the high range [110] and a cut-ff of 1.25 is usually applied as a threshold of a potentially causal effect [109]. Hence, the very elevated E-values reported in the present study are well above this range and qualify findings to be designated as causal in nature.
Together with highly concordant results reported from several other jurisdictions and the strong theoretical mechanistic basis of these findings the highly significant results reported herein exclude the possibility that the results reported may be due to either chance (p-values) or bias (E-values).
In part this effects happens directly but it is also partly epigenomically mediated. Hence, for example genes such as GLI3, (Gli family zinc finger 3), MEGF8 (multiple EGFlike domains 8), TMEM107 (transmembrane protein 107) and BMP4 (bone morphogenetic protein 4) which are frequently cited in the epigenome-wide screen as being altered in cannabis dependence and withdrawal, are known to function in or around or antago-nize sonic hedgehog signalling which is a critically important controller of embryonic morphogenesis across numerous tissues and organs.
This link between cannabidiol and prostate cancer has been confirmed in recent epidemiological reports [1][2][3].
This link between cannabidiol and ovarian cancer was confirmed in a recent epidemiological study [1][2][3].

Nephrogenic Genes in Renal Organoids
A recent detailed study of renal morphogenesis has been most revealing in terms of disclosing, and confirming genes critically involved in renal morphogenesis [125]. Inhibition of ROCK (Rho associated coiled-coil containing protein kinase) proved to be a key step in dramatically increasing the yield of embryonic stem cell differentiation into nephroblastic progenitor cells. Many genes involved in the induction of ciliopathies were also implicated in grossly abnormal nephrogenesis. Furthermore, the key notch ligand Jag1 (Jagged canonical notch ligand) was shown to potently inhibit differentiation into renal tubular cells. Table 17 lists several of the genes from the notch pathway which were identified in the epigenomic EWAS of Schrott and colleagues [55]. RBPJ (Recombination Signal Binding Protein for Immunoglobulin Kappa J Region) and PSENEN (Presenilin enhancer, gamma secretase subunit, a key transmembrane effector of notch ligand cleavage and thus signal transduction) were prominently identified as shown.
Sonic hedgehog (shh) is another key morphogen involved in nephrogenesis [125,126]. Table 18 lists some of the EWAS hits from the Schrott study where shh pathway members were identified as being disrupted by cannabis exposure. Patched is the canonical shh receptor. SUFU (SUFU negative regulator of hedgehog signalling) is an inhibitor, PSENEN also cleaves the shh ligand at the cell membrane and Gli3 (GLI family zinc finger 3) is a key DNA binding transcription factor of the shh pathway. As there were 185 EWAS hits for Gli3 only a sample can been displayed in the table. Table 19 provides a further list of key genes involved in nephrogenesis. The CCDC genes are coiled-coil domain containing proteins. CCDC170 was noted to be key for renal tubule cell formation as was MYH7 (Myosin heavy chain 7) and EPCAM (epithelial cell adhesion molecule), the bone morphogenetic proteins (BMP's) and the Wnt's. Indeed the Schrott data included 191 hits for BMP's and 203 hits for Wnt family members.     Table 20 and their very high level of statistical significance is evident.   TGFβ was also shown to be important in the renal organoid GWAS screen and induced interstitial fibrosis in line with its well known activities in many other tissues. Eight EWAS hits from the Schrott data for TGFβ are shown in Table 22. Downstream TGFβ signals are transduced by a cascade of SMAD (mothers against decapentaplegic homologue) intermediates and transcription factors. SMAD4 is the final transcription factor to enter the nucleus from this pathway. Cannabis exposure EWAS hits for this series of mediators is also shown in the table.

Exponential Genotoxic Effects
A wide variety of studies have documented the exponential dose-response relationship of mutagenicity and various indices of cannabinoids exposure [43,[65][66][67][68][69][70][71][72][73][74][75][76][77][78]. This applies both to directly genotoxic studies and to studies of the intermediary mitochondrial metabolism upon which genomic and epigenomic reactions are based both as substrates and for energy supply. Moreover, several reports confirm that the exponential effects seen in the laboratory are also documented in patterns of human disease, particularly by modelling studies and by showing an abrupt discontinuous step in disease rates from the fourth to the faith quintile of cannabinoid exposure. This exponential cannabinoid dose-response relationship where disease rates occur at higher levels of cannabinoid exposure is very important when interpreting epidemiological patterns of UCAR's. This is especially so for Europe where increases in cannabis use prevalence, cannabis use daily intensity, and THC concentration in cannabis herb and resin are all going on at the same time creating an exponential cannabinoid exposure curve. This in turn implies that whole nations can be relatively suddenly moved up into the genotoxic dose-response zone where major and severe genotoxic outcomes become commonplace. As was noted this may account for severe and apparently abrupt teratological outbreaks such as the limblessness phenomena described earlier.

Strengths and Limitations
This study has many strengths including the use of one of the largest UCAR databases in the world and the very comprehensive EMCDDA drug dataset. Moreover, advanced statistical modelling has been utilized in performing this analysis including inverse probability weighting, lagged models, the formal tools of quantitative causal inference and geospatiotemporal regression modelling. Panelled maps and graphs have been used to show all covariates across time at a single glace. Bivariate maps have also been used, which is unusual. Moreover, a sophisticate mechanistic understanding has been presented which not only explains the pattern of UCAR's observed but contributes pivotal information into the key issue of causal assignment. Study limitations include the lack of individual participant exposure information, a limitation which is shared in common with many epidemiological investigations. Moreover, missing data was an issue form some covariates particularly daily cannabis use and this must be borne in mind in considering the present results.

Generalizability
Results of this study are widely generalizable for several reasons. The analysis is based on one of the largest UCAR datasets in the world, and concordant with similar published results presented from other jurisdictions. Further, results are at high levels statistical significance and sensitivity analysis using E-values indicates that they are robust to extraneous explanations. Moreover, the analysis has been conducted within a formal quantitative causal inferential framework so that the effects demonstrated may be designated as causal effects. For all of these reasons we feel that these results are widely generalizable to other situations wherever data of sufficient quality exist to make the relevant assessments.

Conclusions
In conclusion this analysis confirms previous reports showing that many UCA's are related to various metrics of cannabis exposure. To those anomalies which have already been described this report adds congenital hydronephrosis, multicystic renal disease and urinary anomalies overall as newly cannabis associated congenital anomalies. This analysis has been conducted within a causal inferential framework including the use of inverse probability weighting and E-values which formally demonstrate that the quantitative criteria of formal causal inference are fulfilled. Particular concern applies to the relatively rapid rise of many parts of Europe into the higher dose cannabinoid exposure zone including by contamination of the food chain where major genotoxic outcomes become commonplace. As discussed there is some evidence that this may be occurring in parts of France and Germany already. Another major and related corollary is real concern for multigenerational damage to the epigenome. Evidence presented indicates that epigenomic explanations may well underlie many of the effects observed and the potential for transgenerational transmission of epigenomic effects is very real and of great concern. Such issues, together with the robust and causal relationships demonstrated in the present investigation indicate that community cannabinoid penetration should be carefully and closely constrained and restricted and appropriate protections must be put in place for protection of the food chain and the genomes and epigenomes of generations yet to come.

Supplementary Materials:
The following supporting information can be downloaded at: https:// www.mdpi.com/article/10.3390/ijerph192113769/s1, Figure S1: Geospatial Links; Table S1: Overall  Study Profile; Table S2, Daily Cannabis Use-Raw Data; Table S3: Daily Cannabis Use-Interpolated Data; Table S4: All regression slopes and results from bivariate analyses; Table S5: Variable Importance  Tables from Ranger random forest regression-Uronephrological anomalies; Table S6: Variable Importance Tables from Ranger random forest regression-Multicystic renal disease; Table S7: Variable Importance Tables from Ranger random forest regression-Bilateral renal agenesis; Table S8: Variable  Importance Tables from Ranger random forest regression-Hydronephrosis; Table S9: Variable Importance Tables from Ranger random forest regression-Congenital posterior urethral valve; Table S10: Inverse probability weighted panel regression results-Uronephrological anomalies; Table S11: Inverse probability weighted panel regression results-Multicystic renal disease; Table S12: Inverse probability weighted panel regression results-Bilateral renal agenesis; Table S13: Inverse probability weighted panel regression results-Hydronephrosis; Table S14: Inverse probability weighted panel regression results-Congenital posterior urethral valve.
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. Informed Consent Statement: Patient consent was waived due to all data being publicly available, anonymous and grouped. No identifiable data was accessed used, employed or published.

Institutional Review Board
Data Availability Statement: 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 has been made publicly available on the Mendeley Database Repository and can be accessed from these URL's: https://data.mendeley.com/datasets/c6psrbr34j/2 (accessed on 1 February 2022) and https://data.mendeley.com/datasets/vd6mt5r5jm/1 (accessed on 1 February 2022).