Quantitative modeling of radioactive cesium concentrations in large omnivorous mammals after the Fukushima nuclear power plant accident

Large quantities of radionuclides released by the Fukushima nuclear power plant accident entered terrestrial and marine ecosystems. The resulting radioactive contamination of large omnivorous wild mammals such as wild boar (Sus scrofa) and Asian black bear (Ursus thibetanus) varied greatly depending on location, season, and time after the accident. Quantitative modeling of how such factors influence radionuclide burdens in these species is important for enhancing current knowledge of chronic radionuclide exposure consequences in mammalian populations, and for assessing potential human risks from consumption of contaminated animal meat. Here we modeled the time course of radioactive cesium (134Cs + 137Cs) concentrations in boar and black bears from Fukushima Prefecture over ~ 7 years after the accident, using nonlinear robust and quantile regressions and mixed-effects modeling. To estimate predictive performance, models fitted to the full data set were compared with those fitted only to the first 3.5 years of data, and tested on the last 3.5 years of data. Ecological half-lives for radioactive cesium, and magnitudes and phase shifts for sinusoidal seasonal oscillations in cesium burdens, were estimated by each analysis method for each species. These results can improve the understanding and prediction of radionuclide concentrations in large mammals that inhabit radioactively contaminated areas.

We ln-transformed the cesium activity values to bring their distribution closer to normal, creating the variable LnCsTot. To facilitate regression analyses (described below), we removed instances of missing data and cesium levels below detection: 20 samples (1.4%) for boar and 15 samples (3.3%) for bears. The time when each sample was taken (labeled "Day of collection" in the Fukushima Prefecture Government data set) was converted to years since the Fukushima accident (since March 11, 2011), assuming that 1 year = 365.25 days. This time of sample collection in years was called variable T.
Since for each sample some time passed between sample collection and radioactivity measurement (labeled "Result found Date", called T r in our notation), we needed to correct the reported LnCsTot values for physical decay over this time, which was different for different samples. The procedure used to perform this correction is described in Supplementary methods. The data with corrected total cesium values (LnCs c ) are provided in Supplementary data (Supplementary_Dataset_File_Full).
Mathematical model. To describe the data on ln-transformed total radioactive cesium levels (LnSc c ) in each species as function of time after the accident (T), we developed the following simple mathematical model (Eqs. 1A, 1B): Here the term X represents the estimated average radioactive cesium level in the studied area, based on the intercepts (LnCs134 t0r for 134 Cs and LnCs137 t0r for 137 Cs, respectively) from robust regression discussed in Supplementary methods, and taking into account only physical decay for each isotope (with half-lives of Th Cs134 for 134 Cs and Th Cs134 for 137 Cs, respectively). The terms Q, µ, ν, A and P represent adjustable model parameters. Parameter Q represents the fitted relationship between radioactive cesium levels in the animal (Bq/kg), relative to the external environment (Bq/m 2 ). Parameter µ represents the net rate of radioactive cesium reduction in animal tissues over time due to all processes except physical decay (e.g. decrease in bioavailability due to migration of cesium into deeper soil layers, human-mediated cleanup efforts, etc.). Parameter ν is a potential power dependence for these processes. By default, ν was set to ν = 1, but exploratory calculations using ν = 2 or treating ν as a freely adjustable parameter (≥ 0.1) were performed as well. Parameters A and P in the sine function represent a sinusoidal approximation for seasonal changes in radioactive cesium levels in animal tissues (e.g. due to seasonal variations in diet and life style), where A is the amplitude of the oscillations, P is the phase shift, and the period is set to 1 year. For simplicity, these parameters were assumed to be the same for both studied cesium isotopes. The descriptions of each parameter are also presented in Table 1.
Model fitting approaches. Initially, we used nonlinear ordinary least squares (OLS) regression (nls R function) to fit the model (Eq. 1A, 1B) to the data. To find the global optimum fit, we repeated the fitting procedure 2000 times with slightly different random initial parameter values and recorded the solution with the smallest root mean squared error (RMSE). Diagnostics on this regression included checking of convergence criteria and analyses of residuals (by scatter plot and histogram, regressing residuals as function of T, visualizing the QQ plot, autocorrelation and partial autocorrelation functions with 95% confidence intervals, performing the Shapiro-Wilk normality test, and calculating skewness and kurtosis). For boar data, diagnostics revealed problems with convergence (both X-convergence and relative convergence) and non-normality of residuals: e.g. Shapiro-Wilk p-value = 1.476 × 10 -7 , skewness = − 0.37, kurtosis = 3.50. For black bear data similar problems (1A) www.nature.com/scientificreports/ occurred with convergence, but residuals were closer to the normal distribution (perhaps due to smaller sample size): e.g. Shapiro-Wilk p-value = 0.0526, skewness = − 0.058, kurtosis = 2.45. Due to these issues, we used robust nonlinear regression (nlrob R package) to reduce the effects of "outlier" data points. To find the global optimum, we repeated the fitting procedure 2000 times with slightly different random initial parameter values and selected the solution with the smallest absolute value of median residuals. The best-fit parameters for OLS and robust regressions were somewhat different for both boar and bear data. For boar data, the minimum robustness weight was 0.339 and the median was 0.762, and the corresponding values for black bear data were 0.557 and 0.821, respectively.
For each species, we compared the performances of model variants with different assumptions about parameter ν: (1) The default case with ν = 1, which represents an exponential rate of radioactive cesium decrease due to processes other than physical decay. (2) The case with ν = 2, which represents quadratic decay. (3) The case with ν being freely adjustable (≥ 0.1). The comparisons were based on Akaike information criterion (AIC) 22,23 . The purpose of these calculations was to better assess the shape of the time course for non-physical factors involved in radioactive cesium level decline in animal tissues over time after the accident.
In addition to analyzing the full data set for each species, we also performed separate analyses on subsets of data from specific locations: from those districts of Fukushima Prefecture where the mean radioactive cesium levels in animal samples were the highest, and where a sufficiently large number of samples was present. For wild boar, the two selected districts for this subset analysis were Soso and Kenpoku (819 samples), and for black bear they were Kenpoku and Kenchu (163 samples).
To further assess the sensitivity of model results to geographical and temporal factors, we also constructed a separate subset of data for each species. This subset excluded data from the Aizu and Minamiaizu districts, which are separated by mountains from the Fukushima Daiichi Nuclear Power Plant, and excluded data collected ≤ 6 months after the accident. These restrictions were intended to determine model performance on data collected in a more geographically contiguous area after the initial abrupt changes in contamination levels were completed and the system entered the phase of more stable kinetics. The purpose of all these analyses was to assess whether the time course of radioactive cesium levels in the bodies of each species differed between locations with high contamination vs. those with lower contamination, and as function of time after the accident.
We were interested in quantifying not only the center of the distribution of radioactive cesium values in each species over time, but also in assessing the lower and upper tails of this distribution. For this purpose, we fitted the model (Eq. 1A, 1B) for each species using quantile regression (nlrq function in quantreg R package) for the median (50th percentile), and also for the 25th and 75th percentiles. Initial parameter estimates for the quantile regressions were taken from best-fit parameters from robust regression described above. The 25th and 75th percentiles were selected instead of more extreme values (e.g. 5th and 95th) because the latter resulted in poor quality fits due to limited amounts of data at the fringes of the distribution.
To assess the variability of model parameters by location in more detail, we used mixed effects modeling (nlme R package) on the data from each species. Since original OLS fits suggested substantial deviations of residuals from the normality assumption, we performed mixed effects modeling on data with some outlier data points removed. The OutlierDetection package in R removed 43 boar samples and 5 bear samples. These outliers are listed in the Supplementary_outlier_data_points file. The remaining samples were used for mixed effects model fitting, but model performance metrics like coefficient of determination (R 2 ) and RMSE were assessed on the full data set (with outliers included) for each species.
Since the Fligner-Killeen test of homogeneity of variances by district generated low p-values for both species (4.6 × 10 -14 for boar and 0.018 for black bear), we allowed modelled variances to differ by district (using the weights option in nlme). We investigated several random effects structures for some or all model parameters, with randomness by district only, or by district and municipality within district. Model diagnostics were the same as for fixed effects OLS modeling described above, and also included boxplots of model residuals by district. The mixed effects model variants with different random effects structures were compared using the anova function in R, and also by assessing convergence criteria, normality of residuals, skewness, and kurtosis. Consequently, preferred mixed effects model variants were selected for the full data as well as for the subset of two districts with high radioactive cesium levels, separately for each species.
Model extrapolation from training to testing data. To investigate how the robust and quantile regression fits of our model could extrapolate beyond the time range that was used for model fitting, we split the data for each species into "training" (early times) and "testing" (later times) parts. The split was done based on time since the accident (T variable), so that approximately ½ of the samples were assigned to the training and testing sets, respectively. For wild boar data, the training set included times between 0.20 and 3.45 years after the accident, and the testing set included times between 3.45 and 7.03 years. For black bear data, the training set included times between 0.42 and 3.46 years after the accident, and the testing set included times between 3.46 and 6.87 years.
We also evaluated an alternative approach to splitting the data, where the split was done randomly instead of by time. In other words, any data point regardless of time had an equal probability of being assigned to either the training or the testing data set. Both the training and testing data subsets generated by this random split included the complete time range. This approach was implemented in context of the sensitivity analysis described above.
For each species, robust and quantile regressions were fitted to training data, and their predictions were calculated for testing data. For robust regression, RMSE was calculated on testing data for two scenarios: (1) for the model fitted to training data only, and (2) for the model fitted over the entire data range (training + testing). These RMSE values for conditions 1 and 2 were compared to assess the quality of model extrapolation.

Application of the model to wild boar data from the Chernobyl accident area.
To compare the results of our analysis of wild boar contamination with radioactive cesium in the area affected by the Fukushima accident with data from another location, we also analyzed wild boar data from the Chernobyl accident area. These data were published by Gulakov 14 and contain summaries of 137 Cs contamination levels in the muscles of 188 boar collected between 1991 and 2008 (i.e. from 5 to 22 years after the 1986 accident). Sampling was carried out in three zones with different land contamination levels with 137 Cs. This data set provides important information on radioactive cesium contamination in wild boar in the Chernobyl area. Unfortunately, 137 Cs measurements in each sampled boar were not provided by Gulakov 14 , and only summary statistics are available for each zone and year after the accident (Tables 1-3 in reference 14 ): number of animals, mean, minimum and maximum 137 Cs levels.
We could not apply the full model (Eq. 1A, 1B) to these summary data which lacked seasonality information and 134 Cs data. However, we were able to perform a weighted linear regression to quantify the ecological half-life of 137 Cs in Chernobyl boar and the relationship between 137 Cs levels in the animals (Bq/kg), relative to the external environment (Bq/m 2 ). The data used for this analysis, derived from Gulakov 14 , are provided in Supplementary data (Supplementary_Dataset_File_Full). They contain the following variables. Zone = location of sample collection (Alienation, Permanent control or Periodic control). Time = time in years after the Chernobyl accident. LnMeanCs = ln-transformed mean 137 Cs level in boar muscle (Bq/kg). LnMeanCs_c = LnMeanCs − X, where X is ln-transformed 137 Cs land contamination (Bq/m 2 ) in the given zone, corrected for physical decay of 137 Cs. Weight = weighting of each data point used for regression. Weight = number of animals/(ln[maximum 137 Cs level] − ln[minimum 137 Cs level]) 2 . These approximately inverse-variance weights were normalized by the overall mean across all data points, so that the mean weight across all data points was set to 1.
These data were analyzed by weighted linear regression in R, where LnMeanCs_c was allowed to depend on Time and Zone variables. Model variants containing all possible combinations and pairwise interactions between these predictor variables were fitted and their performances were compared using the Akaike information criterion with correction for small sample size (AICc). These calculations were performed using the glmulti R package. Multimodel inference (MMI) was performed on this collection of fitted model variants. It resulted in the calculation of model-averaged parameter estimates, 95% CIs and importance scores, corrected for model selection uncertainty. Of main interest here were the intercept parameter, which is analogous to parameter Q in the full model (Eq. 1A, 1B), and the Time parameter, which is analogous to parameter µ in the full model. The ecological half-life for 137 Cs was calculated based on the Time parameter.

Results
We analyzed the wild boar data with robust regression using variants of our model formalism (Eq. 1A, 1B) with different power dependences on time: with parameter ν = 1, ν = 2, or freely adjustable ν. These variants represent different assumptions about the shape of the time dependence for the reduction in radioactive cesium levels in boar tissues over time after the accident due to all processes other than physical decay. The variants with ν > 1 had better support by AIC on the full data set than the ν = 1 variant, by up to 41.6 AIC units for the freely adjustable ν variant (where best-fit ν was 2.750 with standard error (SE) of 0.796). However, the simplest variant with parameter ν = 1, which assumes that reduction in radioactive cesium levels in boar tissues over time after the accident due to all processes other than physical decay follows simple exponential kinetics, had the best performance in terms of extrapolating from training to testing portions of the data set (i.e. from early times after the accident to later times). The RMSE value on testing data for this model version fitted to the training data only, divided by the RMSE value on testing data for the model version fitted to all the data (training plus testing) was only 1.062. In comparison, much higher RMSE ratios (> 1.4), which indicate poor extrapolation performance, were found for model variants with ν = 2, or with freely adjustable ν.
These results suggest that despite its worse AIC scores on the full data set, the simplest model version with ν = 1 had the best predictive power for extrapolating from short times to longer times. Due to this important property, which is visualized in the plots showing extrapolation of model predictions (Figs. 1, 2) and in the similarity of best-fit parameters for fits to the entire data set vs. fits to training data only (Tables 2, 3), and also to the model's simplicity, we used this ν = 1 assumption for the quantile and mixed effects regression analyses of the boar data. The results of these analyses are shown in Figs. 1, 2 and in Tables 2, 3, 4.
Model extrapolations from training to testing data of course resulted in some variations in estimated model parameters (Tables 2, 3 and Supplementary Tables 1-2). The variations were larger when the training/testing split was done based on time, compared with the sensitivity analysis, where the split was done randomly, as described in Materials and Methods. Although these parameter variations caused clearly visible effects when best-fit predictions were plotted graphically (Figs. 1, 2, 3, 4, Supplementary Figs. 1-2), they were generally not large, compared with corresponding parameter uncertainties (standard errors). This conclusion is supported by calculations of the ratio of parameters for fits to training data only / full data, which are shown in Tables 2, 3  and Supplementary Tables 1-2. The best-fit parameters and visualization of model behaviors show that radioactive cesium levels in boar tissues continuously decreased over time after the accident. This decreasing trend could be reasonably extrapolated using our model from early to late times after the accident. Understandably, extrapolation quality was lower on subsets of the data for specific districts, especially for the tails of the distribution (25th and 75th percentiles). Similar results were obtained by our sensitivity analysis (described in Materials and Methods), where only a subset of data was used for each species (Supplementary Tables 1-2 and Supplementary Figs. 1-2). These findings www.nature.com/scientificreports/ suggest that the radioactive cesium burden in the boar population of Fukushima Prefecture is generally decreasing, but stable or increasing trends are potentially possible in some areas and/or for a fraction of individuals. Perhaps these patterns are caused by variations in bioavailability of cesium by location and/or different food source preferences between boar individuals. Substantial seasonal oscillations in radioactive cesium levels in boar were evident, and they were reasonably modeled by the simple sinusoidal function (including this function in the model removed any significant autocorrelation of residuals). Peak radioactive cesium levels in boar tissues are predicted to occur in winter and early spring, probably due to increased reliance on food sources prone to cesium accumulation during these seasons. These findings are consistent with previous analyses 16 , but extend them to longer times after the accident.
Best-fit parameters for the robust, quantile and mixed effects analyses of boar data are compared in Table 5. For each method and subset of data, we calculated the ecological half-life of radioactive cesium in boar, defined as ln [2]/µ. As shown in Table 5, all regression methods produced similar ecological half-life values ranging from 3.3 to 5.7 years, but with considerable uncertainties. Analysis of the subset of data from specific districts   www.nature.com/scientificreports/ (Soso and Kenpoku) with high radioactive cesium contamination produced somewhat higher estimates for the ecological half-life, from 4.9 to 13 years, but the uncertainties were expectedly larger due to reduced sample size, and included infinity in some cases ( Table 5). The sensitivity analysis also produced similar ecological halflife estimates ( Table 5). The seasonality effects (amplitude and phase shift) were very similar for all regression methods, for the full data set and for the subsets of locations ( Table 5). The parameter estimates shown in Table 5 suggest that radioactive cesium levels in wild boar in the Fukushima accident area vary by several orders of magnitude between boar individuals and are strongly influenced by seasonal effects. Importantly, the overall trends of decreasing radioactive cesium levels over time were detected by the implemented analyses and data selection approaches. The ecological half-life estimates and their uncertainties quantify these trends numerically. Of course, the inherent differences in assumptions and model fitting methodology between different regression methods lead to some variability in ecological half-life estimates. For example, robust regression estimated only the central tendency of the data, whereas quantile regression estimated three percentiles: 25th, 50th and 75th. In both robust and quantile regressions, the effects of geographical location were not explicitly modeled, and the entire area where data were collected was in essence treated as a single homogeneous unit. In contrast, mixed effects modeling accounted for the effects of location in more detail by allowing model parameter estimates to vary between different districts and municipalities within each district. However, compared with mixed effects modeling, robust and quantile regressions were less sensitive to "outlier" data points with extremely low or extremely high radioactive cesium levels. Each method therefore has its strengths and weaknesses, and model parameter estimates generated by each method should be compared and interpreted together to identify common tendencies. Table 4. Parameter values for fitting our model (Eq. 1A, 1B) by mixed-effects regression to wild boar and Asian black bear data. The preferred mixed-effects model structure for each data set was selected based on several criteria described in the main text. Random effects for selected parameters were allowed to vary by district and municipality within district, and the variance was allowed to be different for each district. SE represents standard errors and SD represents standard deviations.

Species
Para-meter  www.nature.com/scientificreports/ The analogous analyses of black bear data also showed that the simplest model variant with parameter ν = 1 had the best extrapolation performance from training to testing data, with an RMSE ratio of only 1.04. The other variants had worse predictive power and worse AIC scores (by up to 2.7 units). The results of analyzing black bear data by robust, quantile and mixed effects regressions with the ν = 1 assumption are shown in Tables 2, 3, 4 and Figs. 3, 4. Best-fit parameters for all these models are compared in Table 6. These findings suggest that the ecological half-life values for radioactive cesium in black bears were not too different from those in wild boar, considering the large uncertainties (which were larger for the bear data due to lower sample size). The amplitudes of seasonal oscillations in cesium levels were also relatively similar for the two species, but the phase shifts tended to be different (Tables 5, 6). Details of mixed-effects model fits to wild boar and black bear data are provided in Supplementary material online (Supplementary_mixed_effects_model_output).
Parameter Q, which is related to accumulation of radioactive cesium from the environment by each species, tended to be lower for black bears than for wild boar (Tables 2, 3). This finding is consistent with previous  Table 5. Summary metrics for robust, quantile and mixed-effects fits of our model (Eq. 1A, 1B) to wild boar data. In this and the following tables, infinite values for the ecological half-life of radioactive cesium indicate instances where the uncertainty on the corresponding model parameter µ extended to negative values. Sensitivity analysis represents analysis of a subset of data collected ≥ 6 months after the Fukushima accident, and where Aizu and Minamiaizu districts were excluded. www.nature.com/scientificreports/ analyses 16 . The differences in best-fit Q values between the 75th and 25th quantiles were similar for both species: 1.565 for boar and 1.573 for black bears. These results suggest that the magnitude of variability of radioactive cesium burdens between individuals was similar in Fukushima boar and black bear populations. We also applied our modeling concepts to 137 Cs contamination data in wild boar from a completely different location: the Chernobyl accident area. As described in "Materials and methods", these data were published by Gulakov 14 , and we analyzed them by weighted linear regression and multimodel inference (MMI). The bestsupported weighted linear regression model variant, based on sample size corrected Akaike information criterion score, contained Time and Zone variables as predictors, without an interaction between them. The fitted coefficient for the Zone variable represents the relationship between mean ln-transformed radioactive cesium levels in boar tissues and in the external environment, which is analogous to parameter Q in our full model (Eq. 1A, 1B). Its best-fit values were: Alienation zone = − 5.552 ± 0.234 (standard error); Permanent control zone = 0.955 ± 0.186 units higher; and Periodic control zone = 1.438 ± 0.201 units higher. The fitted coefficient for the Time variable represents the time-dependent decrease in radioactive cesium level in boar tissues due to all processes except physical decay, which is analogous to parameter µ in our full model (Eq. 1A, 1B). Its best-fit value was 0.033 ± 0.020 years −1 , which was not significantly different from zero (p-value = 0.115). The results of the MMI analysis, which included not only a single best-supported model, but averaged across all tested model variants, are shown in Supplementary Table 3 and in Supplementary Fig. 3.
These findings suggest that the relationship between mean ln-transformed radioactive cesium levels in boar tissues and in the external environment (Q) is similar in order of magnitude between Fukushima and Chernobyl data. However, the time-dependent decrease in radioactive cesium level in boar tissues (µ) in Chernobyl had 95% CIs that overlapped zero (Supplementary Table 3). The resulting ecological half-life estimate for 137 Cs based on MMI analysis of Chernobyl data was 113.2 (95% CI: 18.5, ∞) years, and it was ∞ based on the single best-supported model variant mentioned above. This result suggests that ecological processes of radioactive cesium reduction in wild boar appear to be more rapid in Fukushima than in Chernobyl, although there is some overlap in the uncertainties of the estimates from these locations (e.g. the uncertainties of some ecological halflife estimates from Fukushima also included infinity, Table 5).

Discussion
Our modeling analysis of wild boar and Asian black bear data from Fukushima Prefecture over approximately 7 years after the Fukushima Daiichi Nuclear Power Plant accident, using robust and quantile regressions and mixed-effects modeling techniques, generated detailed quantitative estimates for important parameters such as the ecological half-lives and seasonal oscillation magnitudes of radioactive cesium levels in each species (Tables 2,  3 , 4, 5, 6). The mixed-effects models, which included the effects of geographical location, achieved relatively strong performance (R 2 > 0.6 and RMSE < 1 Bq/kg on a natural log scale, Table 4), despite the simplicity of the underlying model structure (Eq. 1A, 1B) and only time and location as predictors.
Wild boars were quite abundant in the contaminated area, and > 1000 individuals killed by hunters provided detailed data on radioactive cesium burdens in this species at different times and locations. Our analysis estimated an ecological half-life for radioactive cesium in this species to be < 10 years, but with considerable uncertainties (Table 5). In contrast, in Europe some boar continued to have very high radioactive cesium concentrations in their tissues even several decades after the Chernobyl disaster, although the average levels of cesium isotopes in www.nature.com/scientificreports/ soil decreased dramatically over this time due to ecological processes and physical decay 14,15,[24][25][26] . The finding of slow decline or even no decline in radioactive cesium levels due to ecological factors in European boar was also supported by our analysis of Chernobyl boar data (Supplementary Table 3 and Supplementary Fig. 3). The persistence of high radioactive cesium levels in European boar populations despite their decline in soil was probably caused by increased availability of these radionuclides in boar diet over time. For example, some fungi (e.g. genus Elaphomyces) strongly accumulate radioactive cesium from soil and can serve as sources for boar contamination because these fungi are eaten by boar, particularly during colder months when green vegetation is not available. This process of boar contamination through consumption of fungi appears to be present, but considerably less pronounced in Japan, compared with Europe 16,19,24 . In addition, the level of radioactivity following the Fukushima accident decreased faster than in Chernobyl due to the effects of runoff and decontamination 27 .
Ecological half-lives of radioactive cesium in Fukushima Asian black bears were also estimated to be < 10 years, although uncertainties were larger than in boar because of smaller sample sizes (Table 6). A direct comparison with Europe is not possible for this species because it is found only in Asia, and has different ecology and diet than the European brown bear (Ursus arctos). Unlike wild boar, which feed actively throughout the year, Asian black bears in the Fukushima area hibernate in dens from November to April. Such periodic dormancy affects the seasonal oscillations in radioactive cesium levels in bear tissues, but our results (e.g. Figs. 1, 2, 3, 4) suggest that the sinusoidal pattern provides a decent description of these patterns. Qualitatively similar conclusions were reached in previous analyses of black bear data from Fukushima 16 .
The magnitudes and phases of seasonal oscillations in radioactive cesium levels in boar and black bears found in this study (Tables 5, 6) were similar to those reported by previous investigations 16 . Our estimates suggest that the seasonal oscillation amplitudes in these species are approximately 1.3 to 2.0 on the natural log scale (Tables 5,  6), which correspond to ~ 3.7 to 7.4-fold changes on the linear scale.
In conclusion, we developed a simple mathematical model of radioactive cesium ( 137 Cs + 134 Cs) levels in animal tissues (Eq. 1A, 1B) and fitted it to the most recent publicly available version of wild boar and Asian black bear data from Fukushima Prefecture, which covers approximately 7 years after the accident. We used several techniques to fit the model (robust regression, quantile regression and mixed effects modeling), which have different assumptions and provide complementary information. Specifically, model parameters estimated by robust regression are not very sensitive to "outlier" data points (i.e. extremely high or low radioactive cesium concentrations in some wild boar or black bear individuals). Parameters estimated by quantile regression are also relatively insensitive to outliers, and not only the central tendency, but several percentiles of the data distribution can be modeled. The strength of mixed effects modeling, as it was used here, is in its detailed treatment of parameter variations by geographical location, but its drawback is greater sensitivity to outliers than the previous methods. Importantly, ecological half-lives for radioactive cesium, and magnitudes and phase shifts for sinusoidal seasonal oscillations in cesium burdens, were estimated by each regression method for each species.
The study has limitations because many potentially important variables (e.g. age and sex of the animals, their home ranges and migration patterns) were not measured and, therefore, could not be included into the model. Aging and life span in wild boar were studied in European populations 28,29 , as well as in Japan 11,30 . In the Fukushima accident area, boars of all ages, from < 3 months to > 4 years, were found, and some studies reported a marginally statistically significant tendency for higher radioactive cesium levels in older individuals 11 . Unfortunately, information on age of the boar or black bears was not available in the data sets analyzed here, so age could not be incorporated as a predictor in the model.
Our approach was intended to describe the general tendencies in radioactive cesium contamination levels (e.g. median, 25th and 75th percentiles) over time in the studied species, but it cannot accurately predict the level in any specific individual. Despite these limitations, we believe that the study results can improve current understanding and prediction of radionuclide concentrations in large omnivorous mammals from radioactively contaminated areas.