Perinatal mortality after the Fukushima nuclear accident: An ecological study

Objective This study continues former studies on perinatal mortality in Japan after the Fukushima Daiichi nuclear power plant (FDNPP) accident in March 2011. An increased study region is chosen, and the study period is extended to 2019. Methods Japanese monthly perinatal mortality data are provided on a prefecture level by the Japanese government. The study region consists of 12 prefectures around the FDNPP; the rest of Japan is used as the control region. A combined non-linear regression of perinatal mortality rates in the study- and control regions is conducted. The regression model allows for a common asymptotic lower limit of perinatal mortality, seasonal variations, and periodic peaks in 2012–2019 in the study region. To determine the dependency of the effect on distance from the FDNPP, the study region is divided into four core prefectures and eight prefectures surrounding the core prefectures. Results Perinatal mortality rates in the study region show a significant 6.4% (95% CI: 1.8%, 13.4%) overall increase in 2012–2019 relative to the trend in preceding years with no attenuation during 2012–19. The increase translates to 590 (165, 1226) excess perinatal deaths (p = 0.016). It is characterized by annual peaks with maxima in April. A 13.6% increase is determined in the four core prefectures and a 4.3% increase in eight prefectures surrounding the core prefectures. Before 2012, there is a peak around April 2011 and a decline in October 2011; another significant peak is detected in November 2012. In the 4 core prefectures, large increases are found in the first quarter of 2018 (+70%) and in May 2019 (+130%). Conclusion This study finds periodic peaks in perinatal mortality in spring 2012–2019 in 12 prefectures of Japan surrounding the FDNPP. In light of massive increases in 2018 and 2019 in the four core prefectures, continued investigation of perinatal mortality in contaminated regions of Japan is recommended.


Introduction
After the accident at the Fukushima Daiichi nuclear power plant (FDNPP) in March 2011, little attention was paid to possible adverse effects on pregnancy outcomes in Fukushima and neighboring prefectures. The 2013 UNSCEAR report on Fukushima stated that prenatal exposures from the accident at FDNPP "were not expected to increase the incidence of spontaneous abortions, miscarriages, perinatal mortality, congenital effects or cognitive impairment" [1], although there were several reports on teratogenic radiation effects after the Chernobyl accident in 1986. An official German study [2], conducted by the Federal Office of Radiation Protection (BfS), investigated perinatal mortality rates in Bavaria which was the German state with the highest mean cesium soil contamination. In the southern part of Bavaria (south of River Danube), the cesium soil contamination was three times higher than north of River Danube. Any increase in perinatal mortality after Chernobyl relative to the expected trend should have been greater in southern Bavaria than in northern Bavaria. In both parts of Bavaria, there was no significant excess in perinatal mortality in any of eight 3-month periods from June 1986 through May 1988. However, a trend analysis of German annual perinatal mortality rates found a statistically significant 5% increase in 1987, one year after the Chernobyl disaster [3]. The result was later confirmed by Scherb et al. [4]. The analysis of the monthly data detected peaks in perinatal mortality rates at the beginning and end of 1987 which were associated with peaks of cesium concentration in pregnant women, lagged by seven months [3]. The effects were interpreted as teratogenic radiation effects.
After the Fukushima nuclear accident in March 2011, a survey of stillbirth rates, pre-term births, low birth weight, and congenital anomalies by the Fukushima Radiation Medical Science Center did not find significant regional differences in stillbirth rates within Fukushima Prefecture [5]. Scherb et al. examined the trend of perinatal mortality before and after the Fukushima nuclear accident in six prefectures near the Fukushima Daiichi nuclear power plant (FDNPP) [6]. Using linear logistic regression, they determined a highly significant upward shift in perinatal mortality rates in January 2012 through December 2014 relative to the expected trend before 2012. Körblein and Küchenhoff [7] conducted a combined regression of annual perinatal mortality rates in a study region (Fukushima Prefecture plus 6 nearby prefectures) and the rest of Japan (the control region). To model the long-term trend of perinatal mortality, they used individual exponential trends with a common constant that represented a natural lower limit of perinatal mortality. They found a significant upward shift in perinatal mortality in 2012-2015, but the effect was substantially smaller than reported in [6]. A subsequent study by the same authors used monthly data of perinatal mortality through December 2017 from a study region consisting of Fukushima plus four adjacent prefectures; the rest of Japan served as the control region [8]. In 2012-2017, periodic peaks in perinatal mortality were detected with maxima in April. In Fukushima Prefecture, the effect was 3-times greater than in the four neighboring prefectures.
In the present study, the same statistical methods as in [8] are applied, but with data up to 2019 and with an expanded study region (12 prefectures). To determine the dependency of the effect size on distance from the FDNPP, the study region is subdivided into four prefectures near the FDNPP (Area A) and eight surrounding prefectures (area B). The enlarged study region allows a more accurate estimate of the overall impact of the Fukushima accident on perinatal mortality in Japan. These 12 prefectures were used in a recent study by the present author reporting a highly significant decrease in live births in December 2011, nine months after the Fukushima accident [9].

Material and methods
The numbers of live births, stillbirths, and early neonatal deaths (first 7 days), 2002 through 2019, are provided on a website of the Japanese government [10]. Stillbirths are defined in Japan as fetal deaths after 22 weeks of pregnancy. Perinatal mortality is defined as the number of stillbirths plus early neonatal deaths, divided by the number of live births plus stillbirths.
The prefectures of Fukushima plus three adjacent prefectures (Miyagi, Ibaraki, and Tochigi) are chosen as the core area (Area A). The surrounding eight prefectures (Iwate, Akita, Yamagata, Niigata, Gunma, Saitama, Tokyo, and Chiba) are used for comparison (area B). The remaining 35 prefectures of Japan constitute the control region (area C). A map of the study region is provided in Fig 1. The division of the study region is based on distance from the FDNPP rather than on radiation exposure. Table 1 contains the numbers of live births (LB), stillbirths (SB), and early neonatal deaths (NEO) in areas A, B, and C for two periods, 2002-2011 and 2012-2019.

Statistical analysis
First, a combined analysis of the data from the study region (12 prefectures) before March 2011 and the data from the rest of Japan (control region) in 2002-2019 is conducted. To model the time trends in the study-and control regions, exponentials with a constant of the form y = α+exp(β 1 +β 2 �t) are used, where t is time, β 1 and β 2 are individual trend parameters, and constant α is a natural lower limit of perinatal mortality which is assumed to be equal in the study-and control region. Under the Null-hypothesis of no effect of the Fukushima accident on perinatal mortality in the study region, Model (1) has the following form: Here, E(y(t)) is the expected value of perinatal mortality y(t); time t is defined as calendar year minus 2000, in fractions of a year (e.g. mid-January 2002 is t = 2+1/24), and dummy variable study denotes the data from the study region. The number of excess perinatal deaths in March 2011 through December 2019 is determined from the difference between observed and predicted deaths. This analysis is free of assumptions about trends in perinatal mortality after the Fukushima accident.
Next, a regression allowing for a level shift in January 2012 is conducted. The regression model is adopted from [8]. An indicator variable cp (change point) is defined with cp = 1 in January 2012 through December 2019 and cp = 0 before January 2012 and in the control region. Model (1.1) has the following form: Parameter β 5 estimates the size of the level shift in January 2012, and β 6 �(t−12) allows for a slope change in January 2012, i.e. an attenuation of the effect during 2012-2019.
Third, the regression model controls for seasonal variations and allows for periodic peaks of perinatal mortality in 2012-2019. Seasonality, as well as periodicity, is modeled by pairs of sine and cosine terms with periods of 6 and 12 months. Model (1.2) has the following form:

Exploratory analysis
The moving average of the residuals after the Fukushima accident exhibits noticeable peaks around April 2011 and November 2012, as well as a trough in October 2011. To avoid distortion of the long-term trend of the data, Model (1.2) is complemented by three bell-shaped excess terms (normal distributions) to fit the two peaks and the trough. Model (1.3) has the following form: To test the significance of the periodic terms in 2012-2019, the deviance obtained with Model (1.3) is compared with the deviance resulting from a regression without the periodic peaks. The p-value is determined by an F-test.

Combined regression of data from areas A, B, and C
To test whether there is a difference in effect size between areas A and B over the 2012-2019 period, data from Areas A, B, and C are analyzed in a combined regression. A change-point model analogous to Model (1.1) is used, with individual trend parameters (intercepts and slopes) in areas A, B, and C, and individual level shifts (cp) in areas A and B, but with a common slope change in January 2012 in areas A and B. Model (2.1) has the following form: Next, Model (2.2) is applied which allows for seasonal variations and periodic peaks in 2012-2019 in areas A and B, analogous to Model (1.2):

Exploratory analysis
The trend of perinatal mortality in Area A exhibits striking peaks in the first quarter of 2018 and in May 2019. To avoid distortion of the long-term trend, a dummy variable cp18 is added to Model (2.2). It is defined as cp18 = 1 in Jan 2018 through Dec 2019 in Area A and cp18 = 0 otherwise. In addition, bell-shaped excess terms are added to allow for the effects in April 2011, October 2011, and November 2012. A common parameter is used for the width of the bell-shaped excess terms since the regression with Model (1.3) had shown that they agreed within the limits of error. Model (2.3) has the following form: Iteratively reweighted non-linear regression with program R, function nls(), is used for data analysis and plotting, and a two-sided p-value <0.05 is considered statistically significant.

Study region
The results of the combined regression of the data from the study-and control regions with Model (1.0), i.e. parameter estimates, standard errors, t-values, and p-values, are listed in S1 Table in S1 File. S1 The regression with Model (1.1)-which includes a level shift (cp) and a slope change in January 2012 in the study region-leads to a significant improvement in the fit compared to a regression without level shift and slope change (p = 0.016, F-test). The estimated number of excess perinatal deaths from 2012-2019 is 590 (95% CI: 165, 1226), an overall increase of 6.4% (1.8%, 13.4%). The regression results are presented in Table 2  To test the significance of the periodic terms in 2012-2019, the deviance obtained with Model (1.3) is compared with the deviance resulting from a regression without the periodic peaks. The deviance is 431.8 (df = 412) without and 422.6 (df = 408) with the periodic terms (p = 0.065, F-test).
Eventually, the standardized residuals (see Fig 3 panel B) are checked for autocorrelation. There is no noticeable autocorrelation in the data from the study region nor from the control region for time lags of 1 to 11 months (see S3 Fig in S1 File).

Combined regression
The results of the combined regression of the data from areas A, B, and C with Model (2.1), which allows for individual level shifts in areas A and B in January 2012, are listed in Table 4

Exploratory analysis
The regression results with Model (2.3) are listed in S5 Table in S1 File. Table 5  An inspection of the list of parameters in S5 Table in S1 File shows that both the intercepts (parameters β 2 and β 3 ) and the slopes (parameters β 5 and β 6 ) in Areas A and B agree within the limits of error; the differences are negligible (p = 0.93 and p = 0.79, respectively). Moreover, there is no notable attenuation of the effect in 2012-1019; the estimate of the slope change is β 2 = 0.002 ± 0.013 (p = 0.89). Therefore, Model (2.3) is replaced by a more parsimonious and robust regression Model (2.4) with common parameters for the intercepts and slopes in areas A and B and without a slope change in January 2012. Regression with Model (2.4) yields slightly lower deviance (652.0) than model (2.3), although it requires three fewer parameters. The overdispersion factor (OD) decreases from 1.054 to 1.048 so the variance is only 5% larger than in a purely random distribution.
To check the statistical significance of periodicity in 2012-2019, the deviance obtained with Model (2.4) is compared with the deviance resulting from a regression without the periodic terms. With deviance = 652.0 (df = 622) and deviance = 663.1 (df = 626) respectively, the effect of periodicity in 2012-2019 is statistically significant (p = 0.032, F-test). Table 6 shows the numbers of excess deaths in areas A and B before and after 2018. In 2012-2017, the ratio of the increases in areas A and B is 2.11 (10.9% vs. 5.2%) whereas, in 2018-2019, there is a large 24% increase in Area A, but a negligible 1% increase in Area B.
The estimated total number of excess perinatal deaths in 2012-2019 with Model (2.3) is 567, somewhat less than resulting from a regression with Model (2.2) with 596 excess cases. As with the data from the study-and control region, the standardized residuals are checked for autocorrelation. For time lags of 1 to 11 months, no significant autocorrelation can be seen either in Area A or in Area B (see S5 Fig in S1 File). Eventually, the method used by Scherb et al. in [6] is applied to the present data from the study region. Linear logistic regression was used in [2]

Discussion
The present study uses a larger study region (12 prefectures   In addition to periodic peaks in April 2012-2019, a peak in perinatal mortality was observed in April 2011, followed by a trough in October 2011. Another significant increase was detected in November 2011. Fitting these effects by bell-shaped excess terms was intended to eliminate distortion of the long-term trend of the data. Possible explanations for the observed effects are suggested below. 1. The peak in April 2011 might be an immediate effect of the earthquake and tsunami. To check this hypothesis, regressions of perinatal mortality rates in the four prefectures of Area A are conducted with a dummy variable for March 2011. The increases in March 2011 are greater in the three coastal prefectures (Miyagi: +85%, Fukushima: 84%, and Ibaraki: +33%) than in the inland prefecture Tochigi (+20%, p = 0.61). The trends of the data in the four prefectures are presented in S11 and S12 Figs in S1 File.
2. The observation of a drop in perinatal mortality rates in October 2011 is unexpected; one would rather expect to see an increase in perinatal mortality seven months after the Fukushima accident. However, the radiation burst in March 2011 might have caused selective spontaneous abortions of the most radiosensitive embryos, resulting in a more radioresistant surviving fraction. Interestingly, the ratio of the effects in October 2011 in Area A (β 26 = -0.321) and Area B (β 27 = -0.128) is the same (2.5) as the ratio of the level shifts in Area A (β 7 = 0.251) and Area B (β 8 = 0.100), suggesting a common cause (e.g. radiation exposure) for these effects. In Japan, spontaneous abortions are only registered after 14 weeks of pregnancy, so the only way to find out about earlier abortions is to look at the trend of live births. In fact, Japanese monthly data of LB show a significant 2.7% drop in October 2011 (p = 0.006). The birth deficit is greater for males (-3.9%, p = 0.003) than for females (-1.5%, p = 0.25) which leads to a significant (p = 0.008) drop in the sex ratio at birth (male/female births) in Oct 2011, see S13-S15 Figs in S1 File. A similar drop in sex ratio was found in Czechia in November 1986, 7 months after the Chernobyl accident [11]. Greater radio-sensitivity of male than female embryos could explain the distorted sex ratio.
3. The peak in the study region in November 2012 is consistent with a peak in perinatal mortality in November 1987 in Germany after the Chernobyl accident (see S16 Fig and S10 Table in S1 File). The November 1987 peak was associated with an increase in cesium burden in pregnant females during the winter of 1986/87 when cows were fed contaminated silage from the summer of 1986 [3]. In Germany, the main path for cesium was the consumption of cow milk and dairy products. In Japan, the average per capita consumption of cow milk is substantially lower than in Germany (36 vs. 91 kg per year in 2002 [12]), but milk consumption may be greater in urban regions than in rural areas. This could explain why the November 2012 peak is smaller in Area A than in Area B which includes Tokyo.

Strengths and limitations
This study uses a larger study area (12 prefectures) and a longer study period (2002-2019) than previous studies of perinatal mortality after Fukushima [5][6][7][8]. In [8], the effect size in Fukushima prefecture was compared with that in four adjacent prefectures whereas the present study uses four core prefectures and eight surrounding prefectures for comparison. Due to much larger case numbers, any dependency of the effect on distance (as a proxy of dose) from the FDNPP can be determined more precisely. The remaining 35 Japanese prefectures are used as a control region. Thus, the study includes all perinatal deaths in Japan during 18 years; it has the greatest possible statistical power. An initial confirmatory analysis uses a predefined regression model while the subsequent exploratory analysis controls for observed short-term deviations from the regression model. Possible explanations for these deviations are suggested.
The main limitation of the study is its ecological design; the data are highly aggregated, so causal interpretations of the results are not possible. Interpretation of the observed increase in perinatal mortality as a radiation effect is complicated by the fact that, according to conventional wisdom, teratogenic effects are not expected to occur below a threshold dose of 100 mSv (see e.g. ICRP Publication 90 [13]). According to the UNSCEAR 2013 report [1], Annex A, Table 5, the mean estimated radiation dose to adults in the first year after the Fukushima accident was 1.0-4.3 mSv in Fukushima prefecture and 0.2-1.4 mSv in six neighboring prefectures (Group 3 prefectures).
To conclude, the present study finds a significant increase in perinatal mortality after the Fukushima accident in 12 prefectures of Japan surrounding the FDNPP. The increase is characterized by periodic peaks with maxima in April 2012-2019. In light of the massive increase in 2018 and 2019, continued investigation of perinatal mortality in contaminated regions of Japan is recommended.