Philippa Kaur

Foster, G. & Rahmstorf, S. Global temperature evolution 1979–2010. Environ. Res. Lett. 6, 044022 (2011).ADS 
Article 

Google Scholar 
Cahill, N., Rahmstorf, S. & Parnell, A. C. Change points of global temperature. Environ. Res. Lett. 10, 084002 (2015).ADS 
Article 

Google Scholar 
Yan, X. H. et al. The global warming hiatus: Slowdown or redistribution?. Earth’s Future 4, 472–482 (2016).ADS 
CAS 
PubMed 
PubMed Central 
Article 

Google Scholar 
Karl, T. R. et al. Possible artifacts of data biases in the recent global surface warming hiatus. Science 348, 5632 (2015).Article 

Google Scholar 
Cohen, J. L., Furtado, J. C., Barlow, M., Alexeev, V. A. & Cherry, J. E. Asymmetric seasonal temperature trends. Geophys. Res. Lett. 39, 04705. https://doi.org/10.1029/2011GL050582 (2012).ADS 
Article 

Google Scholar 
Pepin, N. C. & Lundquist, J. D. Temperature trends at high elevations: patterns across the globe. Geophys. Res. Lett. 35, 14 (2008).Article 

Google Scholar 
Rangwala, I. & Miller, J. R. Climate change in mountains: a review of elevation-dependent warming and its possible causes. Clim. Change 114, 527–547 (2012).ADS 
Article 

Google Scholar 
Wang, Q., Fan, X. & Wang, M. Recent warming amplification over high elevation regions across the globe. Clim. Dyn. 43, 87–101 (2014).CAS 
Article 

Google Scholar 
Fan, X., Wang, Q., Wang, M. & Jiménez, C. V. Warming amplification of minimum and maximum temperatures over high-elevation regions across the globe. PLoS ONE 10, e0140213 (2015).PubMed 
PubMed Central 
Article 

Google Scholar 
Pepin, N. et al. Elevation-dependent warming in mountain regions of the world. Nat. Clim. Chang. 5, 424 (2015).ADS 
Article 

Google Scholar 
Piccarreta, M., Lazzari, M. & Pasini, A. Trends in daily temperature extremes over the Basilicata region (southern Italy) from 1951 to 2010 in a Mediterranean climatic context. Int. J. Climatol. 35, 1964–1975 (2015).Article 

Google Scholar 
Gonzalez-Hidalgo, J. C., Peña-Angulo, D., Brunetti, M. & Cortesi, N. Recent trend in temperature evolution in Spanish mainland (1951–2010): from warming to hiatus. Int. J. Climatol. 36, 2405–2416 (2016).Article 

Google Scholar 
McCullough, I. M. et al. High and dry: high elevations disproportionately exposed to regional climate change in Mediterranean-climate landscapes. Landsc. Ecol. 31, 1063–1075 (2016).Article 

Google Scholar 
Sanz-Elorza, M., Dana, E. D., González, A. & Sobrino, E. Changes in the high-mountain vegetation of the central Iberian Peninsula as a probable sign of global warming. Ann. Bot. 92, 273–280 (2003).PubMed 
PubMed Central 
Article 

Google Scholar 
Peñuelas, J. & Boada, M. A global change induced biome shift in the Montseny mountains (NE Spain). Glob. Change Biol. 9, 131–140 (2003).ADS 
Article 

Google Scholar 
Linares, J. C. & Tíscar, P. A. Climate change impacts and vulnerability of the southern populations of Pinus nigra subsp. salzmannii. Tree Physiol. 30, 795–806 (2010).PubMed 
Article 

Google Scholar 
Giorgi, F., Hurrell, J. W., Marinucci, M. R. & Beniston, M. Elevation dependency of the surface climate change signal: a model study. J. Clim. 10, 288–296 (1997).ADS 
Article 

Google Scholar 
Palazzi, E., Mortarini, L., Terzago, S. & Von Hardenberg, J. Elevation-dependent warming in global climate model simulations at high spatial resolution. Clim. Dyn. 52, 2685–2702 (2019).Article 

Google Scholar 
Poyatos, R., Latron, J. & Llorens, P. Land use and land cover change after agricultural abandonment. Mt. Res. Dev. 23, 362–368 (2003).Article 

Google Scholar 
Mouillot, F., Ratte, J. P., Joffre, R., Mouillot, D. & Rambal, S. Long-term forest dynamic after land abandonment in a fire prone Mediterranean landscape (central Corsica, France). Landsc. Ecol. 20, 101–112 (2005).Article 

Google Scholar 
Zellweger, F. et al. Forest microclimate dynamics drive plant responses to warming. Science 368, 772–775 (2020).ADS 
CAS 
PubMed 
Article 

Google Scholar 
Ríos-Cornejo, D., Penas, Á., Álvarez-Esteban, R. & Del Río, S. Links between teleconnection patterns and mean temperature in Spain. Theor. Appl. Climatol. 122, 1–18 (2015).ADS 
Article 

Google Scholar 
Nogués-Bravo, D., Araújo, M. B., Errea, M. P. & Martinez-Rica, J. P. Exposure of global mountain systems to climate warming during the 21st Century. Glob. Environ. Chang. 17, 420–428 (2007).Article 

Google Scholar 
Vicente-Serrano, S. M., Beguería, S., López-Moreno, J. I., El Kenawy, A. M. & Angulo-Martínez, M. Daily atmospheric circulation events and extreme precipitation risk in northeast Spain: Role of the North Atlantic Oscillation, the Western Mediterranean Oscillation, and the Mediterranean Oscillation. J. Geophys. Res. Atmos. 114, D8 (2009).Article 

Google Scholar 
Guzman-Morales, J. & Gershunov, A. Climate change suppresses Santa Ana winds of southern California and sharpens their seasonality. Geophys. Res. Lett. 46, 2772–2780. https://doi.org/10.1029/2018GL080261 (2019).ADS 
Article 

Google Scholar 
Yu, M. & Ruggieri, E. Change point analysis of global temperature records. Int. J. Climatol. 39, 3679–3688 (2019).Article 

Google Scholar 
Giorgi, F. Climate change hot-spots. Geophys. Res. Lett. 33, 08707. https://doi.org/10.1029/2006GL025734 (2006).ADS 
Article 

Google Scholar 
García, M. J. L. Recent warming in the Balearic Sea and Spanish Mediterranean coast: Towards an earlier and longer summer. Atmósfera 28, 149–160 (2015).Article 

Google Scholar 
Toreti, A., Desiato, F., Fioravanti, G. & Perconti, W. Seasonal temperatures over Italy and their relationship with low-frequency atmospheric circulation patterns. Clim. Change 99, 211–227 (2010).ADS 
Article 

Google Scholar 
Scorzini, A. R. & Leopardi, M. Precipitation and temperature trends over central Italy (Abruzzo Region): 1951–2012. Theor. Appl. Climatol. 135, 959–977 (2019).ADS 
Article 

Google Scholar 
Lee, X. et al. Observed increase in local cooling effect of deforestation at higher latitudes. Nature 479, 384–387 (2011).ADS 
CAS 
PubMed 
Article 

Google Scholar 
Juang, J.-Y., Katul, G., Siqueira, M., Stoy, P. & Novick, K. Separating the effects of albedo from eco-physiological changes on surface temperature along a successional chronosequence in the southeastern United States. Geophys. Res. Lett. 34, 21408. https://doi.org/10.1029/2007.GL03129 (2007).ADS 
Article 

Google Scholar 
Boulant, N., Kunstler, G., Rambal, S. & Lepart, J. Seed supply, drought, and grazing determine spatio-temporal patterns of recruitment for native and introduced invasive pines in grasslands. Divers. Distrib. 14, 862–874 (2008).Article 

Google Scholar 
Améztegui, A. Land-use changes as major drivers of mountain pine (Pinus uncinata Ram.) expansion in the Pyrenees. Glob. Ecol. Biogeogr. 19, 632–641 (2010).
Google Scholar 
Rambal, S. Relations entre couverts végétaux des parcours et cycle de l’eau. In L’eau des troupeaux en alpages et sur parcours: une ressource à gérer, aménager, partager (ed. Lepart, J.) 25–37 (Association Française de Pastoralisme et Cardère éditeur, 2015).
Google Scholar 
Fonderflick, J., Lepart, J., Caplat, P., Debussche, M. & Marty, P. Managing agricultural change for biodiversity conservation in a Mediterranean upland. Biol. Conserv. 143, 737–746 (2010).Article 

Google Scholar 
Abadie, J. et al. Forest recovery since 1860 in a Mediterranean region: Drivers and implications for land use and land cover spatial distribution. Landsc. Ecol. 33, 289–305 (2018).Article 

Google Scholar 
Cervera, T., Pino, J., Marull, J., Padró, R. & Tello, E. Understanding the long-term dynamics of forest transition: From deforestation to afforestation in a Mediterranean landscape (Catalonia, 1868–2005). Land Use Policy 80, 318–331 (2019).Article 

Google Scholar 
Wolpert, F., Quintas-Soriano, C. & Plieninger, T. Exploring land-use histories of tree-crop landscapes: a cross-site comparison in the Mediterranean Basin. Sustain. Sci. 15, 1267–1283 (2020).Article 

Google Scholar 
Lasanta-Martínez, T., Vicente-Serrano, S. M. & Cuadrat-Prats, J. M. Mountain Mediterranean landscape evolution caused by the abandonment of traditional primary activities: A study of the Spanish Central Pyrenees. Appl. Geogr. 25, 47–65 (2005).Article 

Google Scholar 
Malandra, F., Vitali, A., Urbinati, C., Weisberg, P. J. & Garbarino, M. Patterns and drivers of forest landscape change in the Apennines range, Italy. Reg. Environ. Change 19, 1973–1985 (2019).Article 

Google Scholar 
Zhang, Q. et al. Reforestation and surface cooling in temperate zones: Mechanisms and implications. Glob. Change Biol. 26, 3384–3401 (2020).ADS 
Article 

Google Scholar 
Rambal, S., Lacaze, B. & Winkel, T. Testing an area-weighted model for albedo or surface temperature of mixed pixels in Mediterranean woodlands. Int. J. Remote Sens. 11, 1495–1499 (1990).Article 

Google Scholar 
Luyssaert, S. et al. Land management and land-cover change have impacts of similar magnitude on surface temperature. Nat. Clim. Change 4, 389–393. https://doi.org/10.1038/nclimate2196 (2014).ADS 
Article 

Google Scholar 
Novick, K. A. & Katul, G. G. The duality of reforestation impacts on surface and air temperature. J. Geophys. Res. Biogeosci. 125, e05543 (2020).Article 

Google Scholar 
Davy, R. & Esau, I. Differences in the efficacy of climate forcings explained by variations in atmospheric boundary layer depth. Nat. Commun. 7, 1–8 (2016).Article 

Google Scholar 
Serafin, S. et al. Exchange processes in the atmospheric boundary layer over mountainous terrain. Atmosphere 9, 102. https://doi.org/10.3390/atmos9030102 (2018).ADS 
CAS 
Article 

Google Scholar 
Perugini, L. et al. Biophysical effects on temperature and precipitation due to land cover change. Environ. Res. Lett. 12, 053002 (2017).ADS 
Article 

Google Scholar 
Visbeck, M. H., Hurrell, J. W., Polvani, L. & Cullen, H. M. The North Atlantic oscillation: Past, present, and future. Proc. Natl. Acad. Sci. 98, 12876–12877 (2001).ADS 
CAS 
PubMed 
PubMed Central 
Article 

Google Scholar 
Hurrell, J. W. Decadal trends in the North Atlantic oscillation: Regional temperatures and precipitation. Science 269, 676–679 (1995).ADS 
CAS 
PubMed 
Article 

Google Scholar 
Martín, P., Sabatés, A., Lloret, J. & Martin-Vide, J. Climate modulation of fish populations: the role of the Western Mediterranean Oscillation (WeMO) in sardine (Sardina pilchardus) and anchovy (Engraulis encrasicolus) production in the north-western Mediterranean. Clim. Change 110, 925–939 (2012).ADS 
Article 

Google Scholar 
Schwingshackl, C., Hirschi, M. & Seneviratne, S. I. Global contributions of incoming radiation and land surface conditions to maximum near surface air temperature variability and trend. Geophys. Res. Lett. 45, 5034–5044 (2018).ADS 
PubMed 
PubMed Central 
Article 

Google Scholar 
Philipona, R., Behrens, K. & Ruckstuhl, C. How declining aerosols and rising greenhouse gases forced rapid warming in Europe since the 1980s. Geophys. Res. Lett. 36, L02806. https://doi.org/10.1029/2008GL036350 (2009).ADS 
Article 

Google Scholar 
Schwarz, M., Folini, D., Yang, S., Allan, R. P. & Wild, M. Changes in atmospheric shortwave absorption as important driver of dimming and brightening. Nat. Geosci. 13, 110–115 (2020).ADS 
CAS 
Article 

Google Scholar 
Norris, J. R. & Wild, M. Trends in aerosol radiative effects over Europe inferred from observed cloud cover, solar “dimming”, and solar “brightening”. J. Geophys. Res. Atmos. 112, D08214. https://doi.org/10.1029/2006JD007794 (2007).ADS 
Article 

Google Scholar 
Mateos, D. et al. Quantifying the respective roles of aerosols and clouds in the strong brightening since the early 2000s over the Iberian Peninsula. J. Geophys. Res. Atmos. 119, 10–382 (2014).Article 

Google Scholar 
Sanchez-Lorenzo, A. et al. Reassessment and update of long-term trends in downward surface shortwave radiation over Europe (1939–2012). J. Geophys. Res. Atmos. 120, 9555–9569 (2015).ADS 
Article 

Google Scholar 
Kambezidis, H. D., Kaskaoutis, D. G., Kalliampakos, G. K., Rashki, A. & Wild, M. The solar dimming/brightening effect over the Mediterranean Basin in the period 1979–2012. J. Atmos. Solar Terr. Phys. 150, 31–46 (2016).ADS 
Article 

Google Scholar 
Chiacchio, M. & Wild, M. Influence of NAO and clouds on long-term seasonal variations of surface solar radiation in Europe. J. Geophys. Res. Atmos. 115, 0022. https://doi.org/10.1029/2009JD012182 (2010).Article 

Google Scholar 
Wild, M. Decadal changes in radiative fluxes at land and ocean surfaces and their relevance for global warming. Wiley Interdiscipl. Rev. Clim. Change 7, 91–107 (2016).Article 

Google Scholar 
Held, I. M. & Soden, B. J. Water vapor feedback and global warming. Annu. Rev. Energy Environ. 25, 441–475 (2000).Article 

Google Scholar 
Dessler, A. E. & Sherwood, S. C. A matter of humidity. Science 323, 1020–1021 (2009).CAS 
PubMed 
Article 

Google Scholar 
Ruckstuhl, C., Philipona, R., Morland, J. & Ohmura, A. Observed relationship between surface specific humidity, integrated water vapor, and longwave downward radiation at different altitudes. J. Geophys. Res. Atmos. 112(D03302), 2007. https://doi.org/10.1029/2006JD007850 (2007).Article 

Google Scholar 
Parras-Berrocal, I. M. et al. The climate change signal in the Mediterranean Sea in a regionally coupled atmosphere–ocean model. Ocean Sci. 16, 743–765. https://doi.org/10.5194/os-16-743-2020 (2020).ADS 
Article 

Google Scholar 
Reale, M. et al. The regional earth system model RegCM-ES: Evaluation of the Mediterranean climate and marine biogeochemistry. J. Adv. Model. Earth Syst. 12, e001812 (2020).Article 

Google Scholar 
Sen, P. K. Estimates of the regression coefficient based on Kendall’s tau. J. Am. Stat. Assoc. 63, 1379–1389 (1968).MathSciNet 
MATH 
Article 

Google Scholar 
Kelliher, F. M., Leuning, R. & Schulze, E. D. Evaporation and canopy characteristics of coniferous forests and grasslands. Oecologia 95, 153–163 (1993).ADS 
CAS 
PubMed 
Article 

Google Scholar 
Linacre, E. T. Simpler empirical expression for actual evapotranspiration rates-a discussion. Agric. Meteorol. 11, 451–452 (1973).Article 

Google Scholar 
Jones, P. D., Jónsson, T. & Wheeler, D. Extension to the North Atlantic Oscillation using early instrumental pressure observations from Gibraltar and south-west Iceland. Int. J. Climatol. 17, 1433–1450 (1997).Article 

Google Scholar 
Palutikof, J. P. Analysis of Mediterranean climate data: measured and modelled. In Mediterranean Climate: Variability and Trends (ed. Bolle, H. J.) (Springer, 2003).
Google Scholar 
Martin-Vide, J. & Lopez-Bustins, J. A. The western Mediterranean oscillation and rainfall in the Iberian Peninsula. Int. J. Climatol. 26, 1455–1475 (2006).Article 

Google Scholar  More

We selected eleven major wheat production provinces of China for the study area, which comprise the largest winter wheat-sowing fraction: Henan, Shandong, Anhui, Jiangsu, Hebei, Hubei, Shanxi, Shaanxi, Sichuan, Xinjiang, and Gansu (Fig. 1). The wheat planting area is about 22 million ha in these provinces, accounting for more than 93% of the total wheat planting area. The total wheat production in these regions contributes more than 96% of the total wheat production in China, with more than 128 million tons in 201933.We developed a methodological framework for high-resolution AGB mapping. It mainly includes three parts: (1) Data acquisition and processing. (2) The WOFOST model parameterization and calibration. (3) Data assimilation (Fig. 1). Each part is explained in more detail below.Data acquisition and processingMeteorological dataChina Meteorological Forcing Dataset34,35 is used as weather driving data for the WOFOST model. The dataset is based on the internationally existing Princeton reanalysis data, Global Land Data Assimilation System data, Global Energy and Water Cycle Experiment-Surface Radiation Budget radiation data, and Tropical Rainfall Measuring Mission precipitation data. It is made by fusing the conventional meteorological observation data of the China Meteorological Administration. It includes seven elements: near-surface air temperature, air pressure, near-surface total humidity, wind speed, ground downward shortwave radiation, ground downward longwave radiation, and ground precipitation rate. The meteorological drive elements required for WOFOST are daily radiation, minimum temperature, maximum temperature, water vapor pressure, average wind speed, and precipitation. Details of these variables that participated in the WOFOST model can be referred to in Table S1.Soil characteristics measurements and phenology observationsSoil and phenology data were collected at 177 agricultural meteorological stations (AMS) from 2007 to 2015 (Fig. 1). Soil characteristics include soil moisture content at wilting points, field capacity, and saturation. To be consistent with the corresponding units in the crop model, the original data in weight water content was converted into volume water content through the corresponding soil bulk density measurements. Winter wheat phenology observations include the date of emergence (more than 50% of the wheat seedlings in the field show the first green leaves and reached about 2 cm), anthesis (the inner and outer glumes of the middle and upper florets of more than 50% of the wheat ears in the whole field are open, and the anthers loose powder), and maturity (more than 80% of the wheat grains turn yellow, the glumes and stems turn yellow, and only the upper first and second nodes are still slightly green). In most cases, the phenological stage “anthesis” is missing. The anthesis date was calculated by adding seven days to the observed heading date (when more than 50% of the wheat in the whole field exposes the tip of the ear from the sheath of the flag leaf).County-level yield statistics dataThe county-level yield data was collected from city statistical yearbooks of the study area from 2007 to 2015. Since most statistical yearbooks do not directly record per-unit yield data, the county-level yield was obtained by dividing the total yield and planting area. It is worth noting that all yields were calculated in units of metric kilograms per cultivated hectares (kg·ha−1).The winter wheat land cover dataWe used a winter wheat land cover product from a 1 km resolution dataset named ChinaCropArea1km36. This data was derived from GLASS leaf area index products and crop phenology from 2000 to 2015. This dataset is the base map of our data production.The MODIS LAI dataWe used the improved 8-days MODIS LAI products (i.e., 1 km) generated by Yuan et al.32 to assimilate the WOFOST model. The products applied the modified temporal-spatial filter and Savitzky-Golay filter to overcome the spatial-temporal discontinuity and inconsistence of raw MODIS LAI products, which makes them more applicable for the realm of land surface and climate modeling. The products can be accessed via the Land-Atmosphere Interaction Research Group website at Sun Yat-sen University (http://globalchange.bnu.edu.cn/research/lai).The WOFOST model parameterization and calibrationThe WOFOST model introductionThe WOFOST model was initially developed as a crop growth simulation model to evaluate the yield potential of various crops in tropical countries37. In this study, we chose the WOFOST model because the model reaches a trade-off of the complexity of the crop model and is suitable for large-scale simulations3. The WOFOST model is a typical crop growth model that explains crop growth based on underlying processes such as photosynthesis and respiration and their response to changing environmental conditions38. The WOFOST model estimates phenology, LAI, aboveground biomass, and storage organ biomass (i.e., grain yield) at a daily time step39 (Fig. 2).Fig. 2Schematic overview of the major processes implemented in WOFOST. The Astronomical module calculates day length, some variables relating to solar elevation, and the fraction of diffuse radiation.Full size imageZonal parameterizationWe first divided the study area covered by AMS into seamless Thiessen polygon zones. Each Thiessen polygon contains only a single AMS. These zones represent the whole areas where any location is closer to its associated AMS point than any other AMS point. Then, we assigned parameters to the entire zone based on the AMS data, including crop calendar (date of emergence) and soil water retention parameters (soil moisture content at wilting point, field capacity, and saturation). Besides, we also optimized two main crop parameters for controlling phenological development stages, namely TSUM1 (accumulated temperature required from emergence to anthesis) and TSUM2 (accumulated temperature required from anthesis to maturity), by minimizing the cost function of the observational and simulated date corresponding to anthesis and maturity.Parameter calibration within a single zoneWe implemented the calibration of parameters within every single zone, as illustrated in Fig. 3. We calculated the average statistical yield of each county within every single zone from 2007 to 2015, then ranked the counties in descending order and divided them into three groups, namely high, medium, and low-level yield counties, by the 33% quantile and 67% quantile of the average statistical yield. The three counties corresponding to 17% quantile, 50% quantile, and 83% quantile would be used for subsequent calibration and represent the corresponding three yield level groups. We used the statistical yields (converted to dry matter mass based on the standard moisture content of 12.5%) of the three counties for multiple years and a harvest index for each province to convert the county-level yield to AGB for calibration. The harvest index of each province was mainly estimated from the AMS’s dynamic growth records on the biomass composition of the dominant winter wheat varieties of the province and a published literature40. Besides, we collected the maximum LAI observations on all agrometeorological stations in all years in the study area, according to its histogram. We found that the histogram follows a normal distribution with a mean of 6.5 and a standard deviation of 1.5. Finally, we calibrated three sets of parameters corresponding to three yield level groups in each single zone according to the three selected counties.Fig. 3Flow chart of parameter calibration within a single zone.Full size imageWe designed a three-step calibration strategy for a specific yield level group. Firstly, as winter wheat varieties did not change significantly according to information recorded by agrometeorological stations from 2007 to 2015, we assumed the crop parameters of winter wheat remain unchanged every three years to combine three years of observational data to calibrate the parameters of the WOFOST model better. We maximized a log-likelihood function based on the maximum LAI statistics and every three-year county-level yield and AGB data mentioned to optimize selected crop parameters (see Table S2 in the Supplement Materials).The log-likelihood function was constructed as follows:$$log;{{rm{L}}}_{{rm{LAI}}}=-frac{1}{2}left[dlogleft(2pi right)+logleft(left|{Sigma }_{{rm{LAI}}}right|right)+{rm{MD}}{left({{bf{x}}}_{{rm{LAI}}};{mu }_{{rm{LAI}}},{Sigma }_{{rm{LAI}}}right)}^{2}right]$$
(1)
$$log;{{rm{L}}}_{{rm{TWSO}}}=-frac{1}{2}left[dlog(2pi )+logleft(left|{{boldsymbol{Sigma }}}_{{rm{TWSO}}}right|right)+{rm{MD}}{left({{bf{x}}}_{{rm{TWSO}}};{{boldsymbol{mu }}}_{{rm{TWSO}}},{{boldsymbol{Sigma }}}_{{rm{TWSO}}}right)}^{2}right]$$
(2)
$$log;{{rm{L}}}_{{rm{AGB}}}=-frac{1}{2}left[dlog(2pi )+logleft(left|{{boldsymbol{Sigma }}}_{{rm{AGB}}}right|right)+{rm{MD}}{left({{bf{x}}}_{{rm{AGB}}};{{boldsymbol{mu }}}_{{rm{AGB}}},{{boldsymbol{Sigma }}}_{{rm{AGB}}}right)}^{2}right]$$
(3)
$$log;{rm{L}}=log;{L}_{{rm{LAI}}}+log;{L}_{{rm{TWSO}}}+log;{L}_{{rm{AGB}}}$$
(4)
Where log L is the natural logarithm of the likelihood function, d is the dimension, that is, the number of years of joint calibration, which is set to 3 in this study xLAI is the vector composed of the maximum value of the 3-year LAI simulated by WOFOST, μLAI and ΣLAI are the mean vector and error covariance matrix of maximum LAI based on observation statistics. The annual maximum LAI was assumed to be independent, and the mean and standard deviation for each year was set the same as the result of Fig. 3. Similarly, xTWSO and xAGB are the yield vector and AGB vector at maturity of 3 years simulated by WOFOST, and μTWSO, μAGB are their corresponding county-level statistic vector, ΣTWSO and ΣAGB are their corresponding error covariance matrix. In this study, we assumed that the annual yield or AGB was independent, and their corresponding standard deviation was 10% of their statistical value. |Σ| is the determinant of Σ. The expression ({rm{MD}}{({bf{x}};{boldsymbol{mu }},{boldsymbol{Sigma }})}^{2}={({bf{x}}-{boldsymbol{mu }})}^{{rm{T}}}{{boldsymbol{Sigma }}}^{-1}({bf{x}}-{boldsymbol{mu }})), where MD is the Mahalanobis distance between the point x and the mean vector μ.Secondly, we optimized the inter-annual irrigation. We optimized two parameters every year: the critical value of soil moisture (denoted as SMc) and the amount of irrigation (denoted as V). When the soil moisture simulated by WOFOST is lower than SMc, an irrigation event will be triggered, and the irrigation amount is V cm. In this study, we combined three-year data for calibration with six parameters for optimization. The optimization strategy is the same as the previous step by maximizing the log-likelihood function. Finally, we fixed the optimized irrigation parameters and repeated the first step to calibrate the selected crop parameters and obtain the final optimal parameters.Data assimilationConsidering that MODIS LAI is relatively low compared to the actual LAI of winter wheat41, we select a weak-constraint cost function based on the least square of normalized observational and simulated LAI as shown in Eq. (5), which is assimilating the trend information of MODIS LAI into the crop growth model.$$J={sum }_{{rm{t}}=1}^{{rm{n}}}{left(frac{{{rm{LAI}}}_{{rm{MODIS}}}^{{rm{t}}}-{{rm{LAI}}}_{{rm{MODIS}}}^{min}}{{{rm{LAI}}}_{{rm{MODIS}}}^{max}-{{rm{LAI}}}_{{rm{MODIS}}}^{min}}-frac{{{rm{LAI}}}_{{rm{WOFOS}}}^{{rm{t}}}-{{rm{LAI}}}_{{rm{WOFOS}}}^{min}}{{{rm{LAI}}}_{{rm{WOFOS}}}^{max}-{{rm{LAI}}}_{{rm{WOFOS}}}^{min}}right)}^{2}$$
(5)
Where ({{rm{LAI}}}_{{rm{MODIS}}}^{{rm{t}}}) and .. are MODIS LAI and WOFOST simulated LAI of time t. ({{rm{LAI}}}_{{rm{MODIS}}}^{max}) and ({{rm{LAI}}}_{{rm{WOFOS}}}^{max}) are maximum of MODIS LAI and WOFOST simulated LAI. ({{rm{LAI}}}_{{rm{MODIS}}}^{min}) and ({{rm{LAI}}}_{{rm{WOFOS}}}^{min}) are minimum of MODIS LAI and WOFOST simulated LAI. J is the value of the cost function.We reinitialize the day of emergence (IDEM), the life span of leaves growing at 35 °C (SPAN), and thermal time from emergence to anthesis (TSUM1) in the WOFOST model on each 1 km winter wheat pixel according to cost function between WOFOST LAI and MODIS LAI. Besides, we applied the Subplex algorithm from the NLOPT library (https://github.com/stevengj/nlopt) for parameter optimization. More

Virus propagation and enumerationEchovirus-11 (E11, Gregory strain, ATCC VR737) and Coxsackievirus-A9 (CVA9, environmental strain from sewage, kindly provided by the Finnish National Institute for Health and Welfare) stocks were produced by infecting sub-confluent monolayers of BGMK cells as described previously [7]. Viruses were released from infected cells by freezing and thawing the culture flasks three times. To eliminate cell debris, the suspensions were centrifuged at 3000 × g for 5 min. Each stock solution was stored at −20 °C until use. Infectious virus concentrations were enumerated by a most probable number (MPN) infectivity assay as described in the Supplementary Information. The assay limit of detection (LoD), defined as the concentration corresponding to one positive cytopathic effect in the lowest dilution of the MPN assay under the experimental conditions used, corresponding to 2 MPN/mL.Inactivation of enteroviruses by bacterial consortia from lake waterTo study the inactivation of CVA9 and E11 by a bacterial consortium from lake water, four surface water samples were collected from Lake Geneva (Ecublens, Switzerland) during the summer 2021. Each sampling event was conducted on warm and sunny days, to minimize biological variation. Immediately after sampling, large particles of the sample were removed by filtering 500 mL of water on a 8 μm nitrocellulose filter membrane (Merck Millipore, Cork, Ireland). The sample was then filtered through a 0.8 μm nitrocellulose filter membrane (Merck Millipore) to remove large microorganisms such as protists. The resulting water sample corresponds to the bacterial fraction used to study virus inactivation.For inactivation experiments, each virus was spiked into individual 1 mL aliquots of fractionated lake water to a final concentration of 106 MPN/mL, and samples were incubated for 48 h at 30 °C without shaking. Duplicate experiments were conducted for each virus and each lake water sample. Experiments to control for thermal inactivation were conducted using the same procedure but by replacing the fractionated lake water with sterile milliQ water. Viral infectivity at times 0 h and 48 h was determined by MPN as described above. Virus decay was calculated as log10 (C/C0), where C is the residual titer after 48 h of incubation, and C0 is the initial titer. The experimental LoD was approximately 5-log10.These same experiments were conducted for three new water samples in the presence of four protease inhibitors with the following final concentrations: E64—10 μM (E3132, Sigma–Aldrich, Saint-Louis, MO, USA), GM6001—4 μM (CC1010, Sigma–Aldrich), Chymostatin—100 μM (C7268, Sigma–Aldrich), and PMSF—100 μM (P7626, Sigma–Aldrich). Each inhibitor was added to 1 mL of fractionated lake water, vortexed for 30 seconds, and incubated at room temperature for 15 min, before adding the two viral strains under the same conditions as described above.Bacterial isolation, cultivation, and storageBacteria were isolated from two water samples from Lake Geneva’s Ecublens beach, taken in November 2019 (Fall, 89 isolates) and May 2020 (Spring, 47 isolates). Bacteria recovery was performed on R2A agar plate (BD Difco, Franklin Lakes, NJ, USA) as described previously [15]. Briefly, successive dilutions from 10−1 to 10−5 were carried out in sterile water for each sample. For each dilution, a volume of 1 mL was deposited on three separate R2A plates, before being incubated at 22, 30, and 37 °C. After 5 days of incubation, each colony was picked and enriched on a new R2A plate. To ensure purity, each isolate was successively plated five times on R2A plate and incubated at the same temperature as the initial isolation. Each purified isolate was cryopreserved in R2A / 20% glycerol at −80 °C. The isolates were named based on the water body (Lake (L)), isolation temperature, and the isolation order (L-T°C-number).Bacterial identificationThe identification of each isolate was performed by 16 S rRNA gene sequencing using the pair of primers 27 F (5’- AGA GTT TGA TCM TGG CTC AG- 3’, Microsynth AG, Balgach, Switzerland) / 786 R (5’- CTA CCA GGG TAT CTA ATC – 3’, Microsynth AG), following a methodology previously described [15]. The thermocycling conditions and the purification of PCR products are described in the Supplementary Information. The complete list of isolated bacteria and associated accession numbers is given in Supplementary Table 1.Phylogenetic inference and metadata visualizationThe consensus from 16 S rRNA gene sequences of the 136 isolates was aligned using the MUSCLE algorithm [16]. The phylogenetic analysis of 566 bp aligned sequences from the V2-V4 16 S rRNA gene regions (Positions: 152–717) was performed using Molecular Evolutionary Genetics Analysis X software [17]. Phylogeny was inferred by maximum likelihood, with 1000 bootstrap iterations to test the robustness of the nodes. The resulting tree was uploaded and formatted using iTOL [18].Virus incubation with bacterial isolatesFor the preparation of the bacteria before co-incubation, each one was first cultured on R2A agar for 48 h at their initial isolation temperature. Overnight suspensions of each bacterial isolate were grown in R2A broth at room temperature under constant agitation (180 rpm). For co-incubation experiments, 200 μL of each bacterial suspension were mixed with 100 μL of a 105 MPN/mL stock of E11 or CVA9. Then, each condition was supplemented with 600 μL of R2A broth. Incubation was carried out for 96 h at room temperature, without shaking. At the end of the co-incubation, each tube was centrifuged for 15 min at 9000 × g (4 °C) to eliminate bacteria, and the residual infectious viral titer was enumerated by MPN assay as described above [7]. Each co-incubation experiment was carried out in triplicate. Control experiments were performed under the same conditions but using sterile R2A. Virus decay was quantified as log10 (Cexp/Cctrl), where Cexp is the residual titer after a co-incubation for 96 h, and Cctrl is the titer after incubation of the virus in sterile R2A for 96 h. The experimental LoD was 3-log10.Protease activity measurement using casein and gelatin agar platesCasein agar was prepared as follows: 20 g of skim milk (BD Difco), supplemented with 1 g glucose were reconstituted with 200 mL of distilled water. Likewise, a 10% bacteriological agar solution was prepared in a final volume of 200 mL. Finally, a solution consisting of 0.8% NaCl, 0.02% KCl, 0.144% Na2HPO4, and 0.024% KH2PO4 was reconstituted in 600 mL of water. All solutions were autoclaved for 15 min at 110 °C. The solutions were mixed, and 25 mL were poured into each Petri dish. Gelatin agar was composed of 0.4% peptone, 0.1% yeast extract, 1.5% gelatin and 1.5% bacteriological agar. The mixture was autoclaved 15 min at 120 °C, and 25 mL of medium was poured into each Petri dish.For each isolate, an overnight suspension was performed in R2A broth at room temperature, before spotting 15 μL of each suspension at the center of both gelatin and casein agar plates. Each plate was incubated at 22, 30, or 37 °C for 72 h, depending on the initial isolation temperature of the bacteria. Casein-degrading activity (cas), which is exerted by many different protease classes, and gelatin-degrading activity (gel), which is mostly caused by MMPs, were revealed by a hydrolysis halo around the producing bacteria. Hydrolysis diameters were measured in millimeters (mm) to report the extent of the proteolytic effect of each strain on both substrates.Protease activity quantification in cell-free supernatantUsing the same bacterial suspensions as for bacterial/virus co-incubation, 200 μL of each suspension was inoculated into 600 μL of R2A broth and incubated without shaking for 96 h at room temperature. Each culture was centrifuged for 15 min at 9000 × g at 4 °C. The resulting cell-free supernatants (CFS) were stored at −20 °C until use. For each CFS, protease activity was measured using the Protease Activity Assay Kit (ab112152, Abcam, Cambridge, UK), which measures general protease activity (pgen) except MMPs, and the MMP Activity Assay Kit (ab112146, Abcam), which selectively measures MMP activity (mmp). Briefly, for the Protease Activity Assay kit, 50 μL of the substrate was added into each well of a dark-bottom plate containing 50 μL of each CFS. Standard trypsin provided by the kit was used as a positive control. For the MMP Activity Assay kit, 50 μL of each CFS was incubated with 50 μL of 2 mM APMA for 3 h at 37 °C, prior to the activity test. Collagenase I (C0130, Sigma–Aldrich) was used as a positive control. R2A broth was used as a negative control for each assay. Protease activity was measured at time 0 and after 60 min, using a Synergy MX fluorescence reader (BioTek). The excitation and emission wavelengths were set to 485 and 530 nm, respectively. The emitted fluorescence, generated by proteolytic cleavage of the substrate of each kit, was calculated as follows: ∆RFU = RFU (60 min) − RFU (0 min). Proteolytic activity was calculated in mmol/min/μL based on the emitted fluorescence measured for trypsin and collagenase I at known proteolytic activities.Data analysisStatistical analyses to compare inactivation data were performed by one-way t-test or one-way ANOVA with Dunnett’s post-hoc test in GraphPad Prism v.9. An alpha value of 0.05 was used as a threshold for statistical significance. For each dataset we confirmed that data were normally distributed.To analyze a potential correlation between protease activity and viral decay, the decay values for each virus strain was related to the four protease activity tests of this study using a scatterplot combined with a Kernel density estimation. The analyses were performed with R v.3.6.1 using the SmoothScatter function of the R Base package.A Left-Censored Tobit model (CTM) with mixed effects was chosen to investigate interactions between protease activity and the decay measured for each virus strain. Briefly, the CTM with mixed effect was chosen for three reasons: (1) The protocol used to measure viral decay had a limit of quantification of −3-log10, and 152 measurement points reached the detection limit, requiring the use of this value as the left-censored value of the model; (2) The two virus strains used in the study showed distinct responses after exposure to environmental bacteria, preventing the use of a multiple linear regression model; (3) Among biological replicates of co-incubation experiments, inactivation variability was observed, suggesting the concomitant action of random biological effects (e.g., production of other compounds than proteases by bacteria, or differences in protease production rate between replicates for each bacterial isolate). The resulting statistical model was then formulated as follows:$$log left( {frac{{C_{{{{{{mathrm{exp}}}}}}}}}{{C_{{{{{{mathrm{ctrl}}}}}}}}}} right) = ; beta _0 + beta _1;{rm I}_{{{{{{{{mathrm{virus}}}}}}}}_i = 2} + beta _2sqrt {left[ {pgen} right]_i} + beta _3sqrt {left[ {mmp} right]_i} + beta _4sqrt {left[ {cas} right]_i} \ + beta _5sqrt {left[ {gel} right]_i} + beta _6I_{{{{{{{{mathrm{virus}}}}}}}}_i = 2}sqrt {left[ {pgen} right]_i} + beta _7I_{{{{{{{{mathrm{virus}}}}}}}}_i = 2}sqrt {left[ {mmp} right]_i} \ + beta _8I_{{{{{{{{mathrm{virus}}}}}}}}_i = 2}sqrt {left[ {cas} right]_i} + beta _9I_{{{{{{{{mathrm{virus}}}}}}}}_i = 2}sqrt {left[ {gel} right]_i} + alpha _{{{{{{{{mathrm{id}}}}}}}}_i} + varepsilon _i$$$${{{mbox{where}}}}; log left( {frac{{C_{{{{{{mathrm{exp}}}}}}}}}{{C_{{{{{{mathrm{ctrl}}}}}}}}}} right) = left{ {begin{array}{*{20}{c}} { – 3} & {{{{{{{{mathrm{if}}}}}}}};{{{{{{{mathrm{log}}}}}}}}left( {frac{{C_{{{{{{mathrm{exp}}}}}}}}}{{C_{{{{{{mathrm{ctrl}}}}}}}}}} right) le – 3} \ {{{{{{{{mathrm{log}}}}}}}}left( {frac{{C_{{{{{{mathrm{exp}}}}}}}}}{{C_{{{{{{mathrm{ctrl}}}}}}}}}} right)} & {{{{{{{{mathrm{otherwise}}}}}}}}} end{array}} right.$$$$alpha _{{{{{{{{mathrm{id}}}}}}}}_i}sim {{{{{{{mathrm{i}}}}}}}}.{{{{{{{mathrm{i}}}}}}}}.;{{{{{{{mathrm{d}}}}}}}}.;{rm N}left( {0,;sigma _{{{{{{{{mathrm{id}}}}}}}}}^2} right)$$$${{{{{{{mathrm{for}}}}}}}};i in left{ {1,2, ldots } right}$$for which β0 defines the model intercept, (beta _1{rm I}_{{{{{{{{mathrm{virus}}}}}}}}_i = 2}) corresponds to the main effect of the virus factor on the viral decay, (beta _2,;beta _3,;beta _4,;{{{{{{{mathrm{and}}}}}}}};beta _5) corresponds to the main effects of the different protease activity measurements on viral decay, (beta _6I_{{{{{{{{mathrm{virus}}}}}}}}_i = 2},;beta _7I_{{{{{{{{mathrm{virus}}}}}}}}_i = 2},;beta _8I_{{{{{{{{mathrm{virus}}}}}}}}_i = 2},{{{{{{{mathrm{and}}}}}}}};beta _9I_{{{{{{{{mathrm{virus}}}}}}}}_i = 2}) corresponds to the interaction effects between each of these variables and the viral decay, (alpha _{{{{{{{{mathrm{id}}}}}}}}_i}) corresponds to the mixed effect of the model and (varepsilon _i) corresponds to the error term of the model. The selection of the model is further detailed in the Supplementary Information (Supplementary Material and Figs. S1 and S2).The full dataset included in the correlation analysis and the CTM is provided in Supplementary Table 2. A description of the variables used is given in the Supplementary Information. The dataset was analyzed using the censReg package in R [19]. The R code is given in the Supplementary Information. More

Data collectionThe acoustic recordings were collected during 2017 off the Changhua coast (24° 4.283 N/120° 19.102 E) (Fig. 5) by deploying a passive acoustic monitoring (PAM) device from Wildlife Acoustics, which was an SM3M recorder moored at a depth of 18–20 m. The hydrophone recorded continuously with a sampling frequency of 48 kHz and a sensitivity of −164.2 dB re:1 v/µPa. The acoustic files were recorded in the.WAV format with a duration of 60 minutes. The hydrophone setup prior to deployment is shown in Fig. 6. Table 2 contains the details for the monitoring period with the corresponding season and the number of hours of recordings each season used in this study. Studies have shown that the presence of seasonal chorusing at this monitoring site in the frequency range of 500–2500 Hz is caused by two types of chorusing15,38, with chorusing starting in early spring, reaching a peak in summer, and starting to decline late autumn, and silencing in winter6. Previous studies6,15,38 at this monitoring site have derived the details of two types of fish calls responsible for chorusing (Type 1 and Type 2); Supplementary Fig. 1 shows the spectrogram, waveform, and power spectrum density of the individual calls. Supplementary Table 1 tabulated are the acoustic features of the two call types. The monitoring region, Changhua, lies in the Eastern Taiwan Strait (ETS). The ETS is ~350 km in length and ~180 km wide64. The ETS experiences diverse oceanographic and climatic variations influenced by monsoons in summer and winter65 and extreme events caused by tropical storms, typhoons in summer, and wind/cold bursts occurring in winter66,67,68.Fig. 5: Study area located off the Taiwan Strait.Map of the Changhua coast located in Taiwan Strait, Taiwan depicting the deployed passive acoustic monitoring recorder at site A1. The map was produced in Matlab 9.11 (The Mathworks, Natick, MA; http://www.mathworks.com/) using mapping toolbox function geobasemap(). Full global basemap composed of high-resolution satellite imagery hosted by Esri (https://www.esri.com/).Full size imageFig. 6: Setup of the SM3M submersible recorder.SM3M recorder fastened to the steel frame (length and breadth = 1.22 m, height = 0.52 m) with plastic cable zip ties prior to deployment.Full size imageTable 2 Passive acoustic monitoring device specifications and monitoring duration during different seasons.Full size tableAcoustic data analysisThe acoustic data were analyzed using the PAMGuide toolbox in Matlab60. The seasonal spectrograms were computed with an FFT size of 1024 points and a 1 s time segment averaged to a 60 s resolution. The sound pressure levels (SPL) were computed in the frequency band of 500–3500 Hz and programmed to provide a single value every hour, thus resulting in 984, 1344, and 1440 data points in spring, summer, and winter, respectively (Supplementary Data 1).Determining the regularity and complexity with the complexity-entropy planeThe complexity-entropy plane was utilized in this study to quantify the structural statistical complexity and the regularity in the hourly acoustical and seasonal SPL time series data. The C-H plane is a 2D plane representation of the permutation entropy on the horizontal axis that quantifies the regularity, and the vertical axis is represented by the statistical complexity quantifying the correlation structure in the temporal series.For a given time series ({{x(t)}}_{t=1}^{N}), the N’ ≡ N − (m − 1) the values of the vectors for the length m  > 1 are ranked as$${X}_{s}=left({x}_{s-(m-1)},{x}_{s-(m-2)},ldots ,{x}_{s}right),s=1,ldots ,,{N}^{{prime} }$$
(1)
Within each vector, the values are reordered in the ascending order of their amplitude, yielding the set of ordering symbols (patterns) ({r}_{0},{r}_{1},ldots ,{r}_{m-1}) such that$${x}_{s-{r}_{0}}le {x}_{s-{r}_{1}}le ..,..le {x}_{s-{r}_{(m-1)}}$$
(2)
This symbolization scheme was introduced by Bandt and Pompe69. The scheme performs the local ordering of a time series to construct a probability mass function (PMF) of the ordinal patterns of the vector. The corresponding vectors (pi ={r}_{0},{r}_{1},ldots ,{r}_{(m-1)}) may presume any of the m! possible permutations of the set ({{{{{mathrm{0,1}}}}},ldots ,m-1}) and symbolically represent the original vector. For instance, for a given time series {9, 4, 5, 6, 1,…} with length m = 3, provides 3! different order patterns with six mutually exclusive permutation symbols are considered. The first three-dimensional vector is (9, 4, 5), following the Eq. (1), this vector corresponds to ((,{x}_{s-2},{x}_{s-1},{x}_{s})). According to Eq. (2), it yields ({x}_{s-1}le {x}_{s}le {x}_{s-2}). Then, the ordinal pattern satisfying the Eq. (2) will be (1, 0, 2). The second 3-dimensional vector is (4, 5, 6), and (2, 1, 0) will be its associated permutation, and so on.The permutation entropy (H) of order m ≥ 2 is defined as the Shannon entropy of the Brandt-Pompe probability distribution p(π)69$$Hleft(mright)=,-{mathop{sum}limits _{{pi }}}pleft(pi right){{{log }}}_{2}p(pi )$$
(3)
where ({pi }) represents the summation over all possible m! permutations of order m, (p(pi )) is the relative frequency of each permutation (pi), and the binary logarithm (base of 2) is evaluated to quantify the entropy in bits. H(m) attains the maximum ({{log }}m!) for (p(pi )=1/m!). Then the normalized Shannon entropy is given by$$0le H(m)/{{{{{rm{ln}}}}}},m!le 1$$
(4)
where the lower bound H = 0 corresponds to more predictable signals with fewer fluctuations, an either strictly increasing or decreasing series (with a single permutation), and the upper bound H = 1 corresponds to an unpredictable random series for which all the m! possible permutations are equiprobable. Thus, H quantifies the degree of disorder inherent in the time series. The choice of the pattern length m is essential for calculating the appropriate probability distribution, particularly for m, which determines the number of accessible states given by m!70,71. As a rule of thumb, the length T of the time series must satisfy the condition T (gg) m! in order to obtain reliable statistics, and for practical purposes, Bandt and Pompe suggested choosing m = 3,…,7 69.The statistical complexity measure is defined with the product form as a function of the Bandt and Pompe probability distribution P associated with the time series. (Cleft[Pright]) is represented as33$$Cleft[Pright]=frac{J[P,U]}{{J}_{{max }}}{H}_{s}[P]$$
(5)
where ({H}_{s}left[Pright]=Hleft[Pright]/{{log }}m!) is the normalized permutation entropy. (J[P,U]) is the Jensen divergence$$Jleft[P,Uright]=left{Hleft[frac{P+U}{2}right]-frac{H[P]}{2}-frac{H[U]}{2}right}$$
(6)
which quantifies the difference between the uniform distributions U and P, and ({J}_{{max }})is the maximum possible value of (Jleft[P,Uright]) that is obtained from one of the components of P = 1, with all the other components being zero:$$Jleft[P,Uright]=-frac{1}{2}left[frac{m!+1}{m!}{{log }}left(m!+1right)-2{{log }}left(2m!right)+{{log }}(m!)right]$$
(7)
For each value of the normalized permutation entropy (0le {H}_{s}[P]le 1) there is a corresponding range of possible statistical complexity (Cleft[Pright]) values. Thus, the upper (({C}_{{max }})) and lower ((C_{{min }})) complexity bounds in the C-H plane are formed. The periodic sequences such as sine and series with increasing and decreasing (with ({H}_{s}[P]=0)) and completely random series such as white noise (for which (Jleft[P,Uright]=0) and ({H}_{s}[P]=1)) will have zero complexity. Furthermore, for each given value of the (0le {H}_{s}[P]le 1), there exists a range of possible values of the statistical complexity, ({C}_{{min }}le C[P]le {C}_{{max }}). The procedure for evaluating the complexity bounds ({C}_{{min }}) and ({C}_{{max }}) is given in Martin et al.72. We utilized the R package ‘statcomp’73 to evaluate the statistical complexity (C) and the permutation entropy (H) using the command global-complexity() for the order m = 6, and the command limit_curves(m, fun = ‘min/max’) was utilized to evaluate the complexity boundaries ({C}_{{min }}) and ({C}_{{max }}). In this study, we constructed two C-H planes: (1) C and H was computed for each hourly acoustic file during different seasons. The resulting lengths of C and H during spring, summer, and autumn-winter are similar to the number of hours in the particular season (Table 2). (2) C and H was computed every 4–5 days for the seasonal SPL. The resulting length of C and H obtained during spring was 9 points (each value of C and H for every 109 h), and in summer and autumn-winter was 12 points (each value of C and H for every 112 and 120 h).Determining predictability and dynamics (linear/nonlinear) using EDMIn this study, we utilized EDM to quantify the predictability (forecasting) and distinguish between the linear stochastic and nonlinear dynamics in the seasonal soundscape SPL. EDM involves phase-space reconstruction via delay coordinate embeddings to make forecasts and to determine the ‘predictability portrait’ of time series data40. The first step in EDM is to determine the optimal embedding dimension (E), and this is obtained using a method based on simplex projection41. The simplex projection is carried out by dividing the dataset into two equal parts, of which the first part is called the library and the other part is called the target. The library set is used to build a series of non-parametric models (known as predictors) for the one step ahead predictions for the E varying between 1 and 10. Then the model’s accuracies are determined when the model is applied to the target dataset and the prediction skill (⍴) for the actual and predicted datasets is measured. The embedding dimension with the highest predictive skill is the optimal E.For the appropriate optimal E chosen, the predictability profile of the time series data is obtained by evaluating ⍴ at Tp = 1, 2, 3, … steps ahead. The flat prediction profile of the ⍴–Tp curve indicates that the time series is purely random (low ⍴) or regularly oscillating (high ⍴). In contrast, a decreasing ⍴ as Tp increases may reject the possibility of an underlying uncorrelated stochastic process and indicate the presence of low-dimensional deterministic dynamics. However, the concern with the predictability profile is that it may exhibit predictability even if time series are purely stochastic (such as autocorrelated red noise). Hence, a nonlinear test can be performed by using S-maps (sequential locally weighted global linear maps) to distinguish between linear stochastic and nonlinear dynamics in the time series dataset by fitting a local linear map. S-maps similar to simplex projects provide the forecasts in phase-space by quantifying the degree to which points are weighted when fitting the local linear map, which is given by the nonlinear localization parameter θ. When θ = 0, the entire library set will exhibit equal weights regardless of the target prediction, which mathematically resembles the model of a linear autoregressive process. In contrast, if θ  > 0, the forecasts of the library provided by the S-map depend on the local state of the target prediction, thus producing large weights, and the unique local fittings can vary in phase-space to incorporate nonlinear behavior. Consequently, if the (⍴–θ) dynamics profile shows the highest ⍴ at θ = 0 that is reduced as θ increases, it represents linear stochastic dynamics. If the ⍴ achieves the highest value at θ  > 0, then the dynamics are represented by nonlinear dynamics.In this study, the R package “rEDM”74 was used to evaluate the optimal E, prediction profile (⍴–Tp), and dynamics profile (⍴–θ) for the seasonal SPL dataset. While evaluating these entities, the data points are equally into two parts, and sequentially the first half is chosen as the library set and the other as the target set. The length of the library and the target set for spring, summer, and autumn-winter are 492, 672, and 720. The command EmbedDimension() was used to determine the forecast skill for the E ranging from 1 to 10 and the optimal E with the highest forecast skill (Supplementary Fig. 2) was chosen. In this study, we found that for all seasons, the optimal E was 2. The (⍴–Tp) was evaluated for Tp varying between 1 and 100 using the command PredictInterval() and the (⍴–θ) was evaluated using the command PredictNonlinear() for θ = 0, 0.0001, 0.0003, 0.001, 0.003, 0.01, 0.03, 0.1, 0.3, 0.5,0.75, 1.0, 1.5, 2, and 3 to 20.StatisticsThe nonparametric Kruskal–Wallis test, followed by post hoc Bonferroni’s multiple comparisons, was used to test differences in the seasonal H and C that were obtained directly from the hourly acoustic data during chorusing hours, as well as the H and C obtained for the seasonal SPL and the seasonal forecast skill. The statistical calculations were performed using the R package “agricolae” 75. More

IPCC (2013) Climate Change: The Physical Science Basis. Contribution of Working Group I to the Fifth Assessment Report of the Intergovernmental Panel on Climate Change [Stocker TF, Qin D Qin, Plattner G-K, Tignor M, Allen SK, Boschung J, Nauels A, Xia Y, Bex V & Midgley PM (eds.)]. Cambridge University Press, Cambridge, United Kingdom and New York, NY, USA, 1535 pp (2013).Graham, L.P. Projections of Future Anthropogenic Climate Change [in:] Assessment of Climate Change for the Baltic Sea Basin. Regional Climate Studies. Bolle H.J., Menenti M., Rasool I. Series Editors Springer-Verlag Berlin Heidelberg s.133–220 (2008).Früchtenich, E., Bock, J., Feucht, V., Früchtenich W. Reactions of three European oak species ( Q. robur, Q. petraea and Q. ilex ) to repetitive summer drought in sandy soil. Trees, Forests and People 5: 100093 (2021).Gray, S. B. & Brady, S. M. Plant developmental responses to climate change. Dev. Biol. 419, 64–77 (2016).CAS 
Article 

Google Scholar 
Willliams, A. & De Vries, F. T. Plant root exudation under drought: Implications for ecosystem functioning. New Phytol. 225, 1899–1905 (2019).Article 

Google Scholar 
Canarini, A., Merchant, A. & Dijkstra, F. A. Drought effects on Helianthus annuus and Glycine max metabolites: From phloem to root exudates. Rhizosphere 2, 85–97 (2016).Article 

Google Scholar 
De Vries, F. T. et al. Changes in root-exudate-induced respiration reveal a novel mechanism through which drought affects ecosystem carbon cycling. New Phytol. 224, 132–145 (2019).Article 

Google Scholar 
Phillips, R. P., Finzi, A. C. & Bernhardt, E. S. Enhanced root exudation indu ces microbial feedbacks to N cycling in a pine forest under long-term CO2 fumigation. Ecol. Lett. 14, 187–194 (2011).Article 

Google Scholar 
Meier, I. C. et al. Root exudation of mature beech forests across a nutrient availability gradient: The role of root morphology and fungal activity. New Phytol. 226, 583–594 (2020).CAS 
Article 

Google Scholar 
Gianfreda, L. Enzymes of importance to rhizosphere processes. J. Soil Sci. Plant Nutr. 15, 283–306 (2015).
Google Scholar 
Małek, S., Ważny, R., Błońska, E., Jasik, M. & Lasota, J. Soil fungal diversity and biological activity as indicators of fertilization strategies in a forest ecosystem after spruce disintegration in the Karpaty Mountains. Sci. Total Environ. 751, 142335 (2021).ADS 
Article 

Google Scholar 
Zuccarini, P., Asensio, D., Ogaya, R., Sardans, J. & Penuelas, J. Effects of seasonal and decadal warming on soil enzymatic activity in a P-deficient Mediterranean shrubland. Glob. Change Biol. 26, 3698–3714 (2019).ADS 
Article 

Google Scholar 
Sing, S. et al. Soil organic carbon cycling in response to simulated soil moisture variation under field conditions. Sci. Rep. 11, 10841 (2021).ADS 
Article 

Google Scholar 
Sardans, J. & Penuelas, J. Drought decreases soil enzyme activity in a Mediterranean Quercus ilex L. forest. Soil Biol. Biochem. 37, 455–461 (2005).CAS 
Article 

Google Scholar 
Czúcz, B., Gálhidy, L. & Mátyás, C. Present and forecasted xeric climatic limits of beech and sessile oak distribution at low altitudes in Central Europe. Ann. For. Sci. 68, 99–108. https://doi.org/10.1007/s13595-011-0011-4 (2011).Article 

Google Scholar 
Sáenz-Romero, C. et al. Adaptive and plastic responses of Quercus petraea populations to climate across Europe. Glob. Change Biol. 23, 2831–2847 (2018).ADS 
Article 

Google Scholar 
Regulation of the Minister of the Environment. Detailed requirements for the forest reprudactive material (Dz. U. Nr 31, poz. 272) (in Polish) (2004).Phillips, R. P., Erlitz, Y., Bier, R. & Bernhardt, E. S. New approach for capturing soluble root exudates in forest soils. Funct. Ecol. 22, 990–999. https://doi.org/10.1111/j.1365-2435.2008.01495.x (2008).Article 

Google Scholar 
Ostonen, I., Lõhmus, K. & Lasn, R. The role of soil conditions in fine root ecomorphology in Norway spruce (Picea abies (L.) Karst.). Plant Soil 208, 283–292 (1999).CAS 
Article 

Google Scholar 
Pritsch, K. et al. A rapid and highly sensitive method for measuring enzyme activities in single mycorrhizal tips using 4-methylumbelliferone-labelled fluorogenic substrates in a microplate system. J. Microbiol. Methods 58, 233–241 (2004).CAS 
Article 

Google Scholar 
Sanaullah, M., Razavi, B. S., Blagodatskaya, E. & Kuzyakov, Y. Spatial distribution and catalytic mechanisms of β-glucosidase activity at the root-soil interface. Biol. Fert. Soils 52, 505–514 (2016).CAS 
Article 

Google Scholar 
R Core Team (2021). R: A language and environment for statistical computing. R Foundation for Statistical Computing, Vienna, Austria. URL http://www.R-project.org/.Hartmann, H. Will a 385million year-struggle for light become a struggle for water and for carbon?–how trees may cope with more frequent climate change-type drought events. Glob. Change Biol. 17, 642–655 (2011).ADS 
Article 

Google Scholar 
Brunner, I., Herzog, C., Dawes, M. A., Arend, M. & Sperisen, C. How tree roots respond to drought. Front. Plant Sci. 6, 547 (2015).Article 

Google Scholar 
Markesteijn, L. & Poorter, L. Seedling root morphology and biomass allocation of 62 tropical tree species in relation to drought- and shade-tolerance. J. Ecol. 97, 311–325 (2009).Article 

Google Scholar 
Poorter, L. & Markesteijn, L. Seedling Traits Determine Drought Tolerance of Tropical Tree Species. Biotropica 40, 321–331 (2008).Article 

Google Scholar 
Ostonen, I. et al. Specific root length as an indicator of environmental change. Plant Biosyst. 141, 426–442 (2007).Article 

Google Scholar 
Lozano, Y. M., Aguilar-Triqueros, C. A., Flaig, I. C. & Rillig, M. C. Root trait responses to drought are more heterogeneous than leaf trait responses. Funct. Ecol. 34, 2224–2235 (2020).Article 

Google Scholar 
De Vries, F. T., Brown, C. & Stevens, C. J. Grassland species root response to drought: consequences for soil carbon and nitrogen availability. Plant Soil 409, 297–312 (2016).Article 

Google Scholar 
Sell, M. et al. Responses of fine root exudation, respiration and morphology in three early successional ree species to increased air humidity and different soil nitrogen sources. Tree Physiol. 42, 557–569. https://doi.org/10.1093/treephys/tpab118 (2021).Article 

Google Scholar 
Karlowsky, S. et al. Drought-induced accumulation of root exudates supports post-drought recovery of microbes in mountain grassland. Front. Plant Sci. 9, 1593. https://doi.org/10.3389/fpls.2018.01593 (2018).Article 
PubMed 
PubMed Central 

Google Scholar 
Fuchslueger, L., Bahn, M., Fritz, K., Hasibeder, R. & Richter, A. Experimental drought reduces the transfer of recently fixed plant carbon to soil microbes and alters the bacterial community composition in a mountain meadow. New Phytol. 201, 916–927. https://doi.org/10.1111/nph.12569 (2014).CAS 
Article 
PubMed 

Google Scholar 
Gargallo-Garriga, A. et al. Root exudate metabolomes change under drought and show limited capacity for recovery. Sci. Rep. 8, 12696. https://doi.org/10.1038/s41598-018-30150-0 (2018).ADS 
CAS 
Article 
PubMed 
PubMed Central 

Google Scholar 
Zhang, X., Dippold, M. A., Kuzyakov, Y. & Razavi, B. S. Spatial pattern of enzyme activities depends on root exudate composition. Soil Biol. Biochem. 133, 83–93. https://doi.org/10.1016/j.soilbio.2019.02.010 (2019).CAS 
Article 

Google Scholar 
Hommel, R. et al. Impact of interspecific competition and drought on the allocation of new assimilates in trees. Plant Biol. 18, 785–796. https://doi.org/10.1111/plb.12461 (2016).CAS 
Article 
PubMed 

Google Scholar  More

Ellison, W. T., Southall, B. L., Clark, C. W. & Frankel, A. S. A new context-based approach to assess marine mammal behavioral responses to anthropogenic sounds. Conserv. Biol. 26, 21–28 (2012).CAS 
PubMed 
Article 

Google Scholar 
Williams, R., Clark, C. W., Ponirakis, D. & Ashe, E. Acoustic quality of critical habitats for three threatened whale populations. Anim. Conserv. 17, 174–185 (2014).Article 

Google Scholar 
Halliday, W. D., Pine, M. K. & Insley, S. J. Underwater noise and arctic marine mammals: Review and policy recommendations. Environ. Rev. 28, 438–448 (2020).Article 

Google Scholar 
Kvadsheim, P. H. et al. Impact of Anthropogenic Noise on the Marine Environment: Status of Knowledge and Management (Springer, 2020).
Google Scholar 
Weilgart, L. S. & Whitehead, H. Distinctive vocalizations from mature male sperm whales (Physeter macrocephalus). Can. J. Zool. 66, 1931–1937 (1988).Article 

Google Scholar 
Simon, M., Stafford, K. M., Beedholm, K., Lee, C. M. & Madsen, P. T. Singing behavior of fin whales in the Davis Strait with implications for mating, migration and foraging. J. Acoust. Soc. Am. 128, 3200 (2010).ADS 
PubMed 
Article 

Google Scholar 
Alves, D., Amorim, M. C. P. & Fonseca, P. J. Assessing acoustic communication active space in the Lusitanian toadfish. J. Exp. Biol. 219, 1122–1129 (2016).PubMed 

Google Scholar 
Linnenschmidt, M., Teilmann, J., Akamatsu, T., Dietz, R. & Miller, L. A. Biosonar, dive, and foraging activity of satellite tracked harbor porpoises (Phocoena phocoena). Mar. Mamm. Sci. 29, E77–E97 (2013).Article 

Google Scholar 
Giorli, G. & Goetz, K. T. Foraging activity of sperm whales (Physeter macrocephalus) off the east coast of New Zealand. Sci. Rep. 9, 12182 (2019).ADS 
PubMed 
PubMed Central 
Article 
CAS 

Google Scholar 
Baumgartner, M. F. & Fratantoni, D. M. Diel periodicity in both sei whale vocalization rates and the vertical migration of their copepod prey observed from ocean gliders. Limnol. Oceanogr. 53, 2197–2209 (2008).ADS 
Article 

Google Scholar 
Urazghildiiev, I. R. & Van Parijs, S. M. Automatic grunt detector and recognizer for Atlantic cod (Gadus morhua). J. Acoust. Soc. Am. 139, 2532–2540 (2016).ADS 
PubMed 
Article 

Google Scholar 
Ladich, F. Ecology of sound communication in fishes. Fish Fish. 20, 552–563 (2019).Article 

Google Scholar 
Radford, C. A., Stanley, J. A., Simpson, S. D. & Jeffs, A. G. Juvenile coral reef fish use sound to locate habitats. Coral Reefs 30, 295–305 (2011).ADS 
Article 

Google Scholar 
Pierpoint, C. Harbour porpoise (Phocoena phocoena) foraging strategy at a high energy, near-shore site in south-west Wales, UK. J. Mar. Biol. Assoc. UK 88, 1167–1173 (2008).Article 

Google Scholar 
Pijanowski, B. C. et al. Soundscape ecology: The science of sound in the landscape. Bioscience 61, 203–216 (2011).Article 

Google Scholar 
Stanley, J. A., Radford, C. A. & Jeffs, A. G. Location, location, location: Finding a suitable home among the noise. Proc. R. Soc. B Biol. Sci. 279, 3622–3631 (2012).Article 

Google Scholar 
Buscaino, G. et al. Temporal patterns in the soundscape of the shallow waters of a Mediterranean marine protected area. Sci. Rep. 6, 1–13 (2016).Article 
CAS 

Google Scholar 
Gasc, A., Francomano, D., Dunning, J. B. & Pijanowski, B. C. Future directions for soundscape ecology: The importance of ornithological contributions. Auk 134, 215–228 (2017).Article 

Google Scholar 
Putland, R. L., Constantine, R. & Radford, C. A. Exploring spatial and temporal trends in the soundscape of an ecologically significant embayment. Sci. Rep. 7, 5713 (2017).ADS 
CAS 
PubMed 
PubMed Central 
Article 

Google Scholar 
Pieretti, N., LoMartire, M., Farina, A. & Danovaro, R. Marine soundscape as an additional biodiversity monitoring tool: A case study from the Adriatic Sea (Mediterranean Sea). Ecol. Indic. 83, 13–20 (2017).Article 

Google Scholar 
Gillespie, D., Palmer, L., Macaulay, J., Sparling, C. & Hastie, G. Passive acoustic methods for tracking the 3D movements of small cetaceans around marine structures. PLoS ONE 15, e0229058 (2020).CAS 
PubMed 
PubMed Central 
Article 

Google Scholar 
Van Parijs, S. et al. Management and research applications of real-time and archival passive acoustic sensors over varying temporal and spatial scales. Mar. Ecol. Prog. Ser. 395, 21–36 (2009).ADS 
Article 

Google Scholar 
Warren, V. E., McPherson, C., Giorli, G., Goetz, K. T. & Radford, C. A. Marine soundscape variation reveals insights into baleen whales and their environment: a case study in central New Zealand. R. Soc. Open Sci. 8, 201503 (2021).ADS 
PubMed 
PubMed Central 
Article 

Google Scholar 
Ahonen, H. et al. The underwater soundscape in western Fram Strait: Breeding ground of Spitsbergen’s endangered bowhead whales. Mar. Pollut. Bull. 123, 97–112 (2017).CAS 
PubMed 
Article 

Google Scholar 
Hildebrand, J. A. Anthropogenic and natural sources of ambient noise in the ocean. Mar. Ecol. Prog. Ser. 395, 5–20 (2009).ADS 
Article 

Google Scholar 
Ross, D. Ship sources of ambient noise. IEEE J. Ocean. Eng. 30, 257–261 (2005).ADS 
Article 

Google Scholar 
Popper, A. N. & Hawkins, A. The Effects of Noise on Aquatic Life Vol. 730 (Springer, 2012).Book 

Google Scholar 
Hubert, J. et al. Effects of broadband sound exposure on the interaction between foraging crab and shrimp: A field study. Environ. Pollut. 243, 1923–1929 (2018).CAS 
PubMed 
Article 

Google Scholar 
Weilgart, L. The Impact of Ocean Noise Pollution on Fish and Invertebrates (Springer, 2018).
Google Scholar 
Kvadsheim, H., Sivle, L. D., Hansen, R. R. & Karlsen, H. E. Effekter av Menneskeskapt støy på Havmiljø Rapport til Miljødirektoratet om Kunnskapsstatus FFI-RAPPORT. (2017).Parks, S. E., Johnson, M., Nowacek, D. & Tyack, P. L. Individual right whales call louder in increased environmental noise. Biol. Lett. 7, 33–35 (2011).PubMed 
Article 

Google Scholar 
Meh, F. et al. Humpback whales Megaptera novaeangliae alter calling behavior in response to natural sounds and vessel noise. Mar. Ecol. Prog. Ser. 607, 251–268 (2018).Article 

Google Scholar 
Leroy, E. C., Royer, J.-Y., Bonnel, J. & Samaran, F. Long-term and seasonal changes of large whale call frequency in the Southern Indian ocean. J. Geophys. Res. Ocean. 123, 8568–8580 (2018).ADS 
Article 

Google Scholar 
Parks, S. E., Clark, C. W. & Tyack, P. L. Short- and long-term changes in right whale calling behavior: The potential effects of noise on acoustic communication. J. Acoust. Soc. Am. 122, 3725–3731 (2007).ADS 
PubMed 
Article 

Google Scholar 
Clark, C. et al. Acoustic masking in marine ecosystems: intuitions, analysis, and implication. Mar. Ecol. Prog. Ser. 395, 201–222 (2009).ADS 
Article 

Google Scholar 
PAME. Underwater Noise in the Arctic: A State of Knowledge Report (PAME, 2019).
Google Scholar 
Beszczynska-Möller, A., Woodgate, R., Lee, C., Melling, H. & Karcher, M. A synthesis of exchanges through the main oceanic gateways to the Arctic Ocean. Oceanography 24, 82–99 (2011).Article 

Google Scholar 
Ramm, T. Hungry During Migration? Humpback Whale Movement from the Barents Sea to a Feeding Stopover in Northern Norway Revealed by Photo-ID Analysis. (MSc thesis. UiT The Arctic University of Norway, 2020).Broms, F. et al. Recent research on the migratory destinations of humpback whales (Megaptera novaeangliae) from a mid-winter feeding stop-over area in Northern Norway. in Recent research on the migratory destinations of humpback whales (Megaptera novaeangliae) from a mid-winter feeding stop-over area in Northern Norway (ed. Wenzel, F. W.) (European Cetacean Society Special Publication Series, 2015).Aniceto, A. S. et al. Arctic marine data collection using oceanic gliders: Providing ecological context to cetacean vocalizations. Front. Mar. Sci. 7, 547 (2020).Article 

Google Scholar 
Jourdain, E. & Vongraven, D. Humpback whale (Megaptera novaeangliae) and killer whale (Orcinus orca) feeding aggregations for foraging on herring (Clupea harengus) in Northern Norway. Mamm. Biol. 86, 27–32 (2017).Article 

Google Scholar 
Christiansen, J. S., Mecklenburg, C. W. & Karamushko, O. V. Arctic marine fishes and their fisheries in light of global change. Glob. Chang. Biol. 20, 352–359 (2014).ADS 
PubMed 
Article 

Google Scholar 
Rødland, E. S. & Bjørge, A. Residency and abundance of sperm whales (Physeter macrocephalus) in the Bleik Canyon, Norway. Mar. Biol. Res. 11, 974–982 (2015).Article 

Google Scholar 
Nøttestad, L. et al. Prey selection of offshore killer whales Orcinus orca in the Northeast Atlantic in late summer: spatial associations with mackerel. Mar. Ecol. Prog. Ser. 499, 275–283 (2014).ADS 
Article 

Google Scholar 
Bjørge, A., Aarefjord, H., Kaarstad, S., Kleivane, L. & Øien, N. Harbour porpoise (Phocoena phocoena) in Norwegian waters (Springer, 1991).
Google Scholar 
Gjøseter, H. et al. Fisken og Havet. https://doi.org/10.1111/maec.12075 (2010).Article 

Google Scholar 
ICES. ICES Report on Ocean Climate 2009 No.304. (2010).ICES. Report of the Working Group on Widely Distributed Stocks (WGWIDE). (2010).Rey, F. Phytoplankton: The grass of the sea. In The Norwegian Sea Ecosystem (ed. Skjoldal, H. R.) 97–136 (Academic Press, 2004).
Google Scholar 
Huse, G. et al. Effects of interactions between fish populations on ecosystem dynamics in the Norwegian Sea : Results of the INFERNO project. Mar. Biol. Res. 8, 415–419 (2012).Article 

Google Scholar 
Godø, O. R., Johnsen, S. & Torkelsen, T. The LoVe ocean observatory is in operation. Mar. Technol. Soc. J. 48, 24–30 (2014).Article 

Google Scholar 
Cooke, J. G. Balaenoptera physalus. The IUCN Red List of Threatened Species: e.T2478A50349982 (2018).Leonard, D. & Øien, N. Estimated abundances of Cetacean species in the Northeast Atlantic from Norwegian Shipboard Surveys Conducted in 2014–2018. NAMMCO Sci. Publ. 11, 4694 (2020).
Google Scholar 
Øygard, S. H. Simulations of Acoustic Transmission Loss of Fin Whale Calls Reaching the LoVe Ocean Observatory. (MSc thesis. University of Bergen, 2018).Steiner, L. et al. A link between male sperm whales, Physeter macrocephalus, of the Azores and Norway. J. Mar. Biol. Assoc. UK 92, 1751–1756 (2012).Article 

Google Scholar 
Olafsen, T., Winther, U., Olsen, Y. & Skjermo, J. Verdiskaping Basert på Produktive hav i 2050 1–76 (Springer, 2012).
Google Scholar 
Wenz, G. M. Acoustic ambient noise in the ocean: Spectra and sources. J. Acoust. Soc. Am. 34, 1936–1956 (1962).ADS 
Article 

Google Scholar 
Klinck, H. et al. Seasonal presence of cetaceans and ambient noise levels in polar waters of the North Atlantic. J. Acoust. Soc. Am. 132, 176–181 (2012).Article 

Google Scholar 
Burnham, R. E., Duffus, D. A. & Mouy, X. The presence of large whale species in Clayoquot Sound and its offshore waters. Cont. Shelf Res. 177, 15–23 (2019).ADS 
Article 

Google Scholar 
Romagosa, M. et al. Baleen whale acoustic presence and behaviour at a Mid-Atlantic migratory habitat, the Azores Archipelago. Sci. Rep. 10, 61489 (2020).
Google Scholar 
Tervo, O. Acoustic Behaviour of Bowhead Whales Balaena mysticetus in Disko Bay, Western Greenland. PhD thesis. (2011).Magnúsdóttir, E. E. & Lim, R. Subarctic singers: Humpback whale (Megaptera novaeangliae) song structure and progression from an Icelandic feeding ground during winter. PLoS ONE 14, e0210057 (2019).PubMed 
PubMed Central 
Article 
CAS 

Google Scholar 
Samaran, F. et al. Seasonal and geographic variation of southern blue whale subspecies in the Indian Ocean. PLoS ONE 8, e70 (2013).Article 

Google Scholar 
Norris, T. F., Dunleavy, K. J., Yack, T. M. & Ferguson, E. L. Estimation of minke whale abundance from an acoustic line transect survey of the Mariana Islands. Mar. Mammal Sci. 33, 574 (2017).Article 

Google Scholar 
Marques, T. A. et al. Estimating animal population density using passive acoustics. Biol. Rev. Camb. Philos. Soc. 88, 287–309 (2013).PubMed 
Article 

Google Scholar 
Dunlop, R. A. The effects of vessel noise on the communication network of humpback whales. R. Soc. Open Sci. 6, 190967 (2019).ADS 
PubMed 
PubMed Central 
Article 

Google Scholar 
Christensen, I., Haug, T. & Øien, N. A review of feeding and reproduction in large baleen whales (Mysticeti) and sperm whales Physeter macrocephalus in Norwegian and adjacent waters. ICES J. Mar. Sci. 49, 341–355 (1992).Article 

Google Scholar 
Aniceto, A. S. et al. Monitoring marine mammals using unmanned aerial vehicles: quantifying detection certainty. Ecosphere 9, e02122 (2018).Article 

Google Scholar 
Pedersen, G., Storheim, E., Sivle, L. D., Godø, O. R. & Ødegaard, L. A. Concurrent passive and active acoustic observations of high-latitude shallow foraging sperm whales (Physeter macrocephalus) and mesopelagic prey layer. J. Acoust. Soc. Am. 141, 1–10 (2017).Article 

Google Scholar 
Vogel, E. F. The influence of herring (Clupea harengus) biomass and distribution on killer whale (Orcinus orca) movements on the Norwegian shelf (UiT The Arctic University of Norway, 2020).
Google Scholar 
Williams, R. et al. Impacts of anthropogenic noise on marine life: Publication patterns, new discoveries, and future directions in research and management. Ocean Coast. Manag. 115, 17–24 (2015).Article 

Google Scholar 
Garibbo, S. et al. Low-frequency ocean acoustics: Measurements from the Lofoten-Vesterålen Ocean Observatory, Norway. (2020). https://doi.org/10.1121/2.0001324.Dekeling, R. P. A. et al. Monitoring Guidance for Underwater Noise in European Seas, Part I: Executive Summary (Springer, 2014).
Google Scholar 
Erbe, C. International regulation of underwater noise. Acoust. Aust. 41, 1–10 (2013).
Google Scholar 
Halliday, W. D., Insley, S. J., Hilliard, R. C., de Jong, T. & Pine, M. K. Potential impacts of shipping noise on marine mammals in the western Canadian Arctic. Mar. Pollut. Bull. 123, 73–82 (2017).CAS 
PubMed 
Article 

Google Scholar 
Halliday, W. D. et al. Underwater sound levels in the Canadian Arctic, 2014–2019. Mar. Pollut. Bull. 168, 112437 (2021).CAS 
PubMed 
Article 

Google Scholar 
Ødegaard, L., Pedersen, G. & Johnsen, E. Underwater Noise From Wind At the High North Love Ocean Observatory. UACE 2019 Conf. Proc. 359–366 (2019).Zhang, G., Forland, T. N., Johnsen, E., Pedersen, G. & Dong, H. Measurements of underwater noise radiated by commercial ships at a cabled ocean observatory. Mar. Pollut. Bull. 153, 110948 (2020).CAS 
PubMed 
Article 

Google Scholar 
Maystrenko, Y. P., Olesen, O., Gernigon, L. & Gradmann, S. Deep structure of the Lofoten–Vesterålen segment of the Mid-Norwegian continental margin and adjacent areas derived from 3-D density modeling. J. Geophys. Res. Solid Earth 122, 1402–1433 (2017).ADS 
Article 

Google Scholar 
Gillespie, D. et al. PAMGUARD: Semiautomated, open source software for real-time acoustic detection and localization of cetaceans. J. Acoust. Soc. Am. 125, 2547–2547 (2009).ADS 
Article 

Google Scholar 
Hollander, M. & Wolfe, D. A. Nonparametric Statistical Methods (Wiley, 1973).MATH 

Google Scholar 
Vogel, E. F. et al. Killer whale movements on the Norwegian shelf are associated with herring density. Mar. Ecol. Prog. Ser. 665, 217–231 (2021).ADS 
Article 

Google Scholar 
Garcia, H. A. et al. Temporal-spatial, spectral, and source level distributions of fin whale vocalizations in the Norwegian Sea observed with a coherent hydrophone array. ICES J. Mar. Sci. 76, 268–283 (2019).Article 

Google Scholar 
Davis, G. E. et al. Exploring movement patterns and changing distributions of baleen whales in the western North Atlantic using a decade of passive acoustic data. Glob. Chang. Biol. 26, 4812–4840 (2020).ADS 
PubMed 
PubMed Central 
Article 

Google Scholar 
Risch, D. et al. Minke whale acoustic behavior and multi-year seasonal and diel vocalization patterns in Massachusetts Bay, USA. Mar. Ecol. Prog. Ser. 489, 279–295 (2013).ADS 
Article 

Google Scholar 
Le Tixerant, M., Le Guyader, D., Gourmelon, F. & Queffelec, B. How can Automatic Identification System (AIS) data be used for maritime spatial planning?. Ocean Coast. Manag. 166, 18–30 (2018).Article 

Google Scholar 
Team, R. C. R: A Language and Environment for Statistical Computing. (R Foundation for Statistical Computing, 2020).Sumner, M. D. The Tag Location Problem. 133 (2011).Sumner, M. D., Wotherspoon, S. J. & Hindell, M. A. Bayesian estimation of animal movement from archival and satellite tags. PLoS ONE 4, e7324 (2009).ADS 
PubMed 
PubMed Central 
Article 
CAS 

Google Scholar 
Halliday, W. D. et al. The coastal Arctic marine soundscape near Ulukhaktok, Northwest Territories, Canada. Polar Biol. 43, 623–636 (2020).Article 

Google Scholar 
Ezzet, F. & Pinheiro, J. Linear, generalized linear, and nonlinear mixed effects models. Pharm. Sci. Quant. Pharmacol. 1, 103–135. https://doi.org/10.1002/9780470087978.ch4 (2006).Article 

Google Scholar 
Mazerolle, M. J. AICcmodavg: Model Selection and Multimodel Inference Based on (Q)AIC(c). (2020).Wickham, H. ggplot2: Elegant Graphics for Data Analysis (Springer, 2016).MATH 
Book 

Google Scholar 
Pante, E. & Simon-Bouhet, B. marmap: A package for importing, plotting and analyzing bathymetric and topographic data in R. PLoS ONE 8, e73051 (2013).ADS 
CAS 
PubMed 
PubMed Central 
Article 

Google Scholar 
Wessel, P. & Walter, H. F. S. A global self-consistent, hierarchical, high-resolution shoreline database. J. Geophys. Res. 101, 8741–8743 (1996).ADS 
Article 

Google Scholar 
Sueur, J., Aubin, T. & Simonis, C. Equipment review: Seewave, a free modular tool for sound analysis and synthesis. Bioacoustics 18, 213–226 (2008).Article 

Google Scholar 
The Mathworks Inc. MATLAB (R2019a). (MathWorks Inc., 2019).Merchant, N. D. et al. Measuring acoustic habitats. Methods Ecol. Evol. 6, 257–265 (2015).PubMed 
PubMed Central 
Article 

Google Scholar  More

Plant materialsAn ALS-inhibitor-resistant johnsongrass (resistant to nicosulfuron) obtained from the University of Nebraska-Lincoln (source credit: Dr. John Lindquist) was used as the pollen source (male parent), and the natural johnsongrass population present in the experimental field at the Texas A&M University Farm, Somerville (Burleson County), Texas (30° 32′ 15.4″ N 96° 25′ 49.2″ W) with no history of ALS-inhibitor resistance was used as the pollen recipient (female parent). Prior to the initiation of the field experiment, the susceptibility to nicosulfuron of the natural johnsongrass population was verified by spraying Accent Q at the labeled field rate of 63 g ai ha−1 [mixed with 0.25% v/v Crop Oil Concentrate (COC)] on 10 randomly selected 1 m2 johnsongrass patches across the field area at 15–30 cm tall seedling stage. For this purpose, a CO2 pressurized backpack sprayer was calibrated to deliver 140 L ha−1 of spray volume at an operating speed of 4.8 kmph. The natural johnsongrass population was determined to be completely susceptible to nicosulfuron.During spring 2018, the seeds of AR johnsongrass were planted in pots (14-cm diameter and 12-cm tall) filled with potting soil mixture (LC1 Potting Mix, Sungro Horticulture Inc., Agawam, MA, USA) at the Norman Borlaug Center for Southern Crop Improvement Greenhouse Research Facility at Texas A&M University. The environmental conditions were set at 26/22 °C day/night temperature regime and a 14-h photoperiod. In each pot, 5 seeds were planted and thinned to one healthy seedling at 1-leaf stage. Seedlings were supplied with sufficient water and nutrients (Miracle-Gro Water Soluble All Purpose Plant Food, Scotts Miracle-Gro Products Inc., 14111 Scottslawn Road, Marysville, OH 43041). A total of 50 seedlings were established in the greenhouse and were maintained until they reached about 10 cm tall, at which point they were sprayed with 2× the field rate of nicosulfuron (63 × 2 = 126 g ai ha−1) (mixed with 0.25% v/v COC). The herbicide was applied using a track-sprayer (Research Track Sprayer, DeVries, Hollandale, MN) fitted with a flat fan nozzle (TeeJet XR110015) that was calibrated to deliver a spray volume of 140 L ha−1 at 276 kPa pressure, and at an operating speed of 4.8 kmph. All treated seedlings that survived the herbicide application at 21 days after treatment (DAT) were then used as the pollen donor in the field gene flow experiment. All plant materials were handled in accordance with relevant guidelines and regulations. No permissions or licenses were required for collecting the johnsongrass samples from the experimental fields.Dose–response assaysThe degree of resistance/susceptibility to nicosulfuron of the AR and AS johnsongrass biotypes were determined using a classical dose–response experiment. The assay consisted of seven rates (0, 0.0625, 0.125, 0.25, 0.5, 1, and 2×) for the AS population and nine rates (0, 0.25, 0.5, 1, 2, 4, 8, 16, and 32×) for the AR population [1 × (field recommended rate) = 63 g ai ha−1 of Accent Q]. The experimental units were arranged in a completely randomized design with four replications. Seeds of AR and AS plants were planted in plastic trays (25 × 25 cm) filled with commercial potting-soil mix (LC1 Potting Mix, Sungro Horticulture Inc., Agawam, MA, USA) and maintained at 26/22 °C day/night cycle with a 14-h photoperiod in the greenhouse. Seedlings at 1–2 leaf stage were thinned to 20 seedlings per tray; four replications each of 20 seedlings per dose were considered. The seedlings were watered and fertilized as needed. The assay was conducted twice, thus a total of 160 seedlings were screened for each dose.The established seedlings were sprayed with the appropriate herbicide dose at the 10–15 cm tall seedling stage. The herbicide was applied using a track sprayer calibrated to deliver a spray volume of 140 L ha−1 at 4.8 kmph operating speed. Survival (%) and injury (%) were assessed at 28 DAT. Any plant that failed to grow out of the herbicide impact was considered dead. Plant injury was rated for each plot (i.e. on the 20 seedlings per rep) on a scale of 0–100%, where 0 indicates no visible impact compared to the nontreated control and 100 indicates complete death of all plants in the tray. Immediately after the visual ratings were completed, shoot biomass produced by the 20 plants from each tray was determined by harvesting all the tissues at the soil level and drying them in an oven at 60 °C for 72 h. Seedling mortality data were used for fitting dose–response curves that allowed for determining the lethal dose that caused 100% mortality of the susceptible biotype. This dose was used as a discriminant dose to distinguish between a hybrid (that confers resistance to nicosulfuron as a result of gene flow) and a selfed progeny (susceptible to nicosulfuron) in the field gene flow study.Field experimental location and set-upThe field experiment was conducted across two ENVs in 2018 (summer and fall) and one in 2019 (fall) at the Texas A&M University Farm, Somerville (Burleson County), Texas (30° 32′ 15.4″ N 96° 25′ 49.2″ W). The study site is characterized by silty clay loam soil with an average annual rainfall of 98.2 cm. The field experiment followed the Nelder-wheel design40, i.e. concentric donor-receptor design, a widely used method for gene flow studies, wherein the pollen-donors are surrounded by the pollen-receptors (Fig. 1). In this study, the AR plants (planted in the central block of the wheel) served as the pollen-donors, whereas the AS plants (present in the spokes) served as the pollen-receptors.Figure 1Aerial view of the experimental arrangement that was used to quantify pollen-mediated gene flow from ALS-inhibitor resistant (AR) to -susceptible (AS) johnsongrass at the Texas A&M University Research Farm near College Station, Texas. AR johnsongrass plants were transplanted in the pollen-donor block of 5 m diameter at the center of the field. The surrounding pollen-receptor area was divided into four cardinal (N, E, S, W) and four ordinal (NE, SE, SW, NW) directional blocks where naturally-existing AS johnsongrass plants were used as the pollen-recipients. AS panicles exhibiting flowering synchrony with AR plants were tagged at specific distances (5–50 m, at 5 m increments) along the eight directional arms. A tall-growing biomass sorghum border was established in the perimeter of the experimental site to prevent pollen inflow from outside areas.Full size imageThe center of the wheel was 5 m in diameter, and each spoke was 50 m long starting at the periphery of the central circular block. Thirty AR plants (pollen-donors) were transplanted in four concentric rings of 1, 5, 9, and 15 plants in the 5 m diameter central block, surrounded by the pollen-receptors (i.e. AS plants) (Fig. 1). The AR plants were contained within the central block during the 2 years of the field experiment by harvesting and removing all mature seeds and removing any expanding rhizomatous shoots. Further, field cultivation was completely avoided in the central block throughout the study period. Any newly emerging johnsongrass plants (seedling/rhizomatous) other than the transplanted AR plants in the central block were removed periodically by manual uprooting.The wheel consisted of eight spokes (i.e. directional blocks) arranged in four cardinal (N, E, S, W) and four ordinal (NE, SE, NW, SW) directions (Fig. 1). The plots to quantify gene flow frequency were arranged at 0 (border of the central block), 5, 10, 15, 20, 25, 30, 35, 40, 45, and 50 m distances from the central block in all eight directions (Fig. 1). Each plot measured 3 × 2 m and the area surrounding the plots was shredded prior to the booting stage with a Rhino® RC flail shredder (RHINOAG, INC., Gibson City, IL 60936).A tall-growing biomass sorghum border (6 m wide) was established surrounding the experimental area in all directions to prevent potential inflow of pollen from other Sorghum spp. in the nearby areas. Additionally, prevailing weather conditions, specifically wind direction, wind speed, relative humidity, and air temperature measured at 5-min intervals were obtained from a nearby weather station located within the Texas A&M research farm (http://afs102.tamu.edu/). The field did not require any specific agronomic management in terms of irrigation, fertilization, or pest management.Flowering synchrony, tagging, and seed harvestingAt peak flowering, when  > 50% of the plants in the AR block started anther dehiscence (i.e., pollen shedding), ten AS panicles (five random plants × 2 panicles per plant) that showed flowering synchrony with AR plants and displayed protruded, receptive stigma were tagged using colored ribbons at each distance and direction. At seed maturity, the tagged AS panicles were harvested separately for each distance and direction. Panicles were threshed, seeds were cleaned manually, and stored under room conditions until used in the herbicide resistance screening to facilitate after-ripening and dormancy release.Resistance screeningThe hybrid progeny produced on AS plants as a result of outcrossing with AR plants would be heterozygous for the allele harboring nicosulfuron resistance, and would exhibit survival upon exposure to the herbicide applied at the discriminant dose at which all wild type (AS) plants would die. The discriminant dose was determined using the dose–response study described above. Thus, the frequency of resistant plants in the progeny would represent outcrossing/gene flow (%).To effectively detect the levels of gene flow from AR to AS biotypes especially at low frequencies, the minimum sample size required for resistance screening was determined based on the following formula (Eq. 1)41:$${text{N }} = {text{ ln}}left( {{1} – P} right)/{text{ln}}left( {{1} – p} right),$$
(1)
where P is the probability of detecting a resistant progeny in the least frequent class and p is the probability of the least frequent class. Based on this formula, a minimum of 298 to as high as 916 plants were screened for each distance within each direction, allowing for a 1% detection level (p = 0.01) with a 95% (P = 0.95) confidence interval.Approximately one-year old progeny seeds harvested from the AS plants were scarified using a sandpaper for 15–20 s to release dormancy. The seeds for each distance within each direction were planted in four replicates of plastic trays (50 × 25 cm) filled with potting soil mixture (LC1 Potting Mix, Sungro Horticulture Inc., Agawam, MA, USA). The plants were raised at the Norman Borlaug Center for Southern Crop Improvement Greenhouse Research Facility at Texas A&M University. The greenhouse was maintained at 28/24 °C day/night temperature regime and a 14-h photoperiod. About 10–15 cm tall seedlings were sprayed with the discriminant dose of the ALS-inhibitor nicosulfuron (Accent Q, 95 g ai ha−1) using a spray chamber (Research Track Sprayer, DeVries, Hollandale, MN) fitted with a flat fan nozzle (TeeJet XR110015) that was calibrated to deliver a spray volume of 140 L ha−1 at 276 kPa pressure, operating at a speed of 4.8 kmph. At 28 DAT, percent seedling survival was determined based on the number of plants that survived the herbicide application out of the total number of plants sprayed. The number of plants in each tray was counted before spraying.Molecular confirmation of hybridsLeaf tissue samples were collected from thirty random surviving plants (putative resistant) in the herbicide resistance screening study for each of the three field ENVs, thus totaling 90 samples. Genomic DNA was extracted from 100 mg of young leaf tissue using the modified CTAB protocol42. The concentration of DNA was determined using a Nanodrop 1000 UV–Vis spectrophotometer (DeNovix DS-II spectrophotometer, DeNovix Inc., Wilmington, DE 19810, USA). DNA was then diluted to a concentration of 20 ng/µl for PCR assay. The nicosulfuron-resistant johnsongrass from Nebraska used in this study possessed the Trp574Leu mutation39. Hence, single nucleotide polymorphism (SNP) primers targeting a unique short-range haplotype of Inzen® sorghum (Val560Ile + Trp574Leu) were performed using the PCR Allele Competitive Extension (PACE) platform to confirm the resistant plants43. The SNP primers and the PACE genotyping master mix were obtained from Integrated DNA Technologies (IDT) Inc. (Coralville, IA) and 3CR Bioscience (Harlow CM20 2BU, United Kingdom), respectively. In addition to the two no-template controls (NTCs), two nicosulfuron-resistant johnsongrass, one wild-type johnsongrass, and one Inzen® sorghum were also used in the PCR.The PCR was performed according to the manufacturer’s protocol (Bio-Rad Laboratories, Inc., Hercules, CA), with denaturation for 15 min at 94 °C, followed by 10 cycles of denaturation at 94 °C for 20 s, annealing and extension at 65–57 °C for 60 s, 30 cycles of denaturation for 20 s at 94 °C, and annealing and extension for 60 s at 57 °C. Fluorescence of the reaction products were detected using a BMG PHERAStar plate reader that uses the FAM (fluorescein amidite) and HEX (hexachloro-fluorescein) fluorophores.Data analysisFor the dose–response assay, three-parameter sigmoidal curves (Eq. 2) were fit on the seedling mortality data for the AS and AR biotypes (with log of herbicide doses), using SigmaPlot version 14.0 (Systat Software Inc., San Jose, CA).$$y=b/[1+{exp}^{left(-(x-eright)/c)}],$$
(2)
where, y is the mortality (%), x is the herbicide dose (g ai ha−1), b is the slope around e, c is the lower limit (theoretical minimum for y normalized to 0%), and e = LD50 (inflection point, mid-point or estimated herbicide dose when y = 50%). Windrose plots that represented wind speed and frequency during the flowering window in each of the eight directions were created using a macro in Microsoft Excel. Progeny seedling survival (%) that represents gene flow (%) was determined using Eq. (3).$${text{PMGF }}left( {text{%}} right){ } = { }left( frac{X}{Y} right)_{{i,j{ }}} times { }100,$$
(3)
where, X is the number of plants that survived the herbicide application, Y is the total number of plants sprayed for ith distance in jth direction.To test whether gene flow frequencies varied among the directions, ANOVA was conducted using JMP PRO v.14 (SAS Institute, Cary, NC, USA), based on the average gene flow frequency values in each direction; ENVs were considered as replicates in this analysis. A non-linear regression analysis for gene flow rate, describing an exponential decay function (Eq. 4), was fit using SigmaPlot based on the gene flow frequencies observed at different distances pooled across the directions and ENVs.$$y=y0+left[atimes {exp}^{left(-btimes xright)}right],$$
(4)
where, y is the PMGF (%), x is the distance (m) from pollen source, y0 is the lower asymptote (theoretical minimum for y normalized to 0%), a is the inflection point, mid-point or estimated distance when y = 50%, and b is the slope around a.A Pearson correlation analysis was conducted to determine potential association between PMGF [overall PMGF, short-distance PMGF (5 m), and long-distance PMGF (50 m)] and the environmental parameters temperature, relative humidity, and dew point. Further, a correlation analysis was also conducted to understand the association between PMGF frequencies and specific wind parameters such as wind frequency, wind speed, and gust speed. The molecular data were analyzed using KlusterCaller 1.1 software (KBioscience). More

Frelich, L. E. Forest Dynamics and Disturbance Regimes, Study from Every Green and Deciduous Temperate Forest 287 (Cambridge University Press, 2002).Book 

Google Scholar 
Hadley, K. S. The role of disturbance, topography, and forest structure in the development of a montane forest landscape. J. Torrey Bot. Soc. 121(1), 47–61 (1994).Article 

Google Scholar 
Gracia, M., Montane, F., Pique, J. & Retana, J. Overstory structure and topographic gradients determining diversity and abundance of understory shrub species in temperate forests in central Pyrenees (NE Spain). For. Ecol. Manag. 242, 391–397 (2007).Article 

Google Scholar 
Scheller, R. M. & Mladenoff, D. J. Understory species patterns and diversity in old-growth and managed northern hardwood forests. Ecol. Appl. 12(5), 1329–1343 (2002).Article 

Google Scholar 
Sagheb-Talebi, K., Sajedi, T. & Pourhashemi, M. Forest of Iran, a Treasure from the Past, a Hope for the Future 145 (Springer, 2014).
Google Scholar 
Homami Totmaj, L., Alizadeh, K., Giahchi, P., Darvishi Khatooni, J. & Behling, H. Late Holocene Hyrcanian forest and environmental dynamics in the mid-elevated highland of the Alborz Mountains, northern Iran. Rev. Palaeobot. Palynol. 295, 104507 (2021).Article 

Google Scholar 
Vakili, M. et al. Resistance and resilience of Hyrcanian mixed forests under natural and anthropogenic disturbances. Front. For. Glob. Change 4, 98 (2021).Article 

Google Scholar 
Aguirre, O., Hui, G., von Gadow, K. & Jiménez, J. An analysis of spatial forest structure using neighbourhood-based variables. For. Ecol. Manag. 183(1–3), 137–145 (2003).Article 

Google Scholar 
Li, Y., Hui, G., Zhao, Z., Hu, Y. & Ye, S. Spatial structural characteristics of three hardwood species in Korean pine broad-leaved forest—Validating the bivariate distribution of structural parameters from the point of tree population. For. Ecol. Manag. 314, 17–25 (2014).Article 

Google Scholar 
Condit, R. et al. Spatial patterns in the distribution of tropical tree species. Science 288(5470), 1414–8 (2000).ADS 
CAS 
PubMed 
Article 

Google Scholar 
Lü, C. et al. Population structure and spatial patterns of Haloxylon ammodendron population along the northwestern edge of Junggar basin. J. Desert Res. 32, 380–387 (2012).ADS 

Google Scholar 
Fazlollahi Mohammadi, M., Jalali, S. G., Kooch, Y. & Theodose, T. A. The influence of landform on the understory plant community in a temperate Beech forest in northern Iran. Ecol. Res. 30, 385–394 (2015).Article 

Google Scholar 
Fazlollahi Mohammadi, M., Jalali, S. G., Kooch, Y. & Said-Pullicino, D. Slope Gradient and Shape Effects on Soil Profiles in the Northern Mountainous Forests of Iran. Euras. Soil Sci. 49(12), 1366–1374 (2016).ADS 
Article 

Google Scholar 
Fazlollahi Mohammadi, M., Jalali, S. G., Kooch, Y. & Said-Pullicino, D. The effect of landform on soil microbial activity and biomass in a Hyrcanian oriental beech stand. CATENA 149, 309–317 (2017).CAS 
Article 

Google Scholar 
Fazlollahi Mohammadi, M., Jalali, S. G., Kooch, Y. & Theodose, T. A. Tree species composition biodiversity and regeneration in response to catena shape and position in a mountain forest. Scand. J. For. Res. 32(1), 80–90 (2017).Article 

Google Scholar 
Harms, K. E., Condit, R., Hubbell, S. P. & Foster, R. B. Habitat association of tree and shrubs in a 50-ha neotropical forest plot. J. Ecol. 89, 947–959 (2001).Article 

Google Scholar 
Gunatilleke, C. V. S. et al. Species-habitat associations in a Sri Lank and ipterocap forest. J. Trop. Ecol. 22, 371–378 (2006).Article 

Google Scholar 
Rubino, D. L. & McCarthy, B. C. Evaluation of coarse woody debris and forest vegetation across topographic gradients in a southern Ohio forest. For. Ecol. Manag. 183, 221–238 (2003).Article 

Google Scholar 
Mohsennezhad, M., Shokri, M., Zal, H. & Jafarian, Z. The effects of soil properties and physiographic factors on plant communities distribution in Behrestagh Rangeland. Rangeland 4(2), 262–275 (2010).
Google Scholar 
Sefidi, K., Esfandiary Darabad, F. & Azaryan, M. Effect of topography on tree species composition and volume of coarse woody debris in an Oriental beech (Fagus orientalis Lipsky) old growth forests, northern Iran. IFOREST Biogeosci. For. 9, 658–665 (2016).Article 

Google Scholar 
Valipour, A. et al. Relationships between forest structure and tree’s dimensions with physiographical factors in Armardeh forests (Northern Zagros). Iran. J. For. Poplar Res. 21(1), 30–47 (2013).
Google Scholar 
Clark, P. J. & Evans, F. C. Distance to nearest neighbor as a measure of spatial relationships in populations. Ecology 35, 445–453 (1954).Article 

Google Scholar 
Naqinezhad, A. et al. The combined effects of climate and canopy cover changes on understorey plants of the Hyrcanian forest biodiversity hotspot in northern Iran. Glob. Change Biol. 28(3), 1103–1118 (2022).Article 

Google Scholar 
Pelissaria, A. L. et al. Geostatistical modeling applied to spatiotemporal dynamics of successional tree species groups in a natural Mixed Tropical Forest. Ecol. Indic. 78, 1–7 (2017).Article 

Google Scholar 
Pretzsch, H. & Zenner, E. K. Toward managing mixed-species stands: From parametrization to prescription. For. Ecosyst. 4, 19 (2017).Article 

Google Scholar 
Yousefi, S. et al. Spatio-temporal variation of throughfall in a hyrcanian plain forest stand in Northern Iran. J. Hydrol. Hydromech. 66(1), 97–106 (2018).Article 

Google Scholar 
Soil Survey Staff. Keys to Soil Taxonomy 12th edn. (USDA-Natural Resources Conservation Service, 2014).
Google Scholar 
Land Info, L. L. C. http://www.landinfo.com/country-iran.html. Accessed (2013).Beven, K. J. & Kirkby, M. J. A. Physically based, variable contributing area model of basin hydrology/Un modèle à base physique de zone d’appel variable de l’hydrologie du basin versant. Hydrol. Sci. J. 24(1), 43–69 (1979).Article 

Google Scholar 
Bourgeron, P. S. Spatial aspects of vegetation structure. In Ecosystems of the World 14A—Tropical Rain Forest Ecosystems, Structure and Function (ed. Golley, F. B.) 29–47 (Elsevier, 1983).
Google Scholar 
Moeur, M. Characterizing spatial patterns of trees using stem-mapped data. For. Sci. 39(4), 756–775 (1993).ADS 

Google Scholar 
Chokkalingam, U. & White, A. Structure and spatial patterns of trees in old-growth northern hardwood and mixed forests of northern Maine. Plant Ecol. 156(2), 139–160 (2001).Article 

Google Scholar 
Ferhat, K. A. R. A. Spatial patterns of longleaf pine (Pinus palustris Mill.): A case study. Euras. J. For. Sci. 9(3), 151–159 (2021).
Google Scholar 
Pommerening, A. Approaches to quantifying forest structures. Forestry 75(3), 305–324 (2002).Article 

Google Scholar 
Pommerening, A. & Särkkä, A. What mark variograms tell about spatial plant interactions. Ecol. Model. 251, 64–72 (2013).Article 

Google Scholar 
Goovaerts, P. Geostatistical tools for characterizing the spatial variability of microbiological and physico-chemical soil properties. Biol. Fertil. Soils. 27, 315–334 (1998).CAS 
Article 

Google Scholar 
Landim, P. M. B. & Sturaro, J. R. Krigagem indicativa aplicada à elaboração de mapas probabilísticos de riscos. Geomatematica, Texto didático, 6. DGA, IGCE, Universidade Estadual de São Paulo (UNESP), Rio Claro, São Paulo, Brazil. Available at: http://www.rc.unesp.br/igce/aplicada/textodi.html. Accessed 25/05/13 (2002).Deutsch, C. V. & Journel, A. G. GSLIB: Geostatistical Software Library and User’s Guide 119 (Oxford University Press, 1992).
Google Scholar 
Oliver, M. A. & Webster, R. Combining nested and linear sampling for determining the scale and form of spatial variation of regionalized variables. Geogr. Anal. 18, 227–242 (1986).Article 

Google Scholar 
Zhao, Z., Ashraf, M. I. & Meng, F. R. Model prediction of soil drainage classes over a large area using a limited number of field samples: A case study in the province of Nova Scotia, Canada. Can. J. Soil Sci. 93(1), 73–83 (2013).Article 

Google Scholar 
Brubaker, S. C., Jones, A. J., Lewis, D. T. & Frank, K. Soil properties associated with landscape position. Soil Sci. Soc. Am. J. 57, 235–239 (1993).ADS 
Article 

Google Scholar 
Bellingham, P. J. & Tanner, E. V. J. The influence of topography on tree growth, mortality, and recruitment in a tropical Montane Forest. Biotropica 32(3), 378–384 (2000).Article 

Google Scholar 
Luizao, R. C. C. et al. Variation of carbon and nitrogen cycling processes along a topographic gradient in a Central Amazonian forest. Glob. Change Biol. 10, 592–600 (2004).ADS 
Article 

Google Scholar 
Beaty, R. M. & Taylor, A. H. Spatial and temporal variation of fire regimes in a mixed conifer forest landscape, southern cascades, California, USA. J. Biogeogr. 28, 955–966 (2001).Article 

Google Scholar 
Castilho, C. V. et al. Variation in aboveground tree live biomass in a central Amazonian Forest: Effects of soil and topography. For. Ecol. Manag. 234, 85–96 (2006).Article 

Google Scholar 
Swanson, F. J., Kratz, T. K., Caine, N. & Woodmansee, R. G. Landform effects on eco-system patterns and processes. Biol. Sci. 38, 92–98 (1988).
Google Scholar 
Kooch, Y., Hosseini, S. M., Mohammadi, J. & Hojjati, S. M. Windthrow effects on biodiversity of natural forest ecosystem in local scale. Hum. Environ. 9(3), 65–72 (2011).
Google Scholar 
Köhl, M. & Gertner, G. Geostatistics in evaluating forest damage surveys: Considerations on methods for describing spatial distributions. For. Ecol. Manag. 95(2), 131–140 (1997).Article 

Google Scholar 
Habashi, H., Hosseini, S. M., Mohammadi, J. & Rahmani, R. Stand structure and spatial pattern of trees in mixed Hyrcanian beech forests of Iran. Iran. J. For. Poplar Res. 15(1), 64–55 (2007).
Google Scholar 
Von Oheimb, G., Westphal, C., Tempel, H. & Härdtle, W. Structural pattern of a near-natural beech forest (Fagus sylvatica) (Serrahn, North-east Germany). For. Ecol. Manag. 212, 253–263 (2005).Article 

Google Scholar 
Kunstler, G., Curt, T. & Lepart, J. Spatial pattern of beech (Fagus sylvatica L.) and oak (Quercus pubescens Mill.) seedlings in natural pine (Pinus sylvestris L.) woodlands. Eur. J. For. Res. 123(4), 331–337 (2004).Article 

Google Scholar 
Mosandl, R. & Kleinert, A. Development of oaks (Quercus petraea (Matt.) Liebl.) emerged from bird-dispersed seeds under old-growth pine (Pinus sylvestris L.) stands. For. Ecol. Manag. 106, 35–44 (1998).Article 

Google Scholar 
Hosseini, A., Jafari, M. R. & Askari, S. Investigation and recognition of ecological characteristics of sites of Persian oak and pistachio old trees in forests of Ilam province. Wood Sci. Technol. 26(4), 113–128 (2019).
Google Scholar 
Ghalandarayeshi, S., Nord-Larsen, T., Johannsen, V. K. & Larsen, J. B. Spatial patterns of tree species in Suserup Skov—A semi-natural forest in Denmark. For. Ecol. Manag. 406, 391–401 (2017).Article 

Google Scholar 
Petritan, I. C., Marzano, R., Petritan, A. M. & Lingua, E. Overstory succession in a mixed Quercus petraea-Fagus sylvatica old growth forest revealed through the spatial pattern of competition and mortality. For. Ecol. Manag. 326, 9–17 (2014).Article 

Google Scholar 
Watt, A. S. On the ecology of British Beech woods with special reference to their regeneration: Part II, sections II and III. The development and structure of beech communities on the Sussex downs. J. Ecol. 13, 27–73 (1925).Article 

Google Scholar 
Wiegand, T., Gunatilleke, S., Gunatilleke, N. & Okuda, T. Analyzing the spatial structure of a Sri Lankan tree species with multiple scales of clustering. Ecology 88, 3088–3102 (2007).PubMed 
Article 

Google Scholar 
Moradi, M., Marvie Mohadjer, M. R., Sefidi, K., Zobiri, M. & Omidi, A. Over matured beech trees (Fagus orientalis Lipsky.) component of close to nature forestry in northern Iran. J. For. Res. 23(2), 289–294 (2012).Article 

Google Scholar 
Lan, G. Y. et al. Spatial dispersion patterns of trees in a tropical rainforest in Xishuangbanna, southwest China. Ecol. Res. 24, 1117–1124 (2009).ADS 
Article 

Google Scholar 
Lan, G., Hu, Y., Cao, M. & Zhu, H. Topography related spatial distribution of dominant tree species in a tropical seasonal rain forest in China. For. Ecol. Manag. 262(8), 1507–1513 (2011).Article 

Google Scholar 
Menendez, I., Moreno, G., Fernando Gallardo Lancho, J. & Saavedra, J. Soil solution composition in forest soils of sierra de gata mountains, Central-Western Spain: Relationship with soil water content. Arid Land Res. Manag. 9(4), 495–502 (1995).
Google Scholar 
Kopecký, M., Macek, M. & Wild, J. Topographic Wetness Index calculation guidelines based on measured soil moisture and plant species composition. Sci. Total Environ. 757, 143785 (2021).ADS 
PubMed 
Article 
CAS 

Google Scholar 
Delfan Abazari, B., Sagheb-Talebi, K. & Namiranian, M. Development stages and dynamic of undisturbed Oriental beech (Fagus orientalis Lipsky) stands in Kelardasht region (Iran). Iran. J. For. Poplar Res. 12, 307–326 (2004) ((in Persian)).
Google Scholar 
Sagheb-Talebi K., Delfan Abazari B. & Namiranian M. Description of decay stage in a natural Oriental beech (Fagus orientalis Lipsky) forest in Iran, preliminary results. In Natural Forests in the Temperate Zone of Europe – Values and Utilization (eds. Commarmot, B. & Hamor, F.D.), Proceedings of conference in Mukachevo, Oct 13–17, 130–134 (2003).Christensen, M., Emborg, J. & Nielsen, A. B. The forest cycle of Suserup Skov: Revisited and revised. Ecol. Bull. 52, 33–42 (2007).
Google Scholar 
Dobrowolska, D. et al. A review of European ash (Fraxinus excelsior L.): Implications for silviculture. Forestry 84, 133–148 (2011).Article 

Google Scholar 
Akhani, H., Djamali, M., Ghorbanalizadeh, A. & Ramezani, E. Plant biodiversity of Hyrcanian relict forests, N Iran: An overview of the flora, vegetation, palaeoecology and conservation. Pak. J. Bot. 42(1), 231–258 (2010).
Google Scholar 
Pourmajidian, M. R. et al. Effect of shelterwood cutting method on forest regeneration and stand structure in a Hyrcanian forest ecosystem. J. For. Res. 21, 265–272 (2010).Article 

Google Scholar 
Szwagrzyk, J. & Szewczyk, J. Tree mortality and effects of release from competition in an old-growth Fagus-Abies-Picea stand. J. Veg. Sci. 12, 621–626 (2001).Article 

Google Scholar 
Janík, D. et al. Tree spatial patterns of Fagus sylvatica expansion over 37 years. For. Ecol. Manag. 375, 134–145 (2016).Article 

Google Scholar 
Amiri, M. Dynamics of Structural Characteristics of a Natural Unlogged Fagus orientalis Lipsky Stand during a 5-year’s Period in Shast-Kalate Forest, Gorgan, Iran, Ph.D. Dissertation, Gorgan University of Agricultural Sciences and Natural Resources, Gorgan (2013) (in Persian).Soofi, M. Effects of anthropogenic pressure on large mammal species in the Hyrcanian forest, Iran: Effects of poaching, logging and livestock grazing on large mammals (Doctoral dissertation, Dissertation, Göttingen, Georg-August Universität, 2018). More