Supplementary Table ST1 provides an overview of study data. As shown, total CAR data were derived from 14 nations for 24 cardiovascular anomalies for the years 2010–2019. The total CAR includes stillborn and early termination of pregnancy for anomaly rates in addition to babies born normally alive and is a critical metric required to properly assess trends in CARs. The total sample size of the whole Eurocat database for these nations was 77 410 rates. Drug and cannabis use and exposure rates are also shown in this table including various compound metrics. Median household income is also shown.

Daily cannabis use data across the various jurisdictions of Europe were largely incomplete as indicated in Supplementary Table ST2 where only 59 raw data points are shown. These data were completed by linear interpolation, as shown in Supplementary Table ST3, so that 129 points were finally available for analysis.

Figures 1–3 show the rates of the various CAs against the use of the different substances. Eight CAs are shown in each figure which has been split into three parts to accommodate all 24 CAs. From these figures, it is apparent that tobacco exposure is for the most part negatively related to CAR. Regression line slopes for alcohol are mostly flat or negative but generally weak. Regression line slopes for amphetamine are mostly flat or somewhat positive. The lines of best fit for cocaine are generally positive, and some, such as CHD, double outlet right ventricle and mitral valve anomalies, are quite strongly so. The cannabis exposure metric chosen in this graph is daily cannabis use interpolated. For most anomalies, the regression line is flat or weakly positive. However, for some anomalies such as CHD, double outlet right ventricle and mitral valve anomalies, the regression line indicates strong and robust obviously positive associations.

Figures 4–6 chart the bivariate relationships of the different cannabis metrics with the spectrum of CARs. Some of the steepest positively sloped lines are for atrioventricular septal defect (AVSD), severe CHD, double outlet right ventricle, pulmonary valve atresia, transposition of the great vessels and hypoplasia of the right heart which are all seen as a function of cannabis resin THC concentration.

Figure 7 shows a graphical map of the CHD rate across Europe during this decade. Many European nations are darkly shaded, indicating relatively high rates of CHD. Figure 8 depicts the rates of severe CHD which appears to have a more varied pattern. Figures 9–12 illustrate the rates for ASD, VSD, PDA and Ebstein’s anomaly, respectively.

Patterns of daily cannabis use are illustrated in Fig. 13 which indicates increases in Spain, Netherlands, Belgium, France and Norway, reductions in Poland and little change in Bulgaria.

Figure 14 is a bivariate colorplane plot showing the bivariate relationship between CHD and the compound cannabis metric last month cannabis use × resin THC concentration × daily use interpolated. The areas covered by Spain and France are seen to obviously turn from tan and brown in 2010 to purple and crimson in 2019, indicating that they have moved up into the high cannabis metric—high CHD zone of the colorplane (shown in the key at the right hand side of the graph).

A similar change occurs across France and Spain for severe CHD as indicated in Fig. 15.

For the relationship between cannabis herb and CHD, the map is noted to have turned increasingly purple with time as nations more into the higher zone on both covariates (Fig. 16). Similar comments apply to the relationship between severe CHD and cannabis herb THC concentration, as shown in Fig. 17.

Figure 18 shows the bivariate relationship between severe CHD and cannabis resin THC concentration. Here, the map is noted to turn pink and purple across most countries across the decade.

Figure 19 shows the relationship between AVSD and cannabis resin THC concentration. The area for France changes shades of purple across the period. The areas covered by Germany and Norway appear to change in their blue-purple hues across this period.

When the relationship between transposition of the great vessels and cannabis resin THC concentration is considered, the map is also noted to have become more purple generally (Fig. 20).

Based on recent findings reported in major epidemiological reviews [22], it is possible to categorize European nations into increasing and high daily use and low or decreasing daily use groups. A panelled series of these defects across the decade are shown in Fig. 21. For most anomalies, there is a significant overlap between the areas covered by both sets of nations. However, for some anomalies, such as CHD, double outlet right ventricle, hypoplastic left heart syndrome, mitral valve anomalies and tricuspid valve stenosis or atresia, the rate in the nations with increasing daily use is higher than those which do not have this feature. Mixed-effects regression with anomaly as the random effect confirms that the time: daily use rate interaction is significantly higher in the countries with increasing cannabis use (β-Est. = 0.267, t = 2.756, P = 0.059; model AIC = 5317.51, LogLik. = −2652.75).

The slopes of 275 of these regression lines for each anomaly and each substance are provided in Supplementary Table ST4 using a purrr-broom simultaneous multigroup analytical workflow. From this list, the 69 models with positive and significant regression coefficients were selected out and are shown in Table 1. The table is ordered in descending order of minimum E-values (mEV). It is of interest therefore that the first 30 terms in this table all relate to cannabis exposure metrics. Of the 69 terms, 15 relate to cocaine, 5 to amphetamines, 1 to alcohol and 48 (69.6%) to various cannabis metrics. It is noted that the E-value estimates descend from 1.35 × 10¹⁹, and the mEVs decline from 4.78 × 10¹⁰.

Table 1 lists both P-values and E-values. By analogy with the genomic and epigenomic literature, it is therefore possible to present these data for graphical consideration as a ‘volcano plot’ which plots the negative logarithm of the P-value against the fold change in the anomaly rate quantified by the E-value. Figure 22 presents such a plot for the E-value estimate itself, and Fig. 23 presents a similar plot for the mEV. In these figures, strong signals are indicated for CHD and mitral valve anomalies and many other CAs as indicated.

Given these impressive results at bivariate analysis, it is of interest to investigate how all of these variables compare in multivariable regression techniques. Given that there are so many primary and compound cannabis metrics, it is not entirely clear which covariates are most suitable to use, given that the regression equations to be employed have a very finite number of degrees of freedom.

The issue of covariate selection was addressed formally using random Forrest regression in the ranger package in tandem with the variable importance package vip to create tables of relative variable importance. These tables are presented as Supplementary Tables S5–S16 for the 12 CAs in which we were most interested.

As the data are well suited to panel regression, this format was chosen to perform multiple regression analyses on these data. All panel regression models were inverse probability weighted which is an important step that transfers the analysis from merely an observational study into a pseudo-randomized analytical and interpretive context.

Supplementary Table S17 shows the result of panel regression of the CHD dataset against the group of covariates indicated in the table. Additive, interactive and models that lagged temporally to 1 and 2 years are presented. In each case, several terms including cannabis have positive regression coefficients and are highly statistically significant.

This same pattern of results is continued across all of the selected CAs, namely severe CHD, ASD, VSD, AVSDs, patient ductus arteriosus, tetralogy of Fallot, vascular disruption, double outlet right ventricle, transposition of the great arteries, hypoplastic right heart syndrome and mitral valve anomalies, as shown in Supplementary Tables S18–S28. In all cases, all multivariable panel models are inverse probability weighted to allow causal relationships to be formally assessed.

Having made the important above observations, the next issue was whether these relationships would persist when formal space-time modelling was undertaken in a fashion which allows random effects, serial correlation, spatial correlation and spatial autocorrelation in the error structure to be formally considered. Geospatial links were derived, edited and finalized, as indicated in Supplementary Fig. S1, and these formed the basis of the sparse spatial weights matrix which was then employed in the R package spml to conduct formal spatiotemporal regression.

Tables 2–13 present the results of final reduced models for each of the 12 anomalies cited in the preceding paragraphs immediately above. In virtually all models terms including cannabis metrics persist in the final models, have positive regression coefficients and are statistically significant usually at high levels.

Tables 14 and 15 present the E-values to emerge from these positive and significant cannabis-related terms in the panel and spatial regression model sets, respectively.

These lists are then combined into the lists of 149 E-values derived from both sets of models, as shown in Table 16. It was found that 138/149 (92.6%) of the E-value estimates were >9 and so in the high zone, and all 149 exceeded the cut-off for causality at 1.25 [67]. 112/149 (75.2%) mEVs exceeded 9 and so fell in the high zone [68], and 147/149 mEVs (98.7%) exceeded the threshold for causality at 1.25 [67].

Table 17 lists these E-values by CA. Table 18 summarizes these data by CA listed in descending order of median mEV. This table is of interest for several reasons and is considered further in the Discussion section below.

These multivariable data are also similarly suited to illustration by Volcano plots as was done earlier for the bivariate results.

The multivariable data for the E-value estimates are shown in Supplementary Fig. S2. In this plot, the presence of the tetralogy of Fallot far to the right forces all the points to the left side of the graph and makes the remainder of the graph difficult to interpret clearly. A similar problem occurs when the mEVs are considered, as shown in Supplementary Fig. S3.