Study area
Alongside other countries, over the past 100 years, significant loss of natural and semi-natural habitat has occurred in Japan, mainly because of urbanization and agricultural expansion in the course of Japan’s demographic and economic development (Supplementary Fig. 6).
Our research revealed aspects of biodiversity trends relating to depopulation, land use and climate change in WAPU locations across the Japanese archipelago. We focused on these areas because they represent regions of greatest interaction between nature and humans; they are where the greatest human population decline has been occurring in recent decades under national-scale depopulation. Moreover, these regions bear important geographical and ecological similarities to equivalent landscapes in South Korea, Northeast China and beyond, where wet-rice agriculture and its associated human livelihood practices are an embedded feature of local ecosystems.
Macro-level demographic data
Japan is already well known as the first Asian country to encounter long-term below-replacement human fertility, leading to demographic ageing and depopulation. To demonstrate Japan’s role as the DVC for Northeast Asia, we set it within both the world and regional contexts using the United Nations Population Division World Population Prospects 20242. While noting this dataset’s limitations for future fertility projections16, it presents globally comparable and actual national population change to 2023, and models projections from 2024 onwards. Our research measures only changes among selected biodiversity assemblages against changes in human population size, without considering the impact of changes in the population-age structure. We assume that ageing will have a similar, perhaps more muted, effect because of the tendency for older farmers to reduce productivity in labour-intensive agriculture68,69.
Adopting strict inclusion criteria, we narrowed the dataset pool of 233 countries and territories to 11 independent countries showing the earliest onset for indicators of long-term continuous depopulation caused primarily by natural causes (Supplementary Table 1). Among these 11, Japan was the only Asian country. Next, to demonstrate Japan’s position at the vanguard of a region-wide trend of low fertility leading to depopulation, we used the same indicators to compare Japan’s situation in relation to its Northeast and Southeast Asian neighbours. Supplementary Table 2 shows that the whole East Asian region is in the grip of a long-term depopulation trend, with Japan at the vanguard. The table also shows Northeast and Southeast Asian countries as distinct groups, with Northeast Asia depopulating and Southeast Asia growing but with decreasing human fertility leading to depopulation in the future.
Local-level population, land use and climate change data
To account for temporal and local spatial inconsistencies between our demographic and biodiversity data, we combined local-area mesh statistics with the National Census for 1995–2020. Mesh statistics divide a region into local-area grids (meshes) based on latitude and longitude and organize statistical data for each grid. Combining these datasets enables precise observation of demographic and biodiversity change over smaller geographical scales.
First, we aggregated human population data for the areas surrounding the study sites. After determining the latitude and longitude coordinates for each research site, we used R to identify the 1-km mesh code. We selected two scales, 25 km2 (5 × 5 km) and 100 km2 (10 × 10 km), as the basic unit range for the analysis. The process for obtaining data within a 25-km2 range involved identifying the 1-km mesh code, using ArcGIS to open the grid map, taking the confirmed 1-km mesh as the centre, and then counting two grids on the top, bottom, left and right to determine a square area with a side length of 5 km containing twenty-five 1-km grids. We obtained the population size of each 25-km2 range by aggregating the population size of these 25 grids. Similarly, for the 100-km2 range, after confirming the mesh code, we counted four grids to the left and upper sides of the confirmed 1-km grid, and five grids on its right and lower sides, to determine a square area with a side length of 10 km containing one hundred 1-km grids. Finally, we aggregated the population of each grid to obtain the population size within a 100-km2 range.
Next, we included data on land use change from 1997 to 2016, the closest available period. We applied the same scoping method as above, then combined it with national land numerical information data. National land numerical information statistics include data on the area of each land use category (paddy field, forest, wasteland, built site, highway site, river area, lake). Change in land use was divided for later analysis into three groups: urban; forested; and land under cultivation. The proportion of land under cultivation surrounding each survey site in 1997 averaged 21.9% with an s.e. of 1.2%, while forests averaged 42.9% (s.e. = 2.2%) and artificial land averaged 28.0% (s.e. = 1.9%). In 2016, the proportion of land under cultivation decreased to an average of 17.7% (s.e. = 1.2%), while forests increased slightly to 44.4% (s.e. = 2.2%) and artificial land expanded to 31.4% (s.e. = 2.1%).
We then compared the demographic and data on land use and found that they were highly significantly correlated between the 25-km2 and 100-km2 areas (population demography, r = 0.67, P < 0.001; urbanization, r = 0.95, P < 0.001; land under cultivation, r = 0.88, P < 0.001; forestation, r = 0.94, P < 0.001). Considering human movement distances, we used the 100-km2 (10 × 10 km) data in the subsequent analyses.
Finally, when analysing long-term biodiversity change, the impact of climate change cannot be ignored. As a covariate factor, we consequently included annual average temperatures from the Japan Meteorological Agency from 1995 to 2020, using data from the nearest meteorological station to each survey site.
Biodiversity data
We used unpublished biodiversity data for the period spanning 2004 to 2021 from Monitoring Sites 1000 (Appendix E in the Supplementary Information). The data used are from surveys conducted in WAPU landscapes. These landscapes account for approximately 40% of the land area in Japan70. The study sites were selected from agricultural secondary natural areas on land, excluding locations of newly created habitat by humans, such as artificial biotopes. Furthermore, we did not include ‘pristine’ natural environments free from obvious signs of human disturbance in our study. Over the past 50 years, from 1970 through to the present, Japan’s forest cover has remained approximately 67%, with around 40% consisting of plantations, such as Japanese cedar (Cryptomeria japonica) and Japanese cypress (Chamaecyparis obtusa) for timber production. The WAPU areas surveyed in our study include both succession forests and plantations. However, as already stated, our study sites are settled rural and peri-urban areas where agricultural land, forests and artificial land are intermingled in a mosaic pattern.
The survey area for each site was set to encompass a contiguous ecosystem, with a target size of approximately 30–100 hectares. Data from sites with a survey duration of less than 4 years were excluded; analysis was conducted on 158 sites within the latitude and longitude ranges of 24–44° and 123–145°, respectively, and with surveys conducted continuously for over 5 years. Detailed data for each site can be found in Supplementary Table 4. While the survey methods were generally standardized across all sites, we accounted for climate and topography variations at each survey location. Consequently, surveys differed in terms of survey area coverage and the number of surveys conducted annually for each individual site. Hence, we analysed individual sites to understand the cumulative fluctuating trends in biodiversity among selected species across Japan as a whole.
For this study, we focused on reliable information obtained through censuses for specific taxa, including birds (community data), frogs (population data for three species), insects (population data for two species of fireflies and community data for butterfly species) and plants (community data for both native and non-native species). Among the Monitoring Site 1000 data, camera trap data for mammals were included. However, because of changes in camera trap methods over time, direct comparisons with past data are challenging. To avoid drawing incorrect conclusions, we excluded these data from our analysis. Additionally, surveys of habitat area for the harvest mouse were conducted but abundance was not measured for this species. Hence, we excluded it from our analysis. The selected taxa cover a broad range of the ecological pyramid in semi-natural ecosystems, from top predators (birds) to producers (plants). Consequently, they provide a suitable dataset for evaluating the impacts of demographic and environmental change on organisms in semi-natural ecosystems (Supplementary Table 3).
During the breeding (April–August) and wintering (October–March) seasons, observers recorded the names and number of species, and the number of individuals, of birds observed within a 50-m radius along survey routes for subsequent analysis. Surveys were conducted three times per season, totalling six annually; the same routes were traversed to minimize survey biases for each round.
We conducted surveys of frogs, a key indicator species for both terrestrial and aquatic ecosystems, during the breeding season (October–June). Each survey location had a markedly different latitude, so the breeding season for each site was surveyed over a wide range. We counted the total number of egg masses for three frog species (R. japonica, R. ornativentris and R. pirica) because counting egg masses made the survey easier. Adult frogs inhabit areas such as farmland and mountain lowlands, making it preferable to focus on surveying their congregations during the breeding season to quantify their population. Moreover, as R. pirica inhabits a limited number of regions, we analysed it alongside R. ornativentris, which is phylogenetically close71.
Among insects we surveyed butterflies, whose larval host plants and adult nectar sources are located in semi-natural grasslands and forests around farmlands. Additionally, fireflies live close to water resources in the agricultural environment. For butterflies, we recorded the species name and number of individuals observed within a 5-m radius along the survey routes, with surveys conducted one or more times per month. For fireflies, we recorded the number of individuals during the breeding season for two key species, Genji fireflies (N. cruciata) and Heike fireflies (A. lateralis), which serve as indicators for aquatic organisms.
For plant species, once a month we recorded the names of plants with reproductive organs, such as flowers and fruits, observed along the survey routes. Each survey route was set to include a variety of landscape types, such as forests, forest edges, rice fields, fallow fields, wetlands, grasslands and roads; we recorded plant species observed along these routes. The surveyed plants included herbaceous species (seed plants and ferns) and woody species. Assessing plant abundance was not feasible, so we opted to analyse the number of species instead.
For taxonomic groups with available data on individual counts, we conducted analyses of annual variation in abundance using the method described below. Additionally, for birds, butterflies and plants with data on species numbers, we also conducted analyses of annual variation in species richness. We categorized plants into native and non-native species because the impact of population and land use changes may vary for each category. Separate analyses were conducted for the number of species for each group. For birds and butterflies, we analysed only native species because the number of non-native species was seven and one, respectively.
Analysis overview
Our analysis consisted of three components. First, we analysed and mapped human population transitions, change in land use and average temperature change (158 sites). Then, we analysed and mapped fluctuations in species abundance and richness (Supplementary Table 4). Finally, we analysed the relationships between population and change in land use, and changes in species abundance and richness (Supplementary Table 4).
Fluctuation trends in environmental factors
To examine annual change in environmental factors at each site, we used linear models with the population, land use or mean temperature as the response variables and year as the explanatory variable for each study site. In other words, the environmental factors in each sample Ei were defined as:
$${E}_{i} \approx N({\lambda }_{i})$$
Expectations under a Gaussian process λi were:
$${\lambda }_{i}={\beta }_{1}\times {Y}_{i}+{\beta }_{2}$$
where β1 and β2 are the parameters estimated using linear models and Yi is the year. Based on the significance of the slope for the year (α = 0.05), each study site was classified into three environmental change categories: increase (significantly positive); decrease (significantly negative); and no significant trend.
Fluctuating biodiversity trends
To examine annual changes in biodiversity among our selected taxonomic species at each site, we used generalized linear models (GLMs) (family = Poisson; link = log) with abundance or species richness as the response variables and year and month as the explanatory variables for each study site. The exception was where the survey date was included in the bird model as a random effect (generalized linear mixed models (GLMMs)) because bird surveys were performed multiple times per day. In other words, species richness or abundance in each sample Ni were defined as:
$${N}_{i} \approx {\mathrm{Poisson}}({\lambda }_{i})$$
In this case, expectations under a Poisson process λi were:
$${\lambda }_{i}=\left\{\begin{array}{l}\exp \left({\beta }_{1}\times {Y}_{i}+{\beta }_{2}+ {M}_{i}+{d}_{i}\right),\quad{\mathrm{if}}\; {\mathrm{birds}}\\ {{ }}\exp \left({\beta }_{1}\times {Y}_{i}+{\beta }_{2}+ {M}_{i}\right),\,\quad\qquad{\mathrm{other}}\; {\mathrm{organisms}}\end{array}\right.$$
where β1 and β2 are the parameters estimated using GLMs or GLMMs, Yi is the year and di is the random effect of the survey date. The effects Mi of the month on species richness or abundance were:
$${M}_{i}={\beta }_{3}\times {m}_{1,i}+{\beta }_{4}\times {m}_{2,i}+\cdots +{\beta }_{n+1}\times {m}_{n-1,i}$$
where β3, β4 …βn+1 were the parameters estimated using GLMs or GLMMs, and m1,i, m2,i…mn−1,i are the dummy variables of n month.
We excluded from the analysis some sites where the target taxa were not observed in any individuals. Based on the significance of the slope for the year (α = 0.05), each study site was classified into three change categories: increase (significantly positive); decrease (significantly negative); and no significant trend.
Effects of human anthropogenic factors on biodiversity
We evaluated how the population and land use change categorizations (increase, decrease, no significant trend) affect biodiversity. The response variable was the abundance or species richness of each taxon. The explanatory variables were: year, month; changes to population, urbanization, land under cultivation, forest and temperature; and interactions between year and each environmental change. Site identity and spatial autocorrelation were incorporated as random effects in GLMMs (family = Poisson, link = log). The spatial autocorrelation matrix was fitted using the spaMM package72 with the Matérn structure of pairwise correlations between the coordinates of each site location, accounting for spatial dependencies that might otherwise bias the results. In other words, species richness or abundance observed in each sample i was defined as:
$${N}_{i} \approx {\mathrm{Poisson}}({\lambda }_{i}).$$
In this case, expectations under a Poisson process λi were:
$$\begin{array}{l}{\lambda }_{i}=\exp \left(\left({\beta }_{1}\times {P}_{{{\mathrm{inc}}},i}+{\beta }_{2}\times {P}_{{{\mathrm{dec}}},i}+{\beta }_{3}\times {U}_{{{\mathrm{inc}}},i}+{\beta }_{4}\times {C}_{{{\mathrm{dec}}},i}+{\beta }_{5}\times {F}_{{{\mathrm{dec}}},i}\right.\right.\\\qquad\left.+{\beta }_{6}\times {F}_{{{\mathrm{inc}}},i}+{\beta }_{7}\times {T}_{{{\mathrm{inc}}},i}+{\beta }_{8}\right)\times {Y}_{i}+{\beta }_{9}\times {P}_{{{\mathrm{inc}}},i}\\\qquad+{\beta }_{10}\times {P}_{{{\mathrm{dec}}},i}+{\beta }_{11}\times {U}_{{{\mathrm{inc}}},i}+{\beta }_{12}\times {C}_{{{\mathrm{dec}}},i}+{\beta }_{13}\times {F}_{{{\mathrm{dec}}},i}\\\qquad\left.+{\beta }_{14}\times {F}_{{{\mathrm{inc}}},i}+{\beta }_{15}\times {T}_{{{\mathrm{inc}}},i}+{\beta }_{16}+{M}_{i}+{s}_{i}+{\xi }_{i}\right)\end{array}$$
where β1, β2, … β16 are parameters estimated using GLMMs, Yi is the year, Pinc,i and Pdec,i are the dummy variables of the population change categories, Uinc,i is the dummy variable of the urbanization change categories, Cdec,i is the dummy variable of the land under cultivation change categories, Fdec,i and Finc,i are the dummy variables of the forest change categories, Tinc,i is the dummy variable of the temperature change categories and si is the random effect of study site (‘inc’ = increase, ‘dec’ = decrease in all variable labels). ξi is the random effect accounting for spatial autocorrelation. Also, the effects of the month (Mi) on species richness or abundance were:
$${M}_{i}={\beta }_{17}\times {m}_{1,i}+{\beta }_{18}\times {m}_{2,i}+\cdots +{\beta }_{n+15}\times {m}_{n-1,i}$$
where β17, β18, … βn + 15 are the parameters estimated using GLMMs, and m1,i, m2,i…mn-1,i are the dummy variables of n month.
Before the analysis, we calculated the variance inflation factor values for each model to check for issues of multicollinearity (all variance inflation factors were less than 3.3). The best model was selected based on the minimum Akaike information criterion by evaluating all possible combinations of the explanatory variables. We conducted all analyses in R v.4.4.2 (ref. 73) using the fitme function with the spaMM package72.
Reporting summary
Further information on research design is available in the Nature Portfolio Reporting Summary linked to this article.
