The dynamics of the ground radiation dose
The ground radiation dose at the collection sites in the 7 localities was measured at the time of butterfly sample collection in 6 sampling attempts over 3 years. The ground radiation dose reached its maximum level in the spring of 2011 and declined afterwards (Figure 2a, b), likely due to fast-decaying radionuclides. However, even in the fall of 2013, considerable levels were detected, especially from Fukushima and Motomiya. The level of Hirono was also high. It is important to note that the ground radiation doses in these localities in the fall of 2013 were still higher than the doses in other localities in the spring of 2011.
The aAR dynamics of the field-collected adults
We monitored the aAR of field-collected adults from the 7 localities in 6 sampling attempts over 3 years. We first plotted the raw aAR data from the 7 localities on a single graph (Figure 3a). Considering that the normal level of aAR is below 10% (see below and also [12]), the aAR in the spring of 2011 was higher than the normal level in Hirono, Fukushima, and Mito. An increase in the fall of 2011 was then observed in many localities except Tsukuba and Mito. Notably, the localities with high levels of pollution, i.e., Fukushima and Motomiya, showed a clear peak in the fall of 2011, approaching 40%. The localities with intermediate pollution levels, i.e., Iwaki, Hirono, and Takahagi, showed similarly high aARs both in the fall of 2011 and in the spring of 2012. In contrast, the locality with the lowest contamination level, Tsukuba, peaked later, in the spring of 2012. All localities showed normal aAR levels in the fall of 2012 and afterwards. However, we noticed that Mito behaved differently from the rest, at least during the spring and fall of 2011, which was apparent in the cluster analysis (Figure 3b) and the multidimensional scaling plot (Additional file 2: Figure S1). Based on these results, we excluded Mito as an outlier from the subsequent analyses.
After excluding Mito, we found statistically significant differences between the peak time-point (i.e., the fall of 2011) and the latest time-point (i.e., the fall of 2013) (p = 0.0011; Scheffe test), supporting the visual interpretations of the dynamics discussed above (Figure 3a). The cluster analysis revealed that the 6 localities were segregated into two clusters: the Fukushima-Motomiya cluster (FM cluster) and the cluster including Hirono, Takahagi, Iwaki, and Tsukuba (HTIT cluster) (Figure 3b). The FM cluster corresponded to the inland cities relatively far from the ocean (but heavily contaminated) to the northwest of the FNPP, and the HTIT cluster corresponded to the coastal cities (but not affected by Tsunami) along the Pacific Ocean to the south of the FNPP. The cluster analysis also revealed two time-point clusters: the cluster containing the fall of 2011 and the spring of 2012 and the cluster containing the rest of the time points (Figure 3b), suggesting that major biological changes occurred during the fall of 2011 and the spring of 2012 in Fukushima and Motomiya.
Based on the clustering results, we then calculated the aAR values for the FM and HTIT clusters of all the combined samples from the included localities, which were plotted against time (Figure 3c). As expected, the FM cluster showed a higher aAR than the HTIT cluster except in the fall of 2012. The FM cluster showed a steep decrease in the spring of 2012 and also in the fall of 2012, but the HTIT cluster showed a large decrease only in the fall of 2012, after which the aARs were at normal levels in both clusters. The samples from other localities (excluding the major 7 localities) that were minimally contaminated (Additional file 1: Tables S1-S6) were plotted together with the FM and HTIT clusters (Figure 3c). Using these 3 groups, statistically significant differences were found among different time points in multiple comparisons (p = 0.037; Freidman test). Because these other localities were at the minimal contamination level, we found that the aARs of the normal populations were largely below 10%, ranging from 0% to 13.5% with a median of 6.4%. Similar results have been shown in a previous study [12], where the aARs of the field-collected samples had a mean value of 8.2%.
The adult abnormalities were classified into four different categories (wing shape, wing color patterns, appendages, and other parts), and their changes over time were examined in the 6 different localities (Figure 4a-f). Overall, it seemed that abnormalities in the appendages occurred earlier, followed by abnormalities in wing color patterns and wing shape, especially in the relatively highly contaminated localities. The relatively small level of wing-shape abnormalities may be due to the limited mobility of such individuals, resulting in a lower probability of being captured. A cluster analysis on the wing shape abnormalities revealed two clusters: one cluster included Fukushima, Motomiya, and Hirono (FMH cluster), and the other cluster included Iwaki, Takahagi, and Tsukuba (ITT cluster) (Figure 4g). The incidence of abnormality in the FMH cluster peaked in the fall of 2011, whereas abnormality in the ITT cluster peaked in the spring of 2012 (Figure 4h), showing a delay of peaking in the low-contamination localities.
We then plotted the four categories of the abnormalities individually throughout the 3-year period. The incidence of wing shape abnormality was irregular and sporadic over time (Figure 4i). Abnormalities in wing color patterns peaked mainly in the spring of 2012 (Figure 4j), and abnormalities in the appendages peaked mainly in the fall of 2011 (Figure 4k). The abnormalities of the other parts were quite irregular and sporadic (not shown). Wing shape abnormalities may be behaviorally lethal or sublethal, and such abnormal individuals cannot be captured frequently by us in the field, resulting in the sporadic and low-level occurrence. In contrast, abnormalities in appendages and color patterns may be less lethal, and individuals with such abnormalities can still be captured by us in the field, resulting in the peaks at specific time-points. This interpretation is in good agreement with the results of the F1 generation discussed later.
Consistent with these results, the number of individuals caught per minute in the field was low in both the fall of 2011 and in the spring of 2012 (Figure 5a, b). These dynamics should be interpreted with caution because it is likely that the adult density is generally low in the spring and gradually increases toward the summer every year, at least in Tokyo [16] and most likely in the entire Tohoku and Kanto districts. Moreover, the number of individuals caught per minute is also affected by environmental factors including the weather of the day. Nevertheless, the relatively low number of individuals in the fall of 2011 was notable compared with the number in the fall of 2013 (p = 0.041; Scheffe test). For example, in Hirono, we were able to collect only very small number of individuals in the fall of 2011 and in the spring of 2012 despite the fine weather.
In summary, the radiation levels peaked first in the spring of 2011, followed by an increase in the aARs and a decrease in the number of individuals in the fall of 2011. Then, the aAR decreased and the number of individuals increased back to normal levels. This chronological order supports the normalization hypothesis, rejecting the extended mortality hypothesis and the extinction hypothesis (see Introduction).
The tAR and aAR dynamics of the reared offspring
We next monitored the total AR (tAR) of the offspring (F1) generation from the adults caught in the 7 localities (Figure 6a; Additional file 3: Table S7). The tAR includes the deaths in the larval, prepupal, and pupal stages and the abnormalities in the surviving adults, which could be examined by rearing the offspring in the laboratory. The tARs peaked mostly either in the fall of 2011 or in the spring of 2012, similar to the results of the aAR of the parent generation shown above, with a significant difference between the fall of 2011 and the spring of 2013 (p < 0.0001; Scheffe test with Benjamini-Hochberg correction [17]) and between other two time-points. However, we did not obtain clear clusters in the cluster analysis and the multi-dimensional scaling (not shown), suggesting that the tAR patterns over time were similar among the localities monitored. The decrease of the tAR was relatively slow: the tAR levels decreased back to normal in the spring of 2013, taking more time than the aAR of the parent generation. Another noteworthy point was that the levels of the tARs were rather high. These results suggest a transgenerational character of abnormalities, at least from the P to F1 generations, throughout this time period and also suggest that selection against abnormal individuals occurred in the field throughout this time period.
We also plotted the aAR of the F1 generation (Figure 6b). The aARs of the F1 generation had higher values and peaked more sharply in the fall of 2011 than the aAR of the parent generation in all seven localities. For comparison, we integrated all 5 localities (excluding Mito and Tsukuba due to a lack of monitoring data in the fall of 2011) by considering that all samples were collected from a single Tohoku-Kanto area, and we plotted the aAR of the P and F1 generations and the tAR of the F1 generation together (Figure 6c). All 3 ARs showed an overall similar pattern: a peak in the fall of 2011 and a subsequent decrease. However, it is notable that the tAR peak was much higher than the peaks of aARs and that the tAR peak in the fall of 2011 extended to the spring of 2012. The aARs of the P and F1 generations were similar to each other in that they both showed a single peak in the fall of 2011, but they were clearly different in magnitude; the aAR of the F1 generation was much higher. These results confirmed the idea that the wild populations of this butterfly produced many abnormal individuals that were selected out before appearing as viable adults in the field.
As in the previous analysis of the parent generation, the adult abnormalities were classified into 4 categories (wing shape, wing color patterns, appendages, and other parts), and the changes in the aARs over time were examined (Additional file 4: Figure S2). The levels of aARs were especially high for the wing shape in the fall of 2011 (Additional file 4: Figure S2a) compared with the wing color patterns (Additional file 4: Figure S2b), the appendages (Additional file 4: Figure S2c), and the other parts (not shown).
The dynamics of the mortality rate (MR) added further insights into what occurred in the polluted area. When all 5 localities were combined as samples from a single collection area, the MR throughout immature life stages (i.e., larval, prepupal, and pupal stages) showed the highest in the spring of 2012, scoring approximately 60%, and decreased afterwards to the normal level in the spring of 2013 (Figure 7a). When 5 localities were examined separately, Hirono notably scored more than 90% in the fall of 2011 (Figure 7b), suggesting that the large-scale death occurred in Hirono at that time. Takahagi, Mito, Fukushima, and Iwaki scored more than 50% in the spring of 2012.
When the MRs of developing life stages (the larval, prepupal, and pupal MRs) were examined separately in all 5 localities combined, the pupal MR was highest in the spring of 2012 and decreased sharply afterwards (Figure 7c). When the stage-specific MRs were simply divided by the number of days for each stage, prepupal death was notable in the spring of 2011 (Figure 7d). The dying stage may be different between 2011 and 2012. The main selection stage was likely the prepupal stage in the spring of 2011 and then moved to the pupal and larval stages afterwards. Because we reared these larvae, prepupae, and pupae under standard laboratory conditions, where no artificial radionuclides were detected, and in the absence of any predators, the selection pressures against these immature stages were likely to be more severe in the field than in the laboratory.
We then examined the dynamics of the stage-specific MRs in 7 different localities (Additional file 5: Figure S3). Most locations showed peak larval mortality in the spring of 2012, with the exception of Hirono, at which larval mortality peaked in the fall of 2011, and of Motomiya, at which larval mortality peaked in the fall of 2012 (Additional file 5: Figure S3a). At the prepupal stage, the highest level of mortality was found in the spring of 2011 in most localities, again with the exception of Hirono, at which prepupal mortality peaked in the fall of 2011 (Additional file 5: Figure S3b). At the pupal stage, most localities exhibited peak mortality in the spring of 2012, again with the exception of Hirono, at which mortality peaked in the fall of 2011 (Additional file 5: Figure S3c).
Correlations of ARs with the ground radiation level and with the distance
To investigate the possible causal involvement of radioactive pollution released by the Fukushima nuclear accident in the AR dynamics in 2011–2013, we examined correlations of ARs with the ground radiation dose at the collection sites and with the distance from the FNPP, although a similar analysis has already been performed for the 2011 data in previous papers [9,12,14]. Radioactive materials released early, which mainly contained short-lived radionuclides, were dispersed in a concentric fashion, whereas radioactive materials released later, which mainly contained cesium, were dispersed mainly toward the northwestern cities [2]. Therefore, our categorization may help to distinguish the effects.
In the P generation, changes in Pearson correlation coefficients over the 6 time points between the aAR level and the ground radiation level and between the aAR level and the distance were similar to each other, but with a slight delay in the distance (Figure 8a). Correlation was strongest in the fall of 2011 for the ground radiation and in the spring of 2012 for the distance, both of which then decreased toward the fall of 2012 or toward the spring of 2013.
Unexpectedly, the coefficients then increased again in the spring of 2013 for the ground radiation and in the fall of 2013 for the distance. To understand what occurred, scatter plots were examined (Additional file 6: Figure S4). The relatively high correlation of aAR with the ground radiation dose in the fall of 2011 was disrupted mainly because of the fast decrease of the aARs in the high-radiation localities (Additional file 6: Figure S4a). By the fall of 2013, the aARs in the low-radiation localities decreased, recreating relatively high correlation again, but with much lower aARs. Similar dynamics, although less clear, were observed in the distance correlation (Additional file 6: Figure S4b).
To understand the latest state of the aARs, we examined the aARs of the adult samples that we collected in the fall of 2013 from 12 localities (Figure 8b). The aARs were relatively low in most localities but relatively high in Motomiya and Koriyama and very high in Iitate, confirming the origin of the relatively high correlation in the fall of 2013.
Similar dynamics were obtained in the F1 generation with the aARs (Figure 8c) and with the tARs (Figure 8d), confirming the trends obtained in the P generation. Scatter plots also showed a similar trend (Additional files 7, 8: Figure S5, S6). However, the peak time difference, observed in the P generation, between the ground radiation dose and the distance (Figure 8a) was not observed in the F1 generation (Figure 8c, d).