Ecology

Calculating Potential Evapotranspiration using Penman-MonteithAmong several equations used to estimate PET, an implementation of the Penman-Monteith equation originally presented by the Food and Agriculture Organization FAO-561, is considered a standard method3,12,13,49. FAO-561 defined PET as the ET of a reference crop (ET0) under optimal conditions, in this case with the specific characteristics of well-watered grass with an assumed height of 12 centimeters, a fixed surface resistance of 70 seconds per meter and an albedo of 0.231. Less specifically, “reference evapotranspiration”, generally referred to as “ET0”, measures the rate at which readily available soil water is evaporated from specified vegetated surfaces2,13, i.e., from a uniform surface of dense, actively growing vegetation having specified height and surface resistance, not short of soil water, and representing an expanse of at least 100 m of the same or similar vegetations1,13. ET0 is one of the essential hydrological variables used in many research efforts, such as study of the hydrologic water balance, crop yield simulation, irrigation system management and in water resources management, allowing researchers and practitioners to study the evaporative demand of the atmosphere independent of crop type, crop development and management practices2,4,13,49. ET0 values measured or calculated at different locations or in different seasons are comparable as they refer to the ET from the same reference surface. The factors affecting ET0 are climatic parameters, and crop specific resistances coefficients solved for reference vegetation. Other crop specific coefficients (Kc) may then be used to determine the ET of specific crops (ETc), and which can in turn be determined from ET01.As the Penman-Monteith methodology is predominately a climatic approach, it can be applied globally as it does not require estimations of additional site-specific parameters. However, a major drawback of the Penman-Monteith method is its relatively high need for specific data for a variety of parameters (i.e., windspeed, relative humidity, solar radiation). Zomer et al.18 compared five methods of calculating PET with parameters from data available at the time and settled upon using a Modified Hargreaves-Thornton equation50 which required less parametrization to produce the Global-AI_PET_v116,17,18. Several other attempts to produce global PET datasets with concurrently available global datasets came to similar conclusions51,52,53. The Modified Hargreaves-Thornton method required less parameterization with relatively good results, relying on datasets which were available at the time for a globally applicable modeling effort. The Global-AI_PET_v1 used the WorldClim_v1.420 downscaled climate dataset (30 arcseconds; averaged over the period 1960–1990) for input into the global geospatial implementation of the Modified Hargreaves-Thornton equation, applied on a per grid cell basis at approximately 1 km resolution (30 arcseconds). More recently, the UK Climate Research Unit released the “CRU_TS Version 4.04”, which now includes a Penman-Monteith calculated PET (ET0) global coverage, however at a relatively coarse resolution of 0.5 × 0.5 degrees. A number of satellite-based remote sensing datasets22,54,55,56,57 are now available and in use to provide the parameters for ET0 estimates, in some cases providing high spatial and/or temporal resolution and are likely to become increasingly utilized as the historical data record lengthens and sensors improve.The latest 2.0 versions of WorldClim58 (currently version 2.1; released January 2020), in addition to being updated with improved data and analysis, and a revised baseline (1970–2000), includes several additional primary climatic variables, beyond temperature and precipitation, namely: solar radiation, wind speed and water vapor pressure. The addition of these variables allowed that the global data now available was sufficient to effectively parameterize the FAO-56 equation to estimate ET0 globally at the 30 arc seconds scale (~1 km at equator).The FAO-56 Penman-Monteith equation, described in detail below, has been implemented on a per grid cell basis at 30 arc seconds resolution, using the Python programming language (version 3.2). The data to parametrize the various components equations required to arrive at the ET0 estimate were obtained from the Worlclim 2.158 climatological dataset, which provides values averaged over the time period 1970–2000 for minimum, maximum and average temperature; solar radiation; wind speed, and water vapor pressure. Subroutines in the program include calculation of the psychrometric constant (aerodynamic resistance), saturation vapor pressure, vapor pressure deficit, slope of vapour pressure curve, air density at constant pressure, net shortwave radiation at crop surface, clear-sky solar radiation, net longwave radiation at crop surface, net radiation at the crop surface, and the calculation of daily and monthly ET0. This process is described below. Geospatial processing and analysis were done using ArcGIS Pro v 2.9 (ESRI, 2020), Python (ArcPy) programming language (version 3.2), and Microsoft Excel for further data analysis, graphics and presentation.Global Reference Evapotranspiration (Global-ET0)Penman59, in 1948, first combined the radiative energy balance with the aerodynamic mass transfer method and derived an equation to compute evaporation from an open water surface from standard climatological records of sunshine, temperature, humidity and wind speed. This combined approach eliminated the need for the parameter “most difficult” to measure, surface temperature, and allowed for the first time an opportunity to make theoretical estimates of ET from standard meteorological data. Consequently, these estimates could also now be made retrospectively. This so-called combination method was further developed by many researchers and extended to cropped surfaces by introducing resistance factors. Among the various derivations of the Penman equation is the inclusion of a bulk surface resistance term60, with the resulting equation now called the Penman-Monteith equation3, as standardized in FAO-561 and subsequently by the American Society of Civil Engineers – Technical Committee on Standardization of Reference Evapotranspiration12,13,49,61. The FAO-56 Penman-Monteith form of the combination equation to estimate ET0 is calculated as:$$ETo=frac{Delta left({R}_{n}-Gright)+{rho }_{a}{c}_{p}frac{({e}_{s}-{e}_{a})}{{r}_{a}}}{Delta +gamma left(1+frac{{r}_{s}}{{r}_{a}}right)}$$
(1)
WhereET0 is the evapotranspiration for reference crop, as mm day−1Rn is the net radiation at the crop surface, as MJ m−2 day−1G is the soil heat flux density, as MJ m−2 day−1cp is the specific heat of dry airpa is the air density at constant pressurees is the saturation vapour pressure, as kPaea is the actual vapour pressure, as kPaes – ea is the saturation vapour pressure deficit, as kPa(Delta ) is the slope vapour pressure curve, as kPa °C−1(gamma ) is the psychrometric constant, as kPa °C−1rs is the bulk surface resistance, as m s−1ra is the aerodynamic resistance, as m s−1Psychrometric Constant (γ)The Atmospheric Pressure (Pr, [KPa]) is the pressure exerted by the weight of the atmosphere and is thus dependent on elevation (elev, [m]). To a certain (and limited) extent evaporation is promoted at higher elevations:$$Pr=101.3ast {left(frac{293-0.0065ast elev}{293}right)}^{5.26}$$
(2)
Instead, the psychrometric constant, [γ, kPa C−1] is expressed as:$$gamma =frac{{c}_{p}ast Pr}{varepsilon ast lambda }=frac{0.001013ast Pr}{0.622ast 2.45}$$
(3)
Where cp is the specific heat at constant pressure [MJ kg−1 °C−1] and is equal to 1.013 10−3, λ is the latent heat of vaporization [MJ kg−1] and is equal to 2.45, while ε is the molecular weight ratio between water vapour and dry air and is equal to 0.622.Elevation data has been obtained from the Shuttle Radar Topography Mission (SRTM) aggregated to 30 arc-second spatial resolution62 and combined with the USGS GTOPO3063 database for the areas north of 60°N and south of 60°S where no SRTM data was available (available at https://worldclim.org).Air Density at Constant Pressure [ρa]The mean Air Density at Constant Pressure [ρa, Kg m−3] can be represented as:$${rho }_{a}=frac{Pr}{{T}_{Kv}ast R}$$
(4)
While R is the specific heat constant (0.287, KJ Kg−1 K−1), the virtual temperature TKv can be represented as well as:$${T}_{Kv}=1.01ast ({T}_{avg}+273)$$
(5)
With Tavg as the mean daily air temperature at 2 m height [C°].Saturation Vapor Pressure [KPa]Saturation Vapor Pressure [KPa] is strictly related to temperature values (T)$${e}_{s_T}=0.6108ast ex{p}^{left[frac{17.27ast T}{T+237.3}right]}$$
(6)
Values of saturation vapor pressures, as function of temperature, are calculated for both Minimum Temperature [Tmin, C°] and Maximum temperature [Tmax, C°]. Due to nonlinearity of the equation, the mean saturation vapour pressure [es, KPa] is calculated as the average of saturation vapour pressure at minimum [es_min] and maximum temperature [es_max]$${e}_{s}=frac{{e}_{s_Tmax}+{e}_{s_Tmin}}{2}$$
(7)
The actual vapour pressure [ea, KPa] is the vapour pressure exerted by the water in the air and is usually calculated as function of Relative Humidity [RH]. Water vapour pressure is already available as one of the Worldclim 2.1 variables.$${e}_{a}=RH/100,ast ,{e}_{s}$$
(8)
The vapour pressure deficit (es-ea), [KPa] is the difference between the saturation (es) and actual vapour pressure (({e}_{a})).Slope of Saturation Vapor Pressure (Δ)The Slope of Saturation Vapor Pressure [Δ, kPa C−1] at a given temperature is given as function of average temperature:$$Delta =frac{4098ast 0.6108,ex{p}^{left(frac{17.27ast {T}_{avg}}{{T}_{avg}+237.3}right)}}{{left({T}_{avg}+237.3right)}^{2}}$$
(9)
Where Tavg [C°] is the average temperature.Net Radiation At The Crop Surface (R
n)Net radiation [Rn, MJ m−2 day−1] is the difference between the net shortwave radiation [Rns, MJ m−2 day−1] and the net longwave radiation [Rnl, MJ m−2 day−1], and is calculated using solar radiation (Rs). In Worldclim 2.1 solar radiation (Rs) is given as KJ m−2 day−1. Thus, for computation of ET0, its unit should be converted to MJ m−2 day−1 and thus its value should be divided by 1000. The net accounting of either longwave and shortwave radiation sums up the incoming and outgoing components.$${R}_{n}={R}_{ns}-{R}_{nl}$$
(10)
The net shortwave radiation [Rns, MJ m−2 day−1] is the fraction of the solar radiation Rs that is not reflected from the surface. The fraction of the solar radiation reflected by the surface is known as the albedo [α]. For the green grass reference crop, α is assumed to have a value of 0.23. The value of Rns is:$${R}_{ns}={R}_{s},ast ,(1-alpha )$$
(11)
The difference between outgoing and incoming longwave radiation is called the net longwave radiation [Rnl]. As the outgoing longwave radiation is almost always greater than the incoming longwave radiation, Rnl represents an energy loss. Longwave energy emission is related to surface temperature following Stefan-Boltzmann law. Thus, longwave radiation emission is calculated as positive in the outward direction, while shortwave radiation is positive in the downward direction. The net energy flux leaving the earth’s surface is influenced as well by humidity and cloudiness$${R}_{nl}=sigma ast left(frac{{T}_{max,,K}^{4}+{T}_{min,,K}^{4}}{2}right)ast left(0.34-0.14ast sqrt{{e}_{a}}right)ast left(1.35ast frac{{R}_{s}}{{R}_{so}}-0.35right)$$
(12)
Where σ represent the Stefan-Boltzmann constant (4.903 10-9 MJ K−4 m−2 day−1), Tmax,K and Tmin,K the maximum and minimum absolute temperature (in Kelvin; K = C° + 273.16), ea is the actual vapour pressure; Rs the measured solar radiation [MJ m−2 day−1] and Rso is the calculated clear-sky radiation [MJ m−2 day−1]. Rso is calculated as function of extraterrestrial solar radiation [Ra, MJ m−2 day−1] and elevation (elev, m):$${R}_{so}={R}_{a}ast (0.75+0.00002ast elev)$$
(13)
The extraterrestrial radiation, [Ra, MJ m−2 day−1], is estimated from the solar constant, solar declination and day of the year. It requires specific information about latitude and Julian day to accomplish a trigonometric computation of the amount of solar radiation reaching the top of the atmosphere following trigonometric computations as shown in Allen et al.1.Although the soil heat flux is small compared to Rn, particularly when the surface is covered by vegetation, changes of soil heat flux may still be relevant at monthly scale. However, accurate assessments of soil heat flux may require computation of soil heat capacity, related to its mineral composition and water content, which in turn may be rather inaccurate at global scale at resolution of 30 arc sec. Thus, for simplicity, changes in soil heat fluxes are ignored (G = 0).Bulk Surface Resistance (r
s)The resistance nomenclature distinguishes between aerodynamic resistance and surface resistance factors. The surface resistance parameters are often combined into one parameter, the ‘bulk’ surface resistance parameter which operates in series with the aerodynamic resistance. The surface resistance, rs, describes the resistance of vapour flow through stomata openings, total leaf area and soil surface. The aerodynamic resistance, ra, describes the resistance from the vegetation upward and involves friction from air flowing over vegetative surfaces. Although the exchange process in a vegetation layer is too complex to be fully described by the two resistance factors, good correlations can be obtained between measured and calculated evapotranspiration rates, especially for a uniform grass reference surface.A general equation for the bulk surface resistance (rs, [s m−1]) describes a ratio between the bulk stomatal resistance of a well illuminated leaf (rl) and the active sunlit leaf area of the vegetation:$${r}_{s}=frac{{r}_{l}}{LA{I}_{active}}$$
(14)
The stomatal resistance of a single leaf under well-watered conditions has a value of about 100 s m−1. It can be assumed that about half (0.5) of the total LAI is actively contributing to vapour transfer, while it can also be roughly generalized that for short crops there is a linear relation between LAI and crop height (h):$$LAI=24ast h$$
(15)
When the evapotranspiration simulated with the Penman-Monteith method is referred to a specific reference crop, denoted as ET0, a simplified computation of the method can occur that defines a priori specific variables into constant values. In this case, the reference surface is a hypothetical grass reference crop, well-watered grass of uniform height, actively growing and completely shading the ground, with an assumed crop height of 0.12 m, and an albedo of 0.23. The surface resistance for this hypothetical grass can be simplified to the following:$${r}_{s}=frac{100}{0.5ast 24ast h}$$
(16)
For such reference crop the surface resistance is fixed to 70 s m−1 and implies a moderately dry soil surface resulting from about a weekly irrigation frequency.Aerodynamic Resistance (r
a)The aerodynamic resistance [s m−1] verifies the transfer of water vapour and heat from the vegetation surface into the air, and is controlled by both vegetation status but also atmospheric turbulence under theoretical aspect as:$${r}_{a}=frac{lnleft[frac{{z}_{m}-d}{{z}_{om}}right]ast lnleft[frac{{z}_{h}-d}{{z}_{oh}}right]}{{k}^{2}{u}_{z}}$$
(17)
Zm [m] is the height [h] of wind measurements and Zh [m] is the height of humidity measurements. These are normally set at 2 meters height, although several climate models may provide them for higher heights (e.g. 10 m). The zero plane displacement (d [m]) term can be estimated as two thirds of crop height, while Zom is the roughness length governing momentum transfer, and can be calculated as Zom = 0.123 * h.The roughness length governing transfer of heat and vapour, Zoh [m], can be approximated as one tenth of Zom. k is the von Karman’s constant, equal to 0.41, and uz [m s-1] is the wind speed at height z.The reference surface, as stated, is a hypothetical grass reference crop, well-watered grass of uniform height, actively growing and completely shading the ground, with an assumed crop height of 0.12 m, and an albedo of 0.23. For such reference crop the surface resistance is fixed to 70 s m-1 and implies a moderately dry soil surface resulting from about a weekly irrigation frequency.When crop height is equal to 0.12 and wind/humidity measurements are taken at 2 meters height, then the aerodynamic resistance can be simplified as:$${r}_{a}=frac{208}{{u}_{2}}$$
(18)
Reference Evapotranspiration (ET
0)Given the above, and the specific properties of the standard reference crop, the FAO-56 Penman-Monteith method to estimate ET0 then can be calculated as:$$ETo=frac{0.408ast Delta ast left({R}_{n}-Gright)+gamma frac{900}{{T}_{avg}+273}ast {u}_{2}ast left({e}_{s}-{e}_{a}right)}{Delta +gamma left(1+frac{{r}_{s}}{{r}_{a}}right)}$$
(19)
Aridity Index (AI)Aridity is often expressed as a generalized function of precipitation and PET. The ratio of precipitation over PET (or ET0). That is, the precipitation available in relation to atmospheric water demand64 quantifies water availability for plant growth after ET demand has been met, comparing incoming moisture totals with potential outgoing moisture65.Geospatial analysis and global mapping of the AI for the averaged 1970–2000 time period has been calculated on a per grid cell basis, as:$$Al=MA_Prec/MA_E{T}_{0}$$
(20)
where:AI = Aridity IndexMA_Prec = Mean Annual PrecipitationMA_ET0 = Mean Annual Reference EvapotranspirationMean annual precipitation (MA_Prec) values were obtained from the WorldClim v 2.158, as averaged over the period 1970–2000, while ET0 datasets estimated on a monthly average basis by the Global-ET0 (i.e., modeled using the method described above) were aggregated to mean annual values (MA_ET0). Using this formulation, AI values are unitless, increasing with more humid condition and decreasing with more arid conditions.As a general reference, a climate classification scheme for Aridity Index values provided by UNEP64 provides an insight into the climatic significance of the range of moisture availability conditions described by the AI.
Aridity Index Value

Climate Class

0.65

Humid More

Overall characteristics of air pollutantsThe results of previous studies indicated that local pollution is highly important in determining the emissions of air pollutant. Therefore, in this study, we estimated the changes in pollution and the AQI between the pre-COVID and COVID lockdown periods and among the different regions in Ji’nan. A comparison of the different pollutant concentrations analysed in this study shows that the concentrations of almost all pollutants decreased during the COVID lockdown period; only the concentration of O3 increased continuously as the COVID lockdown period progressed (Fig. 1).Figure 1Spatial distributions of the different observation sites and industrial enterprises above a designated size threshold in Ji’nan city. JCE, machine tool factory No. 2; LSX, technical college; JNS, Ji’nan fourth building group; KFQ, economic development zone; KGS, Kegansuo; LWZ, Laiwu memorial hall; NKS, Agricultural Scientific Institute; SZZ, Seed warehouse of Shandong Province; SJC, Ji’nan monitoring station; TXG, Taixing company; CQD, Changqing school. Red circles, red triangles and red squares represent stations in urban, urban-industrial and suburban regions, respectively. The map of Observation site was completed by the geostatistical analysis module of ArcGIS (version 10.3, https://developers.arcgis.com/).Full size imageDuring the observation period, the daily average mass concentrations of PM10, PM2.5, SO2, NO2, CO, and O3 in Ji’nan were 137.09 µg/m3, 101.35 µg/m3, 22.70 µg/m3, 39.77 µg/m3, 1.28 mg/m3, and 71.84 µg/m3, respectively (Fig. 2). The mass concentrations of PM10 and PM2.5 exceeded the daily average Grade I values (50 µg/m3 and 35 µg/m3) of the Ambient Air Quality Standard of China (CAAQS, GB 3095-2012) during the whole observation period. In contrast, the mass concentrations of NO2, SO2, CO and O3 were substantially lower than the daily average Grade I values (80 µg/m3, 50 µg/m3, 4 mg/m3 and 100 µg/m3, respectively) of the CAAQS each day. During the pre-COVID period, the daily average mass concentrations of PM10, PM2.5, SO2, NO2, CO, and O3 in Ji’nan were 177.03 µg/m3, 125.94 µg/m3, 26.39 µg/m3, 54.52 µg/m3, 1.59 mg/m3, and 60.72 µg/m3, respectively. The mass concentrations of all these pollutants, except NO2, CO and O3, exceeded the daily average Grade I values of the CAAQS. The mass concentration trends during the COVID lockdown period were consistent with those during the pre-COVID period, but there were significant differences in the concentrations between the periods. In summary, the air quality in Ji’nan was generally good from January 24 to February 7, 2020, mainly due to the strict prevention and control measures for COVID-19.Figure 2Temporal variations in the mass concentrations of air pollutants (PM10, PM2.5, NO2, SO2, CO and O3) at the urban site in Ji’nan during the observation period.Full size imageEffects of regional differences and lockdown on air pollutantsOur results reveal that the PM10, PM2.5, NO2, SO2, CO and O3 concentrations in the urban, suburban and urban-industrial regions differed significantly between the COVID lockdown and pre-COVID periods (Figs. 3, 4).Figure 3Mean concentrations (± SD, mg/m3) of PM10, PM2.5, NO2, SO2, CO and O3 during the pre-COVID and COVID lockdown periods in 2020; the values were determined by combining the urban, suburban and urban-industrial areas at the regional scale. *, ** and *** represent significant differences between the pre-COVID and COVID lockdown periods in the same region (Duncan test, *p = 0.05; **p = 0.01; ***p = 0.001), with nonsignificant results being excluded.Full size imageFigure 4General reductions in the concentrations of major air pollutants.Full size imageNOx, one of the most important pollutants and a major health hazard, was studied in different countries across the world during COVID-19-related lockdowns. In all three regions studied herein, the highest rate of reduction in NO2 concentrations was observed during the COVID lockdown period (Fig. 4), with the NO2 levels in the COVID lockdown period being 54.02% on average lower than those during the pre-COVID period (53.07% in urban area, 48.31% in the suburban areas and 55.74% in the urban-industrial area) (Fig. 4); this reduction is greater than that reported at other sites by 26–42%11 and 14–38%18 but lower than that (50–62%) in Barcelona and Madrid in Spain33. As shown in Fig. 3E, the NO2 concentrations in the urban, suburban and urban-industrial areas were significantly higher in the pre-COVID period than in the COVID lockdown period, with the pre-COVID the NO2 levels in the urban area being 13.46% and 27.63% higher than those in the suburban and urban-industrial areas, respectively. During the COVID-19 lockdown period, the NO2 levels in urban areas were 4.69% and 31.75% higher than those in the suburban and urban-industrial areas, respectively. Blocking and controlling the air pollution associated with COVID-19 has helped reduce ground NO2 levels34 and this effect might be correlated with the tropospheric NO2 column density27. Among all sources of NO2, automobile emissions and power generation are the most important5. A systematic review confirmed that a short-term increase in the NO2 concentration in urban areas correlates to an increase in the number of pneumonia hospitalizations5,35.The trends in the CO concentration were similar to those in the NO2 level. During the COVID-lockdown period, the average CO mass concentrations in the urban, suburban and urban industrial areas were 1.08 mg/m3, 1.16 mg/m3 and 1.14 mg/m3, respectively, which decreased by 27.78%, 29.46% and 36.61%, respectively, compared with those during the pre-COVID period. The highest levels of PM10 were also observed during the pre-COVID period in the urban, suburban, and urban-industrial areas in Ji’nan (Fig. 4). The reductions in PM2.5 and CO emissions in urban and urban-industrial areas are generally higher than those in suburban areas25, supporting our findings. Notably, PM2.5 and CO are generated mainly by construction activities and from road dust, natural soil dust and dust from urban-industrial activities36. In contrast, the differences in the PM10 concentrations among the three regions were not significant during either the pre-COVID period or the COVID-lockdown period (Fig. 3A), which suggests that particles in Ji’nan are strongly diffused. However, the COVID lockdown period had a significant effect on the PM10 concentrations, with 42.86%, 44.26% and 50.60% differences in the PM10 concentration between the pre-COVID and COVID lockdown periods in the urban, suburban and urban-industrial areas, respectively (average of 44.92%, Fig. 4). The main reasons for the decreases in the concentration of PM were the severe restrictions on vehicle traffic, the cessation of industrial activities, and the stopping of construction projects, which are important sources of floating dust in the urban air37. Despite the overall consistency among the observed changes in all regions for the different air pollutants (except O3), at the regional level, some differences were statistically significant, while others were not due to the variability among stations, with the differences being more pronounced at the urban, suburban and urban-industrial stations.O3 is a secondary pollutant involved in different atmospheric reaction mechanisms and acts as both a source and sink. Generally, the impact of lockdowns on O3 was mixed, with its levels generally falling within ± 20%38, but total O3 levels remained relatively stable18. In this study, by comparing the regional mean concentrations throughout the COVID-19 period, we found that O3 concentrations were higher during the COVID lockdown period than during the pre-COVID period, especially in the urban regions (Fig. 3). Furthermore, the mean O3 concentration at all stations during the COVID lockdown period was 37.42% higher than that during the pre-COVID period (46.84% in the urban areas, 18.27% in the suburban area, and 19.84% in the urban-industrial areas) (Fig. 4); this finding is consistent with the outcomes of other studies, which reported that O3 concentrations increased by (on average) 20% during lockdowns39, potentially due, in part, to atmospheric reactivity37. The higher lockdown O3 concentrations can be attributed to the following three reasons: (1) low PM concentrations can result in more sunlight passing through the atmosphere, encouraging increased photochemical activities and thus higher O3 production40; (2) a reduction in NOx emissions increases O3 formation41; and (3) lower PM2.5 concentrations means their role as a sink for hydroperoxy radicals (HO2) is less effective, which would increase peroxy radical-mediated O3 production42. During the pre-COVID period, the O3 levels were not significantly different among the region, and the same results were observed during the COVID lockdown period. However, in the urban and urban-industrial areas, the O3 levels during the COVID lockdown period were significantly higher than those in the pre-COVID period (p  More

Gamma source and dose modelingThe general literature contains conflicting results on whether the energies of photons interacting with fungi affects the radiotrophic response. As such, we sought to control critical variables while irradiating the fungi with ionizing radiation from a sealed Cs-137 source and a UV source. The Cs-137 source emitted a photon at 662 keV along with other lower energy photons near 30 keV (Table S1).A review of previous studies was conducted to identify the gamma dose rate and total dose that should be targeted for exposure (Table S2). Those dose rates ranged from 600,000 rad delivered in 1.5 h to 0.08 rad delivered in 16 h. Even among studies examining the same fungi attributes, the total dose varied dramatically. For the present study, we used a Health Physics code to target a 50-rad dose over a one-week exposure. This dose was selected as it changes blood count observed in most humans24. We hypothesized that this dose would induce physiological changes in the fungi without causing a high rate of lethality. A MicroShield (Grove Software, Inc.) model was created to identify the quantity of radioactive material and distance between source and sample necessary to achieve the dose of 50 rad in seven days. From a sensitivity analysis of the MicroShield model, it was determined that ~ 350 µCi of Cs-137 would create a dose rate of ~ 50 rad in seven days (Fig. 1; Table S3), if placed 1.8 cm from the surface of the fungi. It should be noted that Microshield values are often conservative and likely underestimate the actual dose on target. In addition, 50 rad falls in the middle of the large range for energies previously reported in the general literature (Table S2).Figure 1Time required on target to achieve an exposure of 50 rad determined in MicroShield and based on an activity of ~ 350 µCi for Cs-137 source and the vertical distance between the source and fungus.Full size imageThe dose from the Cs-137 source on the fungal mycelium is also dependent on the radial growth of the fungus from the center plug used to initiate growth. As the fungus grows away from the source, the leading edge will experience a lower total dose of radiation. Although a uniform dose would have been ideal, a source with activity sufficient to create a uniform radiation field would have initiated a variety of safety controls deemed impractical for this experiment. The background radiation dose at the testing site in Albuquerque is approximately 10 µrem h−1; the dose at the outermost area of the Petri dish was measured at 65,553 µrad h−1. As this dose was primarily from gamma emissions, rad and rem can be considered equivalent. To validate the simulation, a dose rate study was performed using thermoluminescent dosimeters (TLD) placed at varying distances from the center of the source. The TLD placed directly under the source measured ~ 100 rad over the seven-day exposure, which is double the prediction from the simulation (50 rad; Fig. 2A). However, at a radial distance of 3.5 cm, the measured and estimated total dose over seven days were much closer, 12.3 and 11.4 rad, respectively. A comparison of the measured and estimated dose on target demonstrated a non-linear correlation (Fig. 2B), in which the simulation better approximated the dose at larger radial distances from the source.Figure 2(A) The total gamma dose on the fungal mycelial at 7 days as a function of the radial distance from the central mycelium plug based on empirical measurements (-●-) and estimated from simulations (-○-). (B) Observed correlation between the measured and estimated doses at varying radial distances.Full size imageIn order to normalize the energy deposited in the fungi from Cs-137 and the UV lamp sources, the units of MeV g−1 s−1 were selected for additional simulations. Monte Carlo N-Particle transport code (MCNP) simulations were used to determine this quantity for the Cs-137. The materials and geometry of the Petri dish and fungus used for these simulations are shown in Fig. 3. The Cs-137 was simulated as a point source located 1.5 cm from the top of the fungi. The Petri dish was set on a bakelite table. The setup was located in the center of a notional 5 m × 5 m × 5 m room with 30 cm thick concrete walls and filled with air. Leads bricks set on the table surrounded the petri dish and source. The International Commission on Radiological Protection (ICRP) material definitions did not contain data for fungal mycelia. Thus, we selected for skin as the closest approximation of the properties of the fungal mycelium25. This simulation gave a result for the energy deposited per particle as 6.53 × 10–4 MeV g−1, which for a 350 μCi activity, the rate of energy deposition was determined to be 7907 MeV g−1 s−1.Figure 3Top (upper left) and side (upper right) view of the Petri dish and fungi materials and distances used to determine energy deposition rates in MCNP. The overall geometry used for the radiation transport simulations, including the lead bricks, is shown from the top down (lower left) and from the side (lower right).Full size imageUV source and irradiationOur intent was to match the energy absorbed by the fungi to control for all variables except the photon energy difference between the Cs-137 source and UV lamp. The spectrum of energies emitted from the Cs-137 source varied significantly from those of the UV lamp, which in this case was a 30 W deuterium lamp that emitted from 185 to 400 nm (Fig. S1). This wide bandwidth represented photon energies ranging from 3.1 to 6.7 eV. The bandwidth of the UV exposure was limited to 300–350 nm using a 50-nm bandpass filter centered at 325 nm to ensure that incident photons would be in the UV energy range and not form ozone. Because we chose to match the overall energy deposited from the UV source to the gamma source it was necessary to attenuate the beam to the right power level. We assumed that all the UV energy would be absorbed near the surface rather than in the bulk since the fungi were melanized. This simplified the calculations and reduced risk, given the challenge of accurately estimating the absorbance of the fungi. The power deposited by the gamma source was calculated as the rate of energy deposition was determined to be 7907 MeV g−1 s−1 (1.3 nW g−1 s−1). Given the initial size of the plug was 1 cm in diameter, the desired lamp fluence needed to be ~ 2.8 nW cm−2. Across the spectrum of interest, the lamp power was determined to be 3.202 × 10–4 mW, thus requiring an attenuation of 8.7 × 10–9 (OD 8.06), reducing the lamp power to ~ 3 pW cm−2 and achieving a reasonably close power density to the target. Due to the sensitivity of UV detectors, the required power densities could not be measured directly. Alternatively, we measured the neutral density filters to verify the prescription was indeed correct.Response of P. variotii to irradiationUniform plugs (~ 5-mm in diameter) of actively growing mycelia of P. variotii were cut using the end of a Pasteur pipette and transferred a Petri plate containing potato dextrose agar (PDA) one day prior to initiating exposure experiments. The diameter of the mycelium was measured from four images, separated by precisely six hours, over the course of seven days and used to measure the growth rate. Differences in the pigmentation of the fungi under the different conditions was quantified in Fiji26 through analysis of grayscale images collected at day seven, following the method described by Brilhante et al.27 A ratiometric value was derived from the grayscale values and the white background, which corrected for variations in lighting across or between images.Significant differences in the pigmentation but not growth rates of P. variotii were associated with exposure to UV and gamma to irradiation, based on One-Way ANOVA analyses (Fig. 4A; Table S4). P. variotii is a ubiquitous filamentous fungus commonly inhabiting soil, decaying plants, and food products and was reported to be present on the surface of the walls of Unit-4 at ChNPP22,28. P. variotii is also a common food contaminant and is resistant to high temperature and metals29,30, despite being more sensitive to gamma irradiation than other fungi such as Aspergillus fumigatus31. In the present study, we hypothesized that positive radiation-induced effects in P. variotii would result in enhanced growth rates due to gamma irradiation. Across all conditions, the average growth rate of P. variotii was ~ 5.6 ± 0.9 mm d−1 (mean ± standard deviation). While the growth rate of P. variotii exposed to gamma irradiation was greater compared with the control and UV-irradiated samples (Fig. 4A), the difference in the mean growth rates was not significant (P = 0.255) by ANOVA.Figure 4(A) Growth rate and pigmentation of control (orange square), gamma- (blue square), and UV- (red square) irradiated cultures of P. variotti (mean ± standard deviation). (B) Estimated total irradiation dose experienced by the mycelial as a function of the distance from the central source. Exponential decay fit: − 3.6 + 105.7*exp(− 0.75*x); Adjusted R2 = 0.998. (C) Graphical representation of the irradiation dose based on the growth rate and duration of exposure for zones of mycelia as a function of radial distance from the central plug.Full size imageWe also hypothesized that the pigmentation of P. variotii would increase with exposure to gamma and UV irradiation. While P. variotti does not produce melanin, it does produce a pigment, Ywa1, from a polyketide synthesis (PKS) gene cluster and has been shown to protect the fungus against UV-C irradiation28. In some melanized fungi, Ywa1 serves as precursor and can be hydrolyzed to 1,3,6,8-tetrahydroxynaphthalen (T4HN). T4HN may then be converted to 1,8-dihydroxynaphthalene (1,8-DHN) melanin through the DHN pathway32. However, Lim et al.28 concluded that P. variotii does not produce true melanin as the pigmentation was maintained when the DHN-melanin pathway was inhibited. Significant differences in the pigmentation of P. variotii were observed among the three different sample types (P  More

Banks, N. C., Paini, D. R., Bayliss, K. L. & Hodda, M. The role of global trade and transport network topology in the human-mediated dispersal of alien species. Ecol. Lett. 18, 188–199 (2015).
Google Scholar 
Epanchin-Niell, R. et al. Controlling invasive species in complex social landscapes. Front. Ecol. Environ. 8, 210–216 (2009).
Google Scholar 
Charles, H. & Dukes, J. S. in Biological Invasions (ed. Nentwig, W.) 217–237 (Springer, 2007). https://doi.org/10.1007/978-3-540-36920-2_13Gallardo, B., Clavero, M., Sánchez, M. & Vilà, M. Global ecological impacts of invasive species in aquatic ecosystems. Glob. Change Biol. 22, 151–163 (2016).
Google Scholar 
Diagne, C. et al. High and rising economic costs of biological invasions worldwide. Nature 592, 571–576 (2021).CAS 

Google Scholar 
Sardain, A., Sardain, E. & Leung, B. Global forecasts of shipping traffic and biological invasions to 2050. Nat. Sustain. 2, 274–282 (2019).
Google Scholar 
Epanchin-Niell, R. S. & Hastings, A. Controlling established invaders: integrating economics and spread dynamics to determine optimal management. Ecol. Lett. 13, 528–541 (2010).
Google Scholar 
Chades, I. et al. General rules for managing and surveying networks of pests, diseases, and endangered species. Proc. Natl. Acad. Sci. USA 108, 8323–8328 (2011).CAS 

Google Scholar 
Epanchin-Niell, R. S. & Wilen, J. E. Optimal spatial control of biological invasions. J. Environ. Econ. Manag. 63, 260–270 (2012).
Google Scholar 
Epanchin-Niell, R. S. & Wilen, J. E. Individual and cooperative management of invasive species in human-mediated landscapes. Am. J. Agric. Econ. 97, 180–198 (2015).
Google Scholar 
Aadland, D., Sims, C. & Finnoff, D. Spatial dynamics of optimal management in bioeconomic systems. Comput. Econ. 45, 545–577 (2015).
Google Scholar 
Baker, C. M. Target the source: optimal spatiotemporal resource allocation for invasive species control. Conserv. Lett. 10, 41–48 (2017).
Google Scholar 
Bushaj, S., Büyüktahtakın, İ. E., Yemshanov, D. & Haight, R. G. Optimizing surveillance and management of emerald ash borer in urban environments. Nat. Res. Model. 34, e12267 (2021).
Google Scholar 
Fischer, S. M., Beck, M., Herborg, L.-M. & Lewis, M. A. Managing aquatic invasions: optimal locations and operating times for watercraft inspection stations. J. Environ. Manag. 283, 111923 (2021).
Google Scholar 
Büyüktahtakın, İ. E. & Haight, R. G. A review of operations research models in invasive species management: state of the art, challenges, and future directions. Ann. Oper. Res. 271, 357–403 (2018).
Google Scholar 
Epanchin-Niell, R. S. Economics of invasive species policy and management. Biol. Invasions 19, 3333–3354 (2017).
Google Scholar 
Bodin, Ö. et al. Improving network approaches to the study of complex social–ecological interdependencies. Nat. Sustain. 2, 551–559 (2019).CAS 

Google Scholar 
Nowzari, C., Precaido, V. M. & Pappas, G. J. Analysis and control of epidemics: a survey of spreading processes on complex networks. IEEE Control Syst. 36, 26–46 (2016).
Google Scholar 
Newman, M. E. J. Spread of epidemic disease on networks. Phys. Rev. E 66, 016128 (2002).CAS 

Google Scholar 
Kempe, D., Kleinberg, J. & Tardos, E. Maximizing the spread of influence through a social network. In Proc. 9th ACM SIGKDD International Conference on Knowledge Discovery and Data Mining 137–146 (ACM Press, 2003).Pastor-Satorras, R. & Vespignani, A. Immunization of complex networks. Phys. Rev. E 65, 036104 (2002).
Google Scholar 
Pastor-Satorras, R., Castellano, C., Van Mieghem, P. & Vespignani, A. Epidemic processes in complex networks. Rev. Mod. Phys. 87, 925–979 (2015).
Google Scholar 
Holme, P., Kim, B. J., Yoon, C. N. & Han, S. K. Attack vulnerability of complex networks. Phys. Rev. E 65, 056109 (2002).
Google Scholar 
Muirhead, J. R. & Macisaac, H. J. Development of inland lakes as hubs in an invasion network. J. Appl. Ecol. 42, 80–90 (2005).
Google Scholar 
de la Fuente, B., Saura, S. & Beck, P. S. Predicting the spread of an invasive tree pest: the pine wood nematode in southern europe. J. Appl. Ecol. 55, 2374–2385 (2018).
Google Scholar 
Minor, E. S. & Urban, D. L. A graph-theory framework for evaluating landscape connectivity and conservation planning. Conserv. Biol. 22, 297–307 (2008).
Google Scholar 
Morel-Journel, T., Assa, C. R., Mailleret, L. & Vercken, E. Its all about connections: hubs and invasion in habitat networks. Ecol. Lett. 22, 313–321 (2019).
Google Scholar 
Perry, G. L. W., Moloney, K. A. & Etherington, T. R. Using network connectivity to prioritise sites for the control of invasive species. J. Appl. Ecol. 54, 1238–1250 (2017).
Google Scholar 
Kvistad, J. T., Chadderton, W. L. & Bossenbroek, J. M. Network centrality as a potential method for prioritizing ports for aquatic invasive species surveillance and response in the Laurentian Great Lakes. Manag. Biol. Invasions 10, 403 (2019).
Google Scholar 
Haight, R. G., Kinsley, A. C., Kao, S.-Y., Yemshanov, D. & Phelps, N. B. Optimizing the location of watercraft inspection stations to slow the spread of aquatic invasive species. Biol. Invasions 23, 3907–3919 (2021).
Google Scholar 
McEachran, M. C. et al. Stable isotopes indicate that zebra mussels (Dreissena polymorpha) increase dependence of lake food webs on littoral energy sources. Freshw, Biol. 64, 183–196 (2019).CAS 

Google Scholar 
Karatayev, A. Y., Burlakova, L. E. & Padilla, D. K. in Invasive Aquatic Species of Europe. Distribution, Impacts and Management (eds Leppäkoski, E. et al.) 433–446 (Springer, 2002).Prescott, T. H., Claudi, R. & Prescott, K. L. Impact of Dreissenid mussels on the infrastructure of dams and hydroelectric power plants. In Quagga and Zebra Mussels (eds Nalepa, T. F. & Schloesser, D. W.) 243–258 (CRC Press, 2013).Invasive Species of Aquatic Plants and Wild Animals in Minnesota: Annual Report for 2020 (Minnesota Department of Natural Resources, 2020).Kanankege, K. S., Alkhamis, M. A., Phelps, N. B. & Perez, A. M. A probability co-kriging model to account for reporting bias and recognize areas at high risk for zebra mussels and eurasian watermilfoil invasions in Minnesota. Front. Vet. Sci. 4, 231 (2018).
Google Scholar 
Mallez, S. & McCartney, M. Dispersal mechanisms for zebra mussels: population genetics supports clustered invasions over spread from hub lakes in Minnesota. Biol. Invasions 20, 2461–2484 (2018).
Google Scholar 
Kao, S.-Y. Z. et al. Network connectivity of Minnesota waterbodies and implications for aquatic invasive species prevention. Biol. Invasions 23, 3231–3242 (2021).
Google Scholar 
Kleinberg, J. M. Authoritative sources in a hyperlinked environment. In Proc. 9th Annual ACM-SIAM Symposium on Discrete Algorithms 668–677 (1998).McDonald-Madden, E. et al. Using food-web theory to conserve ecosystems. Nat. Commun. 7, 10245 (2016).CAS 

Google Scholar 
Bossenbroek, J. M., Kraft, C. E. & Nekola, J. C. Prediction of long-distance dispersal using gravity models: zebra mussel invasion of inland lakes. Ecol. Appl. 11, 1778–1788 (2001).
Google Scholar 
Leung, B., Bossenbroek, J. M. & Lodge, D. M. Boats, pathways, and aquatic biological invasions: estimating dispersal potential with gravity models. Biol. Invasions 8, 241–254 (2006).
Google Scholar 
Beger, M. et al. Integrating regional conservation priorities for multiple objectives into national policy. Nat. Commun. 6, 8208 (2015).Runting, R. K. et al. Larger gains from improved management over sparing–sharing for tropical forests. Nat. Sustain. 2, 53–61 (2019).
Google Scholar 
Kinsley, A. C. et al. AIS Explorer: prioritization for watercraft inspections. A decision-support tool for aquatic invasive species management. J. Environ. Manage. 314, 115037 (2022).
Google Scholar 
Vander Zanden, M. J. & Olden, J. D. A management framework for preventing the secondary spread of aquatic invasive species. Can. J. Fish. Aquat. Sci. 65, 1512–1522 (2008).
Google Scholar 
Kanankege, K. S. et al. Lessons learned from the stakeholder engagement in research: application of spatial analytical tools in one health problems. Front. Vet. Sci. 7, 254 (2020).
Google Scholar 
Kroetz, K. & Sanchirico, J. The bioeconomics of spatial-dynamic systems in natural resource management. Annu. Rev. Resour. Econ. 7, 189–207 (2015).
Google Scholar 
Cade, B. S. & Noon, B. R. A gentle introduction to quantile regression for ecologists. Front. Ecol. Environ. 1, 412–420 (2003).
Google Scholar 
Koenker, R. in Asymptotic Statistics (eds Mandl, P. & Hušková, M.) 349–359 (Springer, 1994).Ashander, J. Analysis code and data for ‘Guiding large-scale management of invasive species using network metrics’. figshare https://doi.org/10.6084/m9.figshare.14402447 (2021). More

Climate and land cover dataOur study of tropical dry season dynamics required climatic variables with high temporal resolution (i.e., daily) and full coverage of tropic regions. To reduce uncertainties associated with the choice of precipitation (P) and evapotranspiration (Ep or E) datasets, we used an ensemble of eight precipitation products, three reanalysis-based products for Ep, and one satellite-based land E product. These precipitation datasets were derived four gauge-based or satellite observation (CHIRPS58, GPCC59, CPC-U60 and PERSIANN-CDR61), three reanalyses (ERA-562, MERRA-263, and PGF64) and a multi-source weighted ensemble product (MSWEP v2.865). The potential evapotranspiration (Ep) was calculated using the FAO Penman–Monteith equation66 (Eqs. (1, 2)), which requires meteorological inputs of wind speed, net radiation, air temperature, specific humidity, and surface pressure. We derived these meteorological variables from the three reanalysis products (ERA-5, MERRA-2, and GLDAS-2.067). Since PGF reanalysis lacked upward short- and long-wave radiation output and thus net radiation, we used available meteorological outputs from GLDAS-2.0 instead, which was forced entirely with the PGF input data.$${Ep}=frac{0.408cdot triangle cdot left({R}_{n}-Gright)+gamma cdot frac{900}{T+273}cdot {u}_{2}cdot left({e}_{s}-{e}_{a}right)}{triangle +{{{{{rm{gamma }}}}}}cdot left(1+0.34cdot {u}_{2}right)}$$
(1)
$${VPD}={e}_{s}-{e}_{a}=0.6108cdot {e}^{frac{17.27cdot T}{T+237.3}}cdot left(1-frac{{RH}}{100}right)$$
(2)
Where Ep is the potential evapotranspiration (mm day−1). Rn is net radiation at the surface (MJ m−2 day−1), T is mean daily air temperature at 2 m height (°C), ({u}_{2}) is wind speed at 2 m height (m s−1), ((,{e}_{s}-{e}_{a})) is the vapor pressure deficit of the air (kPa), ({RH}) is the relative air humidity near surface (%), ∆ is the slope of the saturation vapor pressure-temperature relationship (kPa °C−1), γ is the psychrometric constant (kPa °C−1), G is the soil heat flux (MJ m−2 day−1, is often ignored for daily time steps G ≈ 0).We derived the daily evapotranspiration data from the Global Land Evaporation Amsterdam Model (GLEAM v3.3a68), which is a set of algorithms dedicated to developing terrestrial evaporation and root-zone soil moisture data. GLEAM fully assimilated the satellite-based soil moisture estimates from ESA CCI, microwave L-band vegetation optical depth (VOD), reanalysis-based temperature and radiation, and multi-source precipitation forcings. The direct assimilation of observed soil moisture allowed us to detect true soil moisture dynamic and its impacts on evapotranspiration. Besides, the incorporation of VOD, which is closely linked to vegetation water content69,70, allowed us to detect the effect of water stress, heat stress, and vegetation phenological constraints on evaporation. Other observation-driven ET products from remote-sensing physical estimation and flux-tower are not included due to their low temporal resolution (i.e., monthly)71 or short duration72,73. ET outputs of reanalysis products are not considered in our analysis, because the assimilation systems lack explicit representation of inter-annual variability of vegetation activities and thus may not fully capture hydrological response to vegetation changes62,63,67.We used land cover maps for the year 2001 from the Moderate-Resolution Imaging Spectroradiometer (MODIS, MCD12C1 C574) based on the IGBP classification scheme to exclude water-dominated and sparely-vegetated pixels (like Sahara, Arabian Peninsula). All climate and land cover datasets mentioned above were remapped to a common 0.25° × 0.25° grid and unified to daily resolution. The main characteristics of the datasets mentioned above are summarized in Supplementary Table 1.Outputs of CMIP6 simulationsTo understand how modeled dry season changes compare with observed changes, we analyzed outputs from the “historical” (1983-2014) runs of 34 coupled models participating in the 6th Coupled Model Inter-comparison Project75 (CMIP6, Supplementary Table 3). We used these models because they offered daily outputs of all climatic variables needed for our analysis, including precipitation, latent heat (convert to E), and multiple meteorological variables for Ep (air temperature, surface specific humidity, wind speed, and net radiation). All outputs were remapped to a common 1.0° × 1.0° grid and unified to daily resolution.Defining dry season length and timingFor each grid cell and each dry season definition (P  More

The management of introduced species, whether kudzu or zebra mussels, is costly and complex. Now, a paper reports a workable, effective solution that harnesses network analyses of ecological phenomena.Invasive species can pose severe economic and environmental problems, costing more than US$1 trillion worldwide since 1970 (ref. 1). Yet managing this human-driven issue is difficult in itself. The regions involved can be vast — entire continents or countries, for instance — while budgets are typically limited. As well, the sites potentially affected and management options can be numerous. Real systems (for example, all the lakes in the United States) can have thousands of locations that could potentially be infested. By contrast, considering just 40 locations means dealing theoretically with over 1 trillion unique combinations (240) where management could be applied (for instance, to reduce the number of invasive species leaving infested areas or entering uninfested ones). Given these constraints, a key problem is how and where to deploy control measures such as invasive-species removal. While sophisticated optimization approaches exist2, which use mathematical rules to exclude most suboptimal combinations and quickly zoom in to which locations should be managed to minimize new invasions, these algorithms are generally unfeasible for very large systems. Now, writing in Nature Sustainability, Ashander et al.3 demonstrate that simpler network metrics revealing linkages between patches can provide solutions that are often comparable to the more complex optimization algorithms. More

Thank you for visiting nature.com. You are using a browser version with limited support for CSS. To obtain
the best experience, we recommend you use a more up to date browser (or turn off compatibility mode in
Internet Explorer). In the meantime, to ensure continued support, we are displaying the site without styles
and JavaScript. More

Gause, G. F. Experimental analysis of Vito Volterra’s mathematical theory of the struggle for existence. Science 79, 16–17 (1934).ADS 
CAS 
PubMed 
Article 

Google Scholar 
Amarasekare, P. Competitive coexistence in spatially structured environments: A synthesis. Ecol. Lett. 6, 1109–1122 (2003).Article 

Google Scholar 
HilleRisLambers, J., Adler, P. B., Harpole, W. S., Levine, J. M. & Mayfield, M. M. Rethinking community assembly through the lens of coexistence theory. Annu. Rev. Ecol. Evol. Syst. 43, 227–248 (2012).Article 

Google Scholar 
Wisz, M. S. et al. The role of biotic interactions in shaping distributions and realised assemblages of species: Implications for species distribution modelling. Biol. Rev. 88, 15–30 (2013).PubMed 
Article 

Google Scholar 
Frey, S., Fisher, J. T., Burton, A. C. & Volpe, J. P. Investigating animal activity patterns and temporal niche partitioning using camera-trap data: Challenges and opportunities. Remote Sens. Ecol. Conserv. 3, 123–132 (2017).Article 

Google Scholar 
Davis, C. L. et al. Ecological correlates of the spatial co-occurrence of sympatric mammalian carnivores worldwide. Ecol. Lett. 21, 1401–1412 (2018).PubMed 
Article 

Google Scholar 
Durant, S. M. Competition refuges and coexistence: An example from Serengeti carnivores. J. Anim. Ecol. 67, 370–386 (1998).Article 

Google Scholar 
Fedriani, J. M., Fuller, T. K., Sauvajot, R. M. & York, E. C. Competition and intraguild predation among three sympatric carnivores. Oecologia 125, 258–270 (2000).ADS 
PubMed 
Article 

Google Scholar 
Kamler, J. F., Ballard, W. B., Gilliland, R. L. & Mote, K. Spatial relationships between swift foxes and coyotes in northwestern Texas. Can. J. Zool. 81, 168–172 (2003).Article 

Google Scholar 
Vanak, A. T. et al. Moving to stay in place: Behavioral mechanisms for coexistence of African large carnivores. Ecology 94, 2619–2631 (2013).PubMed 
Article 

Google Scholar 
Donadio, E. & Buskirk, S. W. Diet, morphology, and interspecific killing in carnivora. Am. Nat. 167, 524–536 (2006).PubMed 
Article 

Google Scholar 
Tsunoda, H. et al. Food niche segregation between sympatric golden jackals and red foxes in central Bulgaria. J. Zool. 303, 64–71 (2017).Article 

Google Scholar 
Palomares, F. & Caro, T. M. Interspecific killing among mammalian carnivores. Am. Nat. 153, 492–508 (1999).CAS 
PubMed 
Article 

Google Scholar 
Linnell, J. D. C. & Strand, O. Interference interactions, co-existence and conservation of mammalian carnivores. Divers. Distrib. 6, 169–176 (2000).Article 

Google Scholar 
Kamler, J. F., Stenkewitz, U., Klare, U., Jacobsen, N. F. & MacDonald, D. W. Resource partitioning among cape foxes, bat-eared foxes, and black-backed jackals in South Africa. J. Wildl. Manag. 76, 1241–1253 (2012).Article 

Google Scholar 
Di Bitetti, M. S., Di Blanco, Y. E., Pereira, J. A., Paviolo, A. & Pírez, I. J. Time Partitioning favors the coexistence of sympatric crab-eating foxes (Cerdocyon thous) and Pampas Foxes (Lycalopex gymnocercus). J. Mammal. 90, 479–490 (2009).Article 

Google Scholar 
Lesmeister, D. B., Nielsen, C. K., Schauber, E. M. & Hellgren, E. C. Spatial and temporal structure of a mesocarnivore guild in Midwestern North America. Wildl. Monogr. 191, 1–61 (2015).Article 

Google Scholar 
Di Bitetti, M. S., De Angelo, C. D., Di Blanco, Y. E. & Paviolo, A. Niche partitioning and species coexistence in a Neotropical felid assemblage. Acta Oecologica 36, 403–412 (2010).ADS 
Article 

Google Scholar 
Monterroso, P., Alves, P. C. & Ferreras, P. Plasticity in circadian activity patterns of mesocarnivores in southwestern Europe: Implications for species coexistence. Behav. Ecol. Sociobiol. 68, 1403–1417 (2014).Article 

Google Scholar 
Tsunoda, H., Ito, K., Peeva, S., Raichev, E. & Kaneko, Y. Spatial and temporal separation between the golden jackal and three sympatric carnivores in a human-modified landscape in central Bulgaria. Zool. Ecol. 28, 172–179 (2018).Article 

Google Scholar 
Tsunoda, H. et al. Spatio-temporal partitioning facilitates mesocarnivore sympatry in the Stara Planina Mountains, Bulgaria. Zoology 141, 125801 (2020).PubMed 
Article 

Google Scholar 
Ramesh, T., Kalle, R., Sankar, K. & Qureshi, Q. Spatio-temporal partitioning among large carnivores in relation to major prey species in Western Ghats. J. Zool. 287, 269–275 (2012).Article 

Google Scholar 
Gómez-Ortiz, Y., Monroy-Vilchis, O. & Castro-Arellano, I. Temporal coexistence in a carnivore assemblage from central Mexico: Temporal-domain dependence. Mammal Res. 64, 333–342 (2019).Article 

Google Scholar 
Ridout, M. S. & Linkie, M. Estimating overlap of daily activity patterns from camera trap data. J. Agric. Biol. Environ. Stat. 14, 322–337 (2009).MathSciNet 
MATH 
Article 

Google Scholar 
Meredith, M. & Ridout, M. Overlap: Estimates of coefficient of overlapping for animal activity patterns. https://cran.r-project.org/web/packages/overlaphttps://cran.r-project.org/web/packages/overlap/index.html (2018).Marinho, P. H., Fonseca, C. R., Sarmento, P., Fonseca, C. & Venticinque, E. M. Temporal niche overlap among mesocarnivores in a Caatinga dry forest. Eur. J. Wildl. Res. 66, 1–13 (2020).Article 

Google Scholar 
Vilella, M., Ferrandiz-Rovira, M. & Sayol, F. Coexistence of predators in time: Effects of season and prey availability on species activity within a Mediterranean carnivore guild. Ecol. Evol. 10, 11408–11422 (2020).PubMed 
PubMed Central 
Article 

Google Scholar 
Zhao, G. et al. Spatio-temporal coexistence of sympatric mesocarnivores with a single apex carnivore in a fine-scale landscape. Glob. Ecol. Conserv. 21, e00897 (2020).Article 

Google Scholar 
Farmer, M. J., Allen, M. L., Olson, E. R., Van Stappen, J. & Van Deelen, T. R. Agonistic interactions and island biogeography as drivers of carnivore spatial and temporal activity at multiple scales. Can. J. Zool. 99, 309–317 (2021).Article 

Google Scholar 
Watabe, R. & Saito, M. U. Diel activity patterns of three sympatric medium-sized carnivores during winter and spring in a heavy snowfall area in northeastern Japan. Mammal Study 46, 69–75 (2021).Article 

Google Scholar 
Lashley, M. A. et al. Estimating wildlife activity curves: comparison of methods and sample size. Sci. Rep. 8, 4173 (2018).ADS 
PubMed 
PubMed Central 
Article 
CAS 

Google Scholar 
Niedballa, J., Wilting, A., Sollmann, R., Hofer, H. & Courtiol, A. Assessing analytical methods for detecting spatiotemporal interactions between species from camera trapping data. Remote Sens. Ecol. Conserv. 5, 272–285 (2019).Article 

Google Scholar 
Karanth, K. U. et al. Spatio-temporal interactions facilitate large carnivore sympatry across a resource gradient. Proc. R. Soc. B Biol. Sci. 284, 20161860 (2017).Article 

Google Scholar 
Cusack, J. J. et al. Revealing kleptoparasitic and predatory tendencies in an African mammal community using camera traps: A comparison of spatiotemporal approaches. Oikos 126, 812–822 (2017).Article 

Google Scholar 
Balme, G. et al. Big cats at large: density, structure, and spatio-temporal patterns of a leopard population free of anthropogenic mortality. Popul. Ecol. 61, 256–267 (2019).Article 

Google Scholar 
Li, Z. et al. Coexistence of two sympatric flagship carnivores in the human-dominated forest landscapes of Northeast Asia. Landsc. Ecol. 34, 291–305 (2019).Article 

Google Scholar 
Lahkar, D., Ahmed, M. F., Begum, R. H., Das, S. K. & Harihar, A. Inferring patterns of sympatry among large carnivores in Manas National Park: A prey-rich habitat influenced by anthropogenic disturbances. Anim. Conserv. 24, 589–601 (2021).Article 

Google Scholar 
Paúl, M. J., Layna, J. F., Monterroso, P. & Álvares, F. Resource partitioning of sympatric African Wolves (Canis lupaster) and side-striped jackals (Canis adustus) in an arid environment from West Africa. Diversity 12, 477 (2020).Article 

Google Scholar 
Prat-Guitart, M., Onorato, D. P., Hines, J. E. & Oli, M. K. Spatiotemporal pattern of interactions between an apex predator and sympatric species. J. Mammal. 101, 1279–1288 (2020).Article 

Google Scholar 
Stone, L. & Roberts, A. The checkerboard score and species distributions. Oecologia 85, 74–79 (1990).ADS 
PubMed 
Article 

Google Scholar 
Griffith, D. M., Veech, J. A. & Marsh, C. J. Cooccur: Probabilistic species co-occurrence analysis in r. J. Stat. Softw. 69, 1–17 (2016).Article 

Google Scholar 
Noor, A., Mir, Z. R., Veeraswami, G. G. & Habib, B. Activity patterns and spatial co-occurrence of sympatric mammals in the moist temperate forest of the Kashmir Himalaya, India. Folia Zool. 66, 231–241 (2017).Article 

Google Scholar 
de Satgé, J., Teichman, K. & Cristescu, B. Competition and coexistence in a small carnivore guild. Oecologia 184, 873–884 (2017).ADS 
PubMed 
Article 

Google Scholar 
Kass, J. M., Tingley, M. W., Tetsuya, T. & Koike, F. Co-occurrence of invasive and native carnivorans affects occupancy patterns across environmental gradients. Biol. Invasions 22, 2251–2266 (2020).Article 

Google Scholar 
Louppe, V., Herrel, A., Pisanu, B., Grouard, S. & Veron, G. Assessing occupancy and activity of two invasive carnivores in two Caribbean islands: implications for insular ecosystems. J. Zool. 313, 182–194 (2020).Article 

Google Scholar 
Proulx, G. et al. World distribution and status of the genus Martes in 20. In Martens and Fishers (Martes) in Human-Altered Environments (eds Harrison, D. J. et al.) 21–76 (Springer, Berlin, 2005). https://doi.org/10.1007/b99487.Chapter 

Google Scholar 
Ohdachi, S. D., Ishibashi, Y., Iwasa, M., Fukuki, D. & Saitoh, T. The Wild Mammals of Japan 2nd edn. (Shokadoh Book Seller, Kyoto, 2015).
Google Scholar 
Kauhala, K. & Saeki, M. Nyctereutes procyonoides. The IUCN Red List of Threatened Species. https://www.iucnredlist.org/species/14925/85658776 (2016).Yamamoto, Y. Comparative analyses on food habits of Japanese marten, red fox, badger and raccoon dog in the Mt. Nyugasa, Nagano Prefecture, Japan. Nat. Environ. Sci. Res. 7, 45–52 (1994) (in Japanese with English summary).
Google Scholar 
Hisano, M. et al. A comparison of visual and genetic techniques for identifying Japanese marten scats enabling diet examination in relation to seasonal food availability in a sub-alpine area of Japan. Zool. Sci. 34, 137–146 (2017).Article 

Google Scholar 
Lindstrom, E. R., Brainerd, S. M., Helldin, J. O. & Overskaug, K. Pine marten-red fox interactions: A case of intraguild predation?. Ann. Zool. Fenn. 32, 123–130 (1995).
Google Scholar 
Waggershauser, C. N., Ruffino, L., Kortland, K. & Lambin, X. Lethal interactions among forest-grouse predators are numerous, motivated by hunger and carcasses, and their impacts determined by the demographic value of the victims. Ecol. Evol. 11, 7164–7186 (2021).PubMed 
PubMed Central 
Article 

Google Scholar 
Watabe, R., Saito, M. U., Enari, H. S. & Enari, H. Mammalian fauna of the Kaminagawa Experimental Forest of Yamagata University detected by camera traps. Tohoku J. For. Sci. 25, 37–40 (2020) (in Japanese).
Google Scholar 
Hofmeester, T. R., Rowcliffe, J. M. & Jansen, P. A. A simple method for estimating the effective detection distance of camera traps. Remote Sens. Ecol. Conserv. 3, 81–89 (2017).Article 

Google Scholar 
Di Bitetti, M. S., Paviolo, A. & De Angelo, C. Camera trap photographic rates on roads vs. off roads: Location does matter. Mastozoología Neotrop. 21, 37–46 (2014).
Google Scholar 
Borcard, D. & Legendre, P. Is the Mantel correlogram powerful enough to be useful in ecological analysis? A simulation study. Ecology 93, 1473–1481 (2012).PubMed 
Article 

Google Scholar 
Oksanen, J. et al. Vegan: community ecology package. https://cran.r-project.org/web/packages/veganhttps://cran.r-project.org/web/packages/vegan/index.html (2019).R Core Team. R: a language environment for statistical computing. r foundation for statistical computing, Vienna, Austria. https://www.r-project.org/https://www.r-project.org/ (2021).Linkie, M. & Ridout, M. S. Assessing tiger-prey interactions in Sumatran rainforests. J. Zool. 284, 224–229 (2011).Article 

Google Scholar 
Watabe, R. & Saito, M. U. Effects of vehicle-passing frequency on forest roads on the activity patterns of carnivores. Landsc. Ecol. Eng. 17, 225–231 (2021).Article 

Google Scholar 
Furukawa, G. genkiFurukawa/rSetDayNightAttr documentation. https://rdrr.io/github/genkiFurukawa/rSetDayNightAhttps://rdrr.io/github/genkiFurukawa/rSetDayNightAttr/ (2019).Mielke, P. W., Berry, K. J. & Johnson, E. S. Multi-response permutation procedures for a priori classifications. Commun. Stat. Theory Methods 5, 1409–1424 (1976).MATH 
Article 

Google Scholar 
Kronfeld-Schor, N. & Dayan, T. Partitioning of time as an ecological resource. Annu. Rev. Ecol. Evol. Syst. 34, 153–181 (2003).Article 

Google Scholar 
Monterroso, P., Alves, P. C. & Ferreras, P. Catch me if you can: Diel activity patterns of mammalian prey and predators. Ethology 119, 1044–1056 (2013).Article 

Google Scholar 
Hendrichsen, D. K. & Tyler, N. J. C. How the timing of weather events influences early development in a large mammal. Ecology 95, 1737–1745 (2014).CAS 
PubMed 
Article 

Google Scholar 
Herfindal, I. et al. Weather affects temporal niche partitioning between moose and livestock. Wildlife Biol. https://doi.org/10.2981/wlb.00275 (2017).Article 

Google Scholar 
Haswell, P. M., Jones, K. A., Kusak, J. & Hayward, M. W. Fear, foraging and olfaction: How mesopredators avoid costly interactions with apex predators. Oecologia 187, 573–583 (2018).ADS 
PubMed 
PubMed Central 
Article 

Google Scholar 
Barrull, J. et al. Factors and mechanisms that explain coexistence in a Mediterranean carnivore assemblage: An integrated study based on camera trapping and diet. Mamm. Biol. 79, 123–131 (2014).Article 

Google Scholar 
Tattersall, E. R., Burgar, J. M., Fisher, J. T. & Burton, A. C. Boreal predator co-occurrences reveal shared use of seismic lines in a working landscape. Ecol. Evol. 10, 1678–1691 (2020).PubMed 
PubMed Central 
Article 

Google Scholar 
Moll, R. J. et al. Humans and urban development mediate the sympatry of competing carnivores. Urban Ecosyst. 21, 765–778 (2018).Article 

Google Scholar 
McCreadie, J. W. & Bedwell, C. R. Patterns of co-occurrence of stream insects and an examination of a causal mechanism: Ecological checkerboard or habitat checkerboard?. Insect Conserv. Divers. 6, 105–113 (2013).Article 

Google Scholar  More