Main

A global ban during the first half of the 20th century saw hemp (Cannabis sativa) production relegated to the sidelines of industry and agriculture, despite its potential as a cheap, ecological and versatile crop1. The recent revision of cannabis regulations has triggered a fast-paced growth in the industrial hemp sector which has brought to the market many novel hemp-derived products1,2. This proliferation raises consumer safety issues because hemp contains cannabinoids, a class of substances that interact with the animal and human endocannabinoid system. Some of them—especially Δ9-tetrahydrocannabinol (Δ9-THC)—exert a psychoactive effect and others are only pharmacologically active, such as cannabidiol (CBD)3. Cultivation of industrial hemp with a maximum Δ9-THC content of 0.2% is allowed in the European Union4. In 2021, an increase to 0.3% was initiated in the European Union5.

A transfer of cannabinoids into foods of animal origin is conceivable when by-products of hemp production and the whole plant are used as feedstuffs6. So far, despite a few case reports, experimental data regarding the transfer of Δ9-THC from feed into the milk of cows are scarce and analytical techniques often failed to differentiate between psychoactive Δ9-THC and its non-psychoactive precursor Δ9-tetrahydrocannabinolic acid (Δ9-THCA)6. Given that Δ9-THCA can be up to nine times more abundant than Δ9-THC in hemp7, this differentiation is crucial to understand transfer processes and subsequently perform risk analysis.

In this article, we focus on the effects of industrial hemp silage feeding in lactating dairy cows with the aim of quantifying the transfer of cannabinoids into milk and determining possible effects on animal health and risks for consumer health. We collected and analysed milk, blood plasma and faeces, measured physiological parameters and observed animal behaviour. We employ a liquid chromatography–tandem mass spectrometry-based analytical technique that ensures differentiation between Δ9-THC and Δ9-THCA in various matrices, and enables quantification of the cannabinoids Δ8-tetrahydrocannabinol (Δ8-THC), Δ9-tetrahydrocannabivarin (Δ9-THCV), CBD, cannabinol (CBN), cannabidivarin (CBDV) and the two Δ9-THC metabolites 11-hydroxy-Δ9-THC (11-OH-THC) and 11-nor-9-carboxy-Δ9-THC (THC-COOH). The resulting data were used to develop a predictive toxicokinetic model, which can be used to simulate other exposure scenarios and to assess the transfer of different cannabinoids into cows’ milk when using industrial hemp as a dietary supplement for dairy cows.

Results and discussion

We conducted a feeding experiment with lactating Holstein Friesian dairy cows, where corn silage in the diet was first partially replaced with hemp silage made from whole plant hemp (very low cannabinoid concentration, hemp silage A) during an adaptation period followed by feeding hemp silage made from leaves, flowers and seeds only (higher cannabinoid concentration, hemp silage E) at two different supplementation levels (group L, low hemp, 0.84 kg dry matter per cow per day; group H, high hemp, 1.68 kg dry matter per cow per day) during an exposure period, and a subsequent hemp-free depuration period (Fig. 1a). Partial replacement of corn silage with hemp silage of different nutritional composition also partially changed the total nutritional composition of the diets; however, the major difference was in total cannabinoid concentration (Table 1 and Supplementary Table 1).

Fig. 1: Study design and results of the feeding experiment.
figure 1

a, Study design and sampling strategy. Grey arrows represent sampling days, striped lines in arrows indicate no sampling for this parameter on that day. aAll animals per group (n = 4). bTwo animals per group. b–e, Feed intake (b), milk yield (c), respiratory rate (d) and heart rate (e) during the feeding experiment. Violin plots represent the density distribution; black error bars are the standard deviations of the means (black dot); red dots represent the median; blue dots (jittered) are raw data. Means followed by a common letter are not significantly different by a Tukey test. P values (with t values and 95% confidence intervals for standardized effect size measures) are from multiple comparisons according to Tukey post hoc tests and are listed in Supplementary Table 2 and Supplementary Fig. 1. To show differences more clearly, standardized effect size measures of animal health parameters were plotted between time periods (pairwise comparisons) with 95% confidence intervals. f, Pictures of cows during exposure period: (1) prolapsed and reddened nictitating membrane; (2) increased salivation with ponding; (3) pronounced tongue play; (4) nasal secretion.

Full size image
Table 1 Cannabinoid concentration (mg per kg dry matter) of the different hemp silages
Full size table

Feeding industrial hemp can affect animal health

The study shows that feeding up to 0.92 kg per cow per day (DM, dry matter) of industrial hemp silage A with very low cannabinoid concentration to dairy cows during the adaptation period had no effect on physiological parameters and health (Fig. 1b–e; see Supplementary Table 2 and Supplementary Fig. 1 for statistical data).

In contrast, feeding of cannabinoid-rich industrial hemp silage E during the exposure period had a significant effect for both group L and group H. A Kruskal–Wallis test showed that feed intake (χ2 = 28.9, P < 0.001, d.f. = 5) and consequently milk yield (χ2 = 44.5, P < 0.001, d.f. = 5) decreased significantly from the second day of exposure period in both experimental groups (Fig. 1b,c). The average ingested doses of ∆9-THC and CBD were 1.6 ± 0.3 and 10.7 ± 1.9 mg per kg body weight for group L, and 3.1 ± 0.7 and 20.4 ± 4.4 mg per kg body weight for group H, whereas the dose of the other cannabinoids was ≤1.1 mg per kg body weight (Supplementary Table 3). Similarly, the respiratory (χ2 = 50.2, P < 0.001, d.f. = 5) and heart rate (χ2 = 77.4, P < 0.001, d.f. = 5) decreased significantly within hours in both groups due to the feeding of hemp silage E (Fig. 1d,e) and in individual animals fell below physiological reference values, allowing classification as either bradypnea or bradycardia. Concomitantly, changes in animal behaviour and appearance were evident (Fig. 1f), such as pronounced tongue play, increased yawning, salivation, nasal secretion formation, prolapse and reddening of the nictitating membrane, and somnolent appearance. Some animals from group H displayed careful, occasionally unsteady gait, unusually long standing and abnormal posture. All changes observed disappeared within two days of discontinuing cannabinoid feeding. Milk constituents (fat, protein, lactose, dry matter, somatic cell count and urea), body temperature and body weight were unaffected in all periods (Supplementary Tables 4 and 5).

For humans, the lowest-observed-adverse-effect-level for ∆9-THC is 0.036 mg per kg body weight6. On average, the cows ingested up to 86 times more ∆9-THC during the exposure period, probably explaining the health effects they exhibited. Another study in ruminants reported no health effects, although this involved feeding only 35 g of industrial hemp, a much lower exposure than the present study8. In animals, only a few studies have reported a reduction in feed consumption due to administration of ∆9-THC and CBD, whereas other studies showed no significant decrease or even an increase in feed consumption9,10,11,12,13. Other factors in the silages, such as fermentation acids or lignin (Supplementary Table 6), may also have led to the drop in feed intake, but were not observed when feeding industrial hemp silage A during the adaptation period. In contrast, hemp silage E was markedly richer in fat, which could also have led to a reduction in feed intake14. Studies on the suitability of industrial hemp components, other than seeds and oil, as animal feed are rare15, complicating the interpretation of our results. The reduction in daily milk yield is presumably caused by the decrease in feed intake. A decrease in water consumption (not measured) due to cannabinoids in the feed may have promoted the decrease16,17,18. Bradypnea and bradycardia are rare symptoms in cows that only occur in the course of serious illnesses or can be pharmacologically induced. In animal studies, administration of ∆9-THC (but not CBD) influenced respiratory and heart rate and caused drowsiness, slow movements, ataxia and salivation12,13,19,20,21. A reddening of the conjunctiva after ∆9-THC intake is known in humans22,23.

Although it is clear that the observed effects of industrial hemp silage feeding on animal health were mainly caused by the cannabinoids, it cannot be clearly defined which cannabinoid was responsible. Due to its high concentration in the cannabinoid-rich silage, ∆9-THC is the most likely cause, but combination effects may also play a role16,17,24,25. Finally, it should be considered that hemp plants produce other phytochemicals (for example, terpenes, flavonoids) that may have an impact on the changes observed26,27. Similarly, an effect due to minor changes in nutrient concentration of the diets by replacement of corn silage with hemp silage (Supplementary Table 1) could not be completely ruled out.

Our study shows that the feeding of industrial hemp silage to dairy cows, even in small amounts, is associated with health consequences. These seem to be dependent on the cannabinoid concentration of the silage, which is influenced by, among other factors, the variety, the parts of the plant from which the silage is derived and the time of harvest7,28. Due to the manifold parameters influencing the cannabinoid concentration of hemp and derived feed, the innocuousness of hemp cannot be reliably assessed without prior cannabinoid analysis.

Cannabinoid transfer from feed to cow’s milk

Feeding hemp silage resulted in measurable levels of ∆9-THC, ∆9-THCA, ∆9-THCV, CBD, CBN and CBDV in cow’s milk at the end of the adaptation period and during the exposure period. Concentrations of up to 316 µg ∆9-THC and 1,174 µg CBD were detected per kg milk (Fig. 2). Maximum values for other cannabinoids in milk were 1.9 (∆9-THCA), 8.0 (∆9-THCV), 2.5 (CBN) and 10.1 µg per kg body weight (CBDV). On the last day of the depuration period, ∆9-THC (group L: 1.4 ± 0.4 µg per kg; group H: 5.0 ± 0.6 µg per kg) and CBD (group L: 7.0 ± 1.9 µg per kg; group H: 16.2 ± 2.6 µg per kg) were still detectable in milk. The other cannabinoids analysed were undetectable in milk during all periods (limit of detection (LOD): 0.1 µg per kg milk for ∆8-THC, 0.5 µg per kg milk for 11-OH-∆9-THC and THC-COOH). At the end of the exposure period, the levels in milk were 6–26 times higher for ∆9-THC, 3–5 times higher for ∆9-THCV and 11–32 times higher for CBD than the corresponding blood plasma levels, pointing to the accumulation potential of the substances in cow’s milk. In contrast, ∆9-THCA did not accumulate in milk and CBDV showed no clear trend (Supplementary Table 7). Since CBN was undetectable in plasma, no ratio could be calculated. Urine was collected during the experiment, but due to analytical challenges, that is, inconsistent results after glucuronidase treatment, the data could not be evaluated.

Fig. 2: Cannabinoid concentration in milk.
figure 2

Mean ± s.d. ∆9-THC and CBD levels. Days 1–6 and 8–13: no milk samples analysed. *Data interpolated. During the experimental periods, milk samples were taken on days 7, 14–24, 26 and 28 from both groups L and H (n = 4 per group) for cannabinoid analysis. To compile a sample for analysis, the evening and the following morning milk of each cow were mixed according to their proportion of daily production. The resulting composite samples were stored at −20 °C until HPLC–MS/MS analysis.

Full size image

In contrast to an earlier study in Holstein calves in which only ∆9-THCA and traces of CBD could be detected in the plasma after hemp feeding8, we were able to detect ∆9-THC, ∆9-THCV, CBD and CBDV, presumably due to the higher dose. Studies on the transfer of cannabinoids from feed into milk are rare and, up to now, only ∆9-THC and THC-COOH have been found in ruminant milk after oral ∆9-THC-exposure29,30. In humans, the transfer of ∆9-THC and CBD into breastmilk after marijuana use has been demonstrated, whereas 11-OH-THC and THC-COOH have not always been found31,32,33. The wide range of reported milk/plasma ratios for ∆9-THC is consistent with an earlier study31.

Toxicokinetic model

We set up two-compartment models, as shown in Fig. 3a, with either two (model A) or one (model B) excretion outputs from the central compartment. Model A was used for Δ9-THC and CBD, for which enough data were available, and the simplified model B for the other cannabinoids with fewer data. The model reproduces milk data with a Pearson coefficient of determination r2 of 0.92 (Supplementary Fig. 2). The model can approximate Δ9-THC and CBD plasma and faeces levels with a Pearson r2 of 0.93 and 0.90, respectively, even though few data points were available. Figure 3b,e shows the relative fate of Δ9-THC and CBD, respectively. Only a small fraction ends up in milk, with most of the substance (77% Δ9-THC and 64% CBD) being ‘eliminated’ in processes such as putative metabolism, biochemical transformations in the gastrointestinal tract and urinary excretion. This hints at a high bioavailability for the compounds, with 20–35% being excreted into faeces. Figure 3c,d (Δ9-THC) and Fig. 3e,f (CBD) show the milk predictive model, obtained using the average of the single model parameters estimates, plotted against the data of all animals in the experiment. The simulation was performed to simulate the same intake scenarios as in the experiment. The predictive models for Δ9-THC and CBD can predict the milk levels in each scenario (L/H) with an average underestimation between 12% and 26% for the exposure period. Nevertheless the model predictions are in line with the observed data range (Pearson r2 between 0.59 and 0.88), such that they appear sufficient for risk assessment. Model B was successfully used to predict cannabinoid milk levels and to provide an estimate for the transfer rates at steady state. The figures showing the results for Δ9-THCA, Δ9-THCV, CBN and CBDV are included in Supplementary Fig. 3.

Fig. 3: Overview of the modelling results.
figure 3

a, Model structure schema. Model A (used for Δ9-THC and CBD) has three output routes: faeces, elimination and milk. Model B (used for the other cannabinoids) is a simplification of model A, with only two outputs: total elimination and milk. b, Modelled fractions of Δ9-THC output. Most Δ9-THC is eliminated, some is excreted with faeces with only a small amount found in milk. c,d, Predictive model of Δ9-THC excretion in milk for group L (c) and group H (d). The plots show the Δ9-THC amount (mg d−1) excreted in milk against the experimental data of the cows in the group, in a simulation of the experiment with low (c) and high (d) dose. e, Modelled fractions of CBD output. Compared to Δ9-THC, CBD is excreted in higher amounts in the faeces, hinting to lower bioavailability. f,g, Predictive model of CBD excretion in milk for group L (f) and group H (g). The plots show the CBD amount (mg d−1) excreted in milk against the experimental data of the cows in the group, in a simulation of the experiment with low (f) and high (g) dose.

Full size image

Both the data and the models show that cannabinoids display biphasic elimination: there is a rapid decline in milk levels following discontinuation of cannabinoid feeding and further decline over a longer period.

Cannabinoids are metabolized in the liver of mammals, where they undergo biotransformation. Our model lumps these processes into ‘elimination’ of the central compartment (model A) and total elimination (model B), without details. Having information on these processes would allow for a more refined model and unify the models for individual cannabinoids into a single one including interconversion. Such processes could be quantified using both in vivo and in vitro techniques. Finally, it should be noted that there is still no information on the fate of cannabinoids in the gastrointestinal tract of the cow prior to absorption.

Predictive model and transfer rates at steady state

The predictive model for each cannabinoid was used to estimate transfer rates at steady state for each substance. The transfer rate is defined as the daily amount of cannabinoids in milk divided by the amount in feed (µg µg−1). For ∆9-THC at steady state, we estimate a milk transfer rate of 0.20%, or around 33–100% higher than reported in previous experiments6,29. No data on cow’s milk transfer rates for the other cannabinoids analysed was available from previous literature. We present here the milk transfer rates for ∆9-THCA (0.015%), Δ9-THCV (0.56%), CBD (0.11%), CBN (0.043%) and CBDV (0.0080%) in cows (Supplementary Table 8). The transfer rate of Δ9-THCA appears much lower than its decarboxylated counterpart Δ9-THC, which could mean lower bioavailability and/or that some of the Δ9-THCA does not reach the intestinal tract of the animal, but instead undergoes chemical reactions in the rumen (decarboxylation, bacterial metabolization). The latter putative processes are theoretically possible for each of the studied cannabinoids, as they have similar chemical structures and can potentially interconvert34. Although Δ9-THCV has the highest transfer rate, the low concentrations found in hemp silages make it unlikely to pose a risk to consumers. Overestimation of its transfer rate through putative rumen processes cannot be ruled out. The transfer rate of CBD is in the same range as that of Δ9-THC. Since the CBD concentration in the silage was much higher, as usual for industrial hemp, this resulted in markedly higher milk levels. Thus, CBD may be as relevant as Δ9-THC when assessing consumer risks. CBN and CBDV showed the lowest transfer rates and very low levels in milk, confirming that Δ9-THC and CBD are the main compounds of interest for animal health and consumer safety.

Human exposure to ∆9-THC via milk may exceed the acute reference dose

The acute reference dose (ARfD) is an estimate of the amount of a substance in food and/or drinking water that can be ingested in a period of up to 24 h without appreciable health risk to the consumer35. To obtain a rough estimate on the potential risk from exposure to ∆9-THC via milk containing levels as measured in the present study, a cursory toxicological assessment was conducted. The EFSA Rapid Assessment of Contaminant Exposure (RACE) software tool35 was used to estimate the exposure to ∆9-THC resulting from consumption of milk and dairy products, considering the measured ∆9-THC levels, and was then compared to the acute reference dose (ARfD) of 1 µg Δ9-THC per kg body weight derived by EFSA6.

For mean and maximum ∆9-THC concentrations in milk from group L, average consumption of milk and dairy products leads to the ARfD being exceeded by up to 14-fold for consumers aged <18 yr, whereas high consumption (95th percentile) exceeds the ARfD up to 57-fold in all population groups. Based on the mean and maximum ∆9-THC levels in milk from group H, the ARfD is exceeded up to 120-fold for all population groups (Table 2). Remarkably, even consumption of milk and dairy products from the adaptation period (lowest Δ9-THC in feed) leads to the ARfD being exceeded by 1.5-fold in infant high consumers (95th percentile). Additionally, it should be taken into account that Δ9-THCA present in milk could be converted into Δ9-THC during processing of the milk (for example, ultrahigh-temperature processing)6.

Table 2 Comparison of the exposure to ∆9-THC and the ARfD for different population groups
Full size table

Thus, the exposure scenario based on the data obtained from the transfer study presented herein shows that exposure to ∆9-THC via milk and dairy products could lead to the ARfD being exceeded in some consumer groups.

Comparison with the Community method for determination of ∆9-THC

We compared our analytical approach to the Community method for the quantitative determination of Δ9-THC content in hemp varieties, which is laid down in Annex III of Regulation (EU) Number 639/2014. The approaches differ substantially in extraction method, extraction time, extraction solvent and analytical determination. The Community method uses hexane containing squalane as an internal standard in an ultrasonic bath for extraction, and without further clean-up steps, the analytical separation and determination is performed isothermally by gas chromatography with flame ionization detection (GC-FID). This method has drawbacks: the extraction of ∆9-THC and ∆9-THCA is not complete, and their separate determination is not possible. Additionally, ∆9-THCA may be converted in the GC injector into ∆9-THC. Moreover, coelution of Δ9-THC with other cannabinoids or coextractives of hemp due to the limited clean-up of the crude sample extract and the isothermal GC conditions cannot be excluded. These drawbacks are overcome by our approach. The use of a Soxhlet apparatus enables an almost exhaustive extraction by repeatedly adding freshly dripping solvent, avoiding saturation. In preliminary tests, acetone was found to be more suitable than n-hexane for extraction of ∆9-THC and ∆9-THCA. By using the Community method, the total THC (the sum of ∆9-THC and ∆9-THCA) was 21 and 186 mg per kg DM, respectively for the two industrial hemp silages. In contrast, our new method yielded 65 and 1,317 mg per kg DM, respectively. Analytical determination using high-performance liquid chromatography coupled to tandem mass spectrometry (HPLC–MS/MS) with isotope-labelled internal standards enables the unequivocal separation and identification not only of Δ9-THC and its precursor but also of further cannabinoids at low limits of quantification. Therefore, we question the sensitivity of the Community method when used to determine the actual total THC concentration in hemp.

Conclusion

Our study shows that feeding cannabinoid-rich industrial hemp silage made from leaves, flowers and seeds leads to a decrease in feed intake and milk yield in dairy cows. Heart rate, respiratory rate and animal behaviour were also negatively affected. Feeding a low-cannabinoid industrial hemp silage made from the whole hemp plant showed no effects on cow health and performance. However, apart from dietary cannabinoids, other factors due to partial replacement of corn silage with industrial hemp silage in this experiment (that is, change in nutrient composition, other secondary plant compounds) could not be completely ruled out. The HPLC–MS/MS-based analysis method we developed enables a precise determination of the cannabinoid concentration, and differentiates between psychoactive Δ9-THC and non-psychoactive Δ9-THCA. The toxicokinetic modelling has shown that the transfer rates of the examined cannabinoids from feed to milk were less than 1%. Nevertheless, due to the high feed intake, cow’s milk reached substantial levels of Δ9-THC such that the exposure might exceed ARfD in some population groups in our exposure scenario based on the transfer properties presented here. For other cannabinoids, in particular for CBD, which was present in high amounts in industrial hemp (and thus also in cow’s milk after feeding), the data are currently insufficient, thereby preventing an assessment of possible health risks.

Methods

Feeding experiment

A feeding experiment was conducted with ten lactating Holstein Friesian dairy cows at the German Federal Institute for Risk Assessment (BfR). The experiment was performed in accordance with relevant guidelines and regulations and was approved by the Regional Office for Health and Social Affairs Berlin (LAGeSo, approval number G 0239/17). For detailed data on the animals at the beginning of the experiment, see Supplementary Table 9.

The animals were divided in two groups (L, low hemp; H, high hemp) and housed in cubicles. A partially mixed ration (PMR) was administered in feeding troughs with transponder-controlled scales allowing recording of individual feed intake. The experiment was divided into four periods: control (7 d), adaptation (7 d), exposure (6 d) and depuration (8 d). In the control period, both groups received a hemp-free PMR, containing 7.56 kg corn silage, 1.03 kg wheat straw, 0.87 kg hay, 1.17 kg fodder beet, 3.60 kg rapeseed meal and 0.88 kg of a supplementary feed rich in vitamins and minerals per cow and day (all data in DM). The PMR covered a milk yield of 15 litres per cow per day. In addition, concentrate was provided through transponder-controlled automated feeders based on individual milk yield (0.5 kg per additional litre of milk yield; maximum 5 kg of concentrate per day). Water was offered ad libitum. In the adaptation period, part of the corn silage was gradually replaced with industrial hemp silage A (0.31–0.92 kg DM per cow), containing cannabinoids at very low concentrations. In the exposure period, groups L and H had the corn silage replaced proportionately with 0.84 and 1.68 kg DM per cow, respectively, of industrial hemp silage E with higher cannabinoid concentration. Different dosages were chosen to determine dose-dependent effects on the transfer of cannabinoids into milk, cow’s performance and health status. In the depuration period, both groups were again fed with hemp-free diets. Energy content of the PMR was 6.4–6.6 MJ net energy for lactation per kg DM, depending on the experimental period. For more information on the composition of the diets, see Supplementary Table 1. Throughout the experiment, animal health parameters and performance (feed intake, milk yield, respiratory rate, heart rate and body temperature) were recorded daily, as were behavioural observations. For this the animals were observed for conspicuous behaviour compared to the control period. The focus was on the parameters ‘muscle tremors’, ‘hypersalivation’ and ‘disturbed movement coordination’ described in the literature. Since only an initial recording of possible effects was intended, the observations were not standardized using an ethogram. The respiratory rate of all animals was determined by counting the respiratory movements of the flank region and expressed as breaths per minute. The heart rate was determined by auscultating the chest wall with a stethoscope and expressed as beats per minute. The body temperature was measured rectally. Animal weight was measured weekly using floor scales. Milk samples were collected for each cow twice a day, while blood and faeces were taken only from two cows per group and on few, selected days.

Two different industrial hemp varieties were used for the experiment. The industrial hemp variety ‘Ivory’ was selected as a representative of particularly low Δ9-THC varieties for the adaptation period of the experiment; the industrial hemp variety ‘Finola’ was selected as a representative of varieties with rather higher Δ9-THC concentration (but still below 0.2% total Δ9-THC) for the exposure period. Both varieties were cultivated at the Leibniz Institute for Agricultural Engineering and Bioeconomy Potsdam (ATB), Germany, in spring 2016. The harvest of the hemp took place at the end of August 2016, at which time the ‘Ivory’ variety was in the seed formation stage. The ‘Finola’ variety already bore the first fully mature seeds. The ‘Ivory’ variety was cut with a harvester (JF FC 800) and the entire plant was chopped to produce whole plant silage (silage A). For the production of a leaf-flower-seed silage (silage E), only leaves, flowers and seeds of the ‘Finola’ variety were harvested. For this purpose, a home-made harvester, constructed by the ATB, was used, which removed the desired plant parts from the stem by means of whips. Immediately after harvest, the hemp was moved into round 510 litre rain barrels, compacted in layers by body weight (tamping) and then sealed airtight. A water-filled hose in the barrel lid allowed gas to escape from the ensiled material without allowing ambient air to flow in. No silage additives were added. In October 2017, the industrial hemp silages were opened, aliquoted in 120 litre bags and stored at −20 °C for the experiment. The cannabinoid concentration of the industrial hemp silages can be found in Table 1.

Sampling

Milk samples for cannabinoid analysis were taken at days 7, 14–24, 26 and 28. The cows were milked in a double tandem milking parlour until the milk flow stopped and the milk was collected in milking pails. Before sampling, milk was stirred to achieve an even distribution of the lipophilic cannabinoids. The resulting milk samples were refrigerated (8 °C) for a maximum of 18 h until mixing. After reheating (37 °C, water bath), samples of evening and the following morning milk of each cow were mixed according to their proportion of daily production. The resulting composite samples were stored in lightproof high-density polyethylene bottles and were frozen at −20 °C until analysis. Milk samples for the determination of milk constituents were taken at the end of each period and were refrigerated at 8 °C until analysis.

Blood plasma and faeces were collected from two animals per group at days 7, 14, 20 and 22. Blood was obtained by puncturing the jugular vein and collected in coagulant tubes with potassium oxalate monohydrate. After cooling for 30 min, blood samples were centrifuged at 1,876g for 15 min and the plasma extracted and stored for analysis. Blood plasma was used instead of whole blood because of ease of handling and because plasma has been shown to contain a large proportion of blood cannabinoids36,37,38. Faeces was collected during spontaneous defecation, when possible. Otherwise, faeces was collected from the rectum. Plasma and faeces samples for the determination of cannabinoids were stored in lightproof high-density polyethylene bottles and were frozen at −20 °C until analysis. Samples of the hemp silage were taken prior to the experiment and were stored in freezer bags at −20 °C until analysis.

Cannabinoid analysis

The cannabinoid analysis was performed at the Chemical and Veterinary Analytical Institute Münsterland-Emscher-Lippe (CVUA-MEL), according to an approach described in ref. 39. The method comprises the following native cannabinoids: Δ9-THC, Δ9-THCA, Δ8-THC, Δ9-THCV, CBD, CBN, CBDV, and the two Δ9-THC-metabolites 11-OH-THC and THC-COOH. For calculation of the recoveries and quantification, the following deuterated internal standards were used: Δ9-THC-D3, CBD-D3, CBN-D3, 11-OH-Δ9-THC-D3, 11-nor-9-carboxy-Δ9-THC-D3 and 11-nor-9-carboxy-Δ9-THC-D9. Except for Δ9-THCA, which was purchased from Fluka, all other compounds were acquired from Cerilliant.

At the beginning of the analysis, the samples were fortified with a mixture of deuterated internal standards. A further deuterated standard was added at the end of the clean-up prior to the analytical determination. The extraction and the subsequent clean-up steps were adapted to the various matrices and differ to some extent. The final analytical separation, identification and quantification by HPLC–MS/MS were the same for all matrices.

Milk extraction and clean-up

The extraction of cannabinoids from milk was performed with acetone in a polypropylene tube by shaking on an orbital shaker. Precipitated proteins and sugars were removed after centrifugation. A subsequent liquid/liquid extraction with hexane in the presence of a saturated NaCl solution followed. After phase separation, the organic phase was subsequently evaporated to dryness under a gentle stream of nitrogen. Following reconstitution with cyclohexane/ethyl acetate (1 + 1), application of gel permeation chromatography (GPC) on Bio-Beads S-X3 served as the final clean-up step. After GPC clean-up, ethylenglycol was added to the extract which after careful evaporation was reconstituted in acetonitrile/water (50/50) for HPLC–MS/MS analysis.

Plasma extraction and clean-up

Plasma samples were subjected to an enzymatic cleavage with β-glucuronidase in phosphate buffer for 3 h at 50 °C. After cooling down, acetic acid and acetone were added to the solution for the extraction of the cannabinoids, which was performed by shaking on an orbital shaker. Plasma samples without cleavage were extracted with acetic acid and acetone by shaking on an orbital shaker. After centrifugation and removal of precipitated coextractives, a liquid/liquid extraction with hexane in the presence of a saturated Na2SO4 solution was performed. Ethylenglycol was added to the organic phase, which was carefully evaporated to dryness and subsequently reconstituted in acetonitrile/water (50/50) for HPLC–MS/MS analysis. Due to unexplained losses of Δ9-THC and other cannabinoids during β-glucuronidase treatment, we decided to use the analysis values without enzymatic cleavage for evaluation and modelling.

Feed and faeces extraction and clean-up

Feed and faeces samples were mixed with Na2SO4, transferred into an extraction thimble and extracted with acetone for 6–8 h in a Soxhlet apparatus. The extract was reduced and, depending on the concentrations of the cannabinoids, an aliquot of the concentrated extract was evaporated under a gentle stream of nitrogen to dryness. After reconstitution in cyclohexane/ethyl acetate (1 + 1), the extract was further purified on GPC as described for milk.

HPLC–MS/MS analysis

HPLC–MS/MS analysis was performed on an Agilent 6460A triple-quadrupole system equipped with an electrospray source operated in positive ionization mode. For separation of the analytes, a ZorBAX SB-C18 column (2.1 × 10 mm, 1.8 µm) at a column temperature of 30 °C was used. A flow rate of 0.2 ml min−1 and a gradient comprising of acetonitrile/0.1% formic acid (eluent A) and water/0.1% formic acid (eluent B) were applied.

For each cannabinoid and metabolite the collision energies in MS/MS mode were optimized to gain at least two intensive transitions. While the first transition was used for quantification based on a seven-point calibration curve, the second transition served as a qualifier; that is, the ratio of the two transitions had to comply with the corresponding reference standard as a second criterion besides the complying retention time for unequivocal identification of the respective analyte.

Depending on the cannabinoid and metabolite, the limits of detection (LOD) were determined as 0.01-0.4 ng ml−1 for milk, 0.1–0.5 ng ml−1 for plasma and 0.2–2.0 ng per g DM for feed and faeces. The corresponding limits of quantification (LOQ) were a factor of three higher.

The validation of the accredited method was performed by spiking experiments of cannabinoid-free matrices at a high and a low level, respectively, which generally revealed recoveries of >85% and coefficients of variation of <15%.

Analysis of milk constituents

The fat, protein, lactose and DM content of milk were analysed using accredited methods. Urea was measured using a Reflotron; somatic cell count was performed using a DeLaval DCC. All analyses were conducted at the BfR. The daily values of each parameter were calculated from the received values of the evening and the following morning milk samples.

Statistics

One animal from group L was excluded from the experiment at the beginning of the adaptation period due to a disease not related to the experiment (lameness). One animal from group H completed the experiment and was subsequently excluded from the statistical calculations due to a metabolic disease (subclinical ketosis) developed during the experiment. The statistical evaluation was conducted with the data from the remaining eight animals. Since our data did not follow a normal distribution, non-parametric tests were used.

Animal health parameters and performance—feed intake, milk yield, respiratory rate, heart rate, body temperature, body weight and milk constituents—were compared between six groups (1, groups L + H: control; 2, groups L + H: adaptation; 3, group L: exposure; 4, group H: exposure; 5, group L: depuration; 6, group H: depuration) applying the non-parametric Kruskal–Wallis test with the function kruskal.test() from the stats package. Post hoc analyses were performed with the non-parametric Tukey test (pairwise comparisons) with the glht function from the R package multcomp40 to find the groups that differ significantly after global the Kruskal–Wallis test rejected the null hypothesis.

To show the distribution of animal health parameters and performance in groups and periods, a violin plot combines the kernel density curve with an error bar41. The error shows the mean and the standard deviation. Different letters on plots indicate statistical significance in pairwise comparisons performed with Tukey post hoc tests. Letters shared in common between six groups (15 pairwise tests) indicate no significant difference. Statistical evaluations and graphs were performed with the statistical software R v.4.0.242. All statistical tests were two-sided.

Modelling

Data processing

As experimental data were available in concentrations and we opted for modelling cannabinoid amounts, we calculated cannabinoid amounts from the experimental data. To do so, we used both experimental data (that is, daily milk yield, body weight) and equations from veterinary literature to infer the missing information. To obtain the daily cannabinoid amount in milk Amilk,daily, the daily milk yield is multiplied with the cannabinoid concentration on each day. We assume each cow to have 65 ml blood per kg body weight (physiological range, 55–75 ml per kg body weight43) and a haematocrit of 33.8% (physiological range, 19.5–42.5%44), so that each cow has 43 ml plasma per kg body weight. The daily dry faecal output Mfaeces,daily was estimated using the mass balance, where

$$M_{\mathrm{faeces,daily}} = \frac{{1 - dDM}}{{0.10}} \times M_{\mathrm{feed}}$$
(1)

with digestible dry matter dDM = 79.36 – 2.85NL, where NL is the nutrition level of the cow, calculated from the cow’s daily feed uptake45. For concentration values under the LOD or LOQ, we set them to the respective LOD/LOQ value (upper bound estimation). For intake data at time t, I(t), we multiplied the feed intake amount with the cannabinoid concentration. Since there were problems with the urine analysis, we could not use the urinary concentration data.

Mathematical models

After analysing different options, we opted for a two-compartment model (Fig. 3a), which provides an accurate mathematical representation of the observed data, compatible with the limitations of the available dataset. We used two different versions of it, model A for Δ9-THC and CBD for which some plasma concentration data are available, and model B for the remaining cannabinoids (Δ9-THCA, Δ9-THCV, CBN and CBDV), for which plasma levels were not consistently detectable.

Model A has two compartments, GI (gastrointestinal) tract and plasma. The GI tract compartment takes the daily amount of cannabinoid ingested by the cow as input I(t), transfers the absorbed cannabinoid amount into the plasma and outputs the amount of cannabinoid excreted with faeces out of the system. The parameters regulating this process are Fa (fraction absorbed) and kp, the latter regulating the rate of passage out of the GI tract, which contains an AGItract amount of cannabinoids. We fixed kp at 1.3 d−1, a value tantamount to a GI tract passage rate of 5% h−1, in accordance with animal physiology literature46. The plasma compartment, which takes its input from the GI tract, outputs either to milk (Amilk, mainly controlled by parameter km) or to elimination (through parameter ke). The latter is mainly controlling the total amount of cannabinoids eliminated Aelim (lumping putative urinary excretion, metabolization and other possible physiological processes). The mathematical representation of the model is a system of ordinary differential equations that assumes zero- or first-order kinetics for all of the modelled processes, so that

$$\frac{{{\mathrm{d}}A_{\mathrm{milk}}}}{{{\mathrm{d}}t}} = A_{\mathrm{plasma}} \times k_{\mathrm{m}}$$
(2)
$$\frac{{{\mathrm{d}}A_{\mathrm{plasma}}}}{{{\mathrm{d}}t}} = A_{\mathrm{GItract}}F_{\mathrm{a}}k_{\mathrm{p}} - A_{\mathrm{plasma}} \times k_{\mathrm{m}} - A_{\mathrm{plasma}} \times k_{\mathrm{e}}$$
(3)
$$\frac{{{\mathrm{d}}A_{\mathrm{faeces}}}}{{{\mathrm{d}}t}} = A_{\mathrm{GItract}} \times (1 - F_{\mathrm{a}}) \times k_{\mathrm{p}}$$
(4)
$$\frac{{{\mathrm{d}}A_{\mathrm{elim}}}}{{{\mathrm{d}}t}} = A_{\mathrm{plasma}} \times k_{\mathrm{e}}$$
(5)
$$\begin{array}{l}\frac{{{\mathrm{d}}A_{\mathrm{GItract}}}}{{{\mathrm{d}}t}} = I\left( t \right) - k_{\mathrm{p}} \times F_{\mathrm{a}} \times A_{\mathrm{GItract}} - \left( {1 - F_{\mathrm{a}}} \right) \times k_{\mathrm{p}} \times A_{\mathrm{GItract}} \\ \kern3pt \quad\qquad = I\left( t \right) - k_{\mathrm{p}} \times A_{\mathrm{GItract}}\end{array}$$
(6)

The amount of cannabinoid transferred into milk is then given by

$$A_{\mathrm{milk}} = \int _0^t A_{\mathrm{plasma}}\left( t \right) \times k_{\mathrm{m}}\,{\mathrm{d}}t$$
(7)
$$A_{\mathrm{milk}{\_}\mathrm{daily}} = A_{\mathrm{milk}{\_}\mathrm{day}_i} - A_{\mathrm{milk}{\_}\mathrm{day}_{i - 1}}$$
(8)

For the complete equation set, see the description of model A in the Supplementary Information.

Model B is a simplification of model A, where the plasma becomes a central compartment with only one output to milk. In this case, the absorption parameter becomes Fb, that is, the absorbed fraction bioavailable for milk excretion. The GI tract output becomes the total elimination, accounting not only for faecal excretion but also for putative physiological elimination processes such as urinary excretion, metabolization and other possible processes. The ordinary differential equation system in this case becomes

$$\frac{{\mathrm{d}}A_{\mathrm{milk}}}{{\mathrm{d}t}} = A_{\mathrm{central}} \times k_{\mathrm{m}}$$
(9)
$$\frac{{{\mathrm{d}}A_{\mathrm{central}}}}{{{\mathrm{d}}t}} = A_{\mathrm{GItract}} \times F_{\mathrm{b}} \times k_{\mathrm{p}} - A_{\mathrm{central}} \times k_{\mathrm{m}}$$
(10)
$$\frac{{{\mathrm{d}}A_{\mathrm{elim,tot}}}}{{{\mathrm{d}}t}} = A_{\mathrm{GItract}} \times (1 - F_{\mathrm{b}}) \times k_{\mathrm{p}}$$
(11)
$$\frac{{{\mathrm{d}}A_{\mathrm{GItract}}}}{{{\mathrm{d}}t}} = I\left( t \right) - k_{\mathrm{p}} \times A_{\mathrm{GItract}}$$
(12)

The amount of cannabinoid transferred into milk is

$$A_{\mathrm{milk}} = \int_0^t A_{\mathrm{central}}\left( t \right) \times k_{\mathrm{m}}\,{\mathrm{d}}t$$
(13)

The complete equation set can be found in the description of model B in the Supplementary Information.

Obtaining optimized model parameters

The models were implemented using Python3 and then fitted with the lmfit module to the data to obtain parameters for the predictive model. The lmfit module (https://lmfit.github.io/lmfit-py/) provides a solid interface to perform nonlinear least-squares minimization. The code is available in the Supplementary Information, as is an executable Python file with the predictive model (PredictiveModel.zip). We used the Levenberg–Marquardt algorithm to solve the minimization problem, with the distance function being the squared difference between the dataset points (amount of cannabinoid excreted with milk/faeces or present in milk on a given day) and the corresponding modelled amount. For each cannabinoid, we fitted one model to each cow dataset, obtaining a range of parameter estimates (Supplementary Fig. 2). Cow H5 model was not included in the final average parameters estimation because of the subclinical ketosis it developed during the experiment. As the available data (non-milk data were available only for two cows in each group) were not sufficient for the solver to find a minimum for the objective function (squared sum of residuals), we placed a boundary on the parameter space. To obtain the final parameter estimation, we used the following procedure:

Model A:

  1. 1.

    Fix kp = 1.3 d−1.

  2. 2.

    Fit model to each single cow dataset (initial estimate: Fa = 0.8, km = 10, ke = 3,000).

  3. 3.

    Calculate \(\overline {F_{\mathrm{a}}}\) using only Fa from models fitted to dataset with non-milk data.

  4. 4.

    Set \(F_{\mathrm{a}} = \overline {F_{\mathrm{a}}}\) for next fit initial estimate.

  5. 5.

    Set Fa not to vary in fit for cows without non-milk data.

  6. 6.

    Fit each single cow dataset to model A with updated initial estimate.

  7. 7.

    Calculate \(\overline {k_{\mathrm{e}}}\) (geometric mean).

  8. 8.

    Set upper bound for \(k_{\mathrm{e}} = \overline {k_{\mathrm{e}}}\).

  9. 9.

    Fit model to each single cow dataset with updated constraints

  10. 10.

    Take the average of each parameter (using the geometric mean where appropriate), excluding cow H5 (subclinical ketosis), as the final model parameters.

Model B:

  1. 1.

    Fix kp = 1.3 d−1.

  2. 2.

    Fix km = 1 d−1.

  3. 3.

    Fit each single cow dataset milk data to model B to obtain Fb (initial estimate: Fb = 0.01).

  4. 4.

    Final model parameter is \(F_{\mathrm{b}} = \overline {F_{\mathrm{b}}}\) (arithmetic average excluding cow H5, subclinical ketosis)

Predictive model

The optimized model parameters were used to create a predictive cannabinoid milk transfer model for each of the selected cannabinoids. Two Python files, containing an implementation of the model, are included in the Supplementary Software (PredictiveModel.zip). To run the model, open the exec.py file in a Python3 console, select a cannabinoid by assigning to the CANNABINOID variable the correct string, set the DOSE, EXPOSURE_DURATION and DEPURATION values and run the code. It will return a plot with the selected cannabinoid amount in milk over the selected exposure scenario.

Assessment of consumer exposure to ∆9-THC using the RACE tool

To evaluate the concentration of ∆9-THC we measured in cow milk, the EFSA RACE software tool was used. The RACE tool estimates acute and chronic exposure from single foods using food consumption information from the EFSA Comprehensive European Food Consumption Database (Comprehensive Database), comparing it to relevant toxicological reference points. The Comprehensive Database summarizes existing national information on food consumption at the individual level and is the most complete and detailed collection of food consumption data currently available for the European Union. Among other things, the output of the tool is the percentage exhaustion of the toxicological reference value for various population groups35. For the assessment, the mean and maximum ∆9-THC concentrations in milk during the adaptation and exposure period were used. Analogous to previous exposure estimates6, it was assumed that this concentration would be found in milk and milk products made from it. Due to the delayed build-up of ∆9-THC in milk, the values of the last day of the adaptation period and the second to sixth day of the exposure period were used to calculate the means. The ARfD of 1 µg per kg body weight derived by the EFSA6 was applied as the toxicological reference value.

Reporting summary

Further information on research design is available in the Nature Portfolio Reporting Summary linked to this article.