Relationship between anuran larvae occurrence and aquatic environment in septentrional east Palearctic landscapes

The presence of amphibian larvae is restricted by both biotic and abiotic variables of the environment. Some of these variables are still undetermined in the septentrional eastern Palearctic where Rana amurensis, Strauchbufo raddei and Dryophytes japonicus are found in large numbers. In this study, we sampled 92 sites across Mongolia, Russia and the Democratic People’s Republic of Korea and measured biotic and abiotic water variables, as well as the height of flooded terrestrial and emergent aquatic vegetation at the breeding site. We determined that the presence of anuran larvae is generally, but not always, linked to pH and temperature. Rana amurensis was not significantly affected by any of the variables measured, while S. raddei was impacted by water conductivity and D. japonicus by pH, temperature and vegetation. Our results highlight a potential risk for these species due to the changes in aquatic variables in response to desertification.


Introduction
All amphibians above 40°N in East Asia rely on water bodies for the development of tadpoles. The characteristics of the aquatic environment is an important factor to the recruitment of populations and a limiting factor to the development of larvae (Altig and McDiarmid 1999;Loman and Lardner 2006). While biotic variables, such as aquatic vegetation, are important for predator avoidance (Kopp et al. 2006;Tavares-Junior et al. 2020), abiotic variables relate to water chemistry and specific physiological adaptations and metabolic processes (Hillman 2008). For instance, ion-exchange and respiratory circulation take place through a permeable skin allowing exchanges between internal and external fluids (Quaranta et al. 2009). Similarly to amphibians worldwide, species in Northeast Asia are also limited by water biotic and abiotic properties (Borzée et al. 2018;Kwon et al. 2021), including the three focal species of this study (Wang et al. 2009;Ko 2012;Huang et al. 2019).
Amphibian species, however, rely on different mechanisms to cope with aquatic environments and show varying degrees of tolerance (Hopkins and Brodie Jr. 2015), with variables, such as water chemistry (Viertel 1999), dissolved oxygen (Yu and Guo 2013), total dissolved solids (Browne et al. 2009) and water temperature (Indermaur et al. 2010), having different impacts on species. For instance, some species have a greater osmotic tolerance to brackish habitats (Wu et al. 2014), with larvae of species, such as Bufotes viridis, found in brackish water bodies with salinity up to 25 ppt (Katz 1973). Similarly, Fejervarya cancrivora larvae can be found in environments with salinity up to 35 ppt and can acclimatise up to 39 ppt (Hopkins and Brodie Jr. 2015). Vegetation is also important for some species as it provides shelter and food (Laurila 1998;Rieger et al. 2004) and, while most biotic variables are also critical, their impact is generally species specific (Altig and Mc-Diarmid 1999).
This study focuses on three anuran species from Northeast Asia. Rana amurensis, S. raddei and D. japonicus are locally abundant species on the eastern Palearctic landscapes of Mongolia, Russia and the Democratic People's Republic of Korea (DPR Korea hereafter). The water biotic and abiotic properties of the area are generally changing due to global change (Kelderman and Batima 2006) and aridification has become a threat (Shnitnikov 1975;Zolotokrylin et al. 2016). In addition, urbanisation, industrialisation and mining activities have contributed to the gradual pollution of waters (Hofmann et al. 2011). In this study, we correlate the presence of anuran tadpoles and the aquatic environment, so that a baseline is available to determine the ecological requirements of the focal species and it can be used in the future to examine the impact of regional and global changes.

Focal species
Our focal species were S. raddei, R. amurensis and D. japonicus. All three species range from northern China, through Mongolia and eastern Russia and into northern DPR Korea. Strauchbufo raddei is affected by genotoxicity depending on the intensity of exposure to metals such as Cu, Pb, Zn and Cd (Wang et al. 2009). Rana amurensis is extremely well adapted to cold and can withstand hypoxia for several months during hibernation (Shekhovtsov et al. 2020). The species with the largest range is D. japonicus and it generally breeds in freshwater (defined as < 1 ppt; Swenson and Baldwin 1965), although tadpoles have been found in 9.8 ppt salinity (Heo 2018).

Field survey
We conducted presence-absence surveys in 2017 and 2018 at 92 independent sites throughout Mongolia (n = 86), Russia (n = 5) and DPR Korea (n = 1). Despite the far distance between the site in DPR Korea and the other sites sampled in this study, we decided to maintain the datapoint in our dataset as the values for the site in DPR Korea were within the 50 percent percentile centred on the mean value for each of the variables. We conducted day-time encounter surveys focused on larvae of our focal species between 13 and 21 June 2017 and between 2 and 17 June 2018 (Fig. 1). All surveys were conducted right after the breeding season of the focal species, when tadpoles are expected to be developing in the water (Kuzmin et al. 2017).
The sites for sampling were selected if visible from the main road travelled for sampling (see Borzée et al. 2019 for details on the sites in Mongolia and Borzée et al. 2021 for details on the site in DPR Korea) and each site had to be at least 1 h drive (ca. 20 km) away from all other sites to ensure independence of sites and visible from the vehicle to be sampled. The surveyed locations were lakes, oxbow lakes, ponds, flooded unused quarries and stream pools. At each site, we conducted multiple observers' surveys along the edges and across the water body. We measured the water parameters and biotic and abiotic properties at each site (PC-STestr 35, OAKTON Instruments, USA): temperature (°C), pH, electrical conductivity (μs), total dissolved solids (TDS; mg/l) and salinity (ppt). The water samples were taken within 1 m from the bank of the water body, between 5 and 20 cm of depth. Proper standardisation, such as controlling for ambient temperature, was not possible in the field, but the data were standardised in that the variables were measured at the periphery of the water body and between 2 and 7 cm from the water surface of the water body. The occurrence of the species was determined by the presence of at least one larva. Tadpoles were captured through net sweep and identified in the field, based on morphological keys, such as general shape and teeth rows. All individuals were released at the site of capture. In addition, we measured the general average height of the vegetation at the site, averaging it along pre-defined values: 5 cm, 15 cm, 40 cm or 60 cm in height. The vegetation included partially flooded terrestrial vegetation such a reed-beds and smaller phragmites, as well as emergent aquatic vegetation, measured at a minimum of three representative points within the quarter of the water body located along the banks. The height of the vegetation was included as earlier observations suggested that it may relate to spawning preferences in the species.

Statistical analysis
We removed TDS from the analyses (n = 92) because the variable was correlated with all other variables (Pearson correlation: 0.21 ≤ R ≤ 0.92, p ≤ 0.041). We also removed salinity as the variable was correlated with pH, conductivity and vegetation (0.17 ≤ R ≤ 0.98, p ≤ 0.005). The remaining variables were not significantly correlated, did not display any outlier (except salinity, but the variable was not included in the analyses) and were generally normally distributed when evaluating the QQ plots.
To analyse the effects of the four remaining variables on the presence of anuran larvae for all species, we followed a two-step approach. First, to identify a general pattern referring to the occurrence of all the studied species, we first used the number of species found per site as the response variable. Thus, we relied on a Generalised Linear Model, with the presence of larvae as the dependent variable (binary encoded), weighted by the number of species (ranging from 1 to 3) found at the site and temperature, conductivity, vegetation and pH as covariates, not considering interactions between variables.
Next, to test for the impact of the variables on each species separately, we performed three independent multinomial logistic regressions with species' presence/ absence as the dependent variable and pH, temperature, conductivity and vegetation as covariates (only main effects). The values for vegetation followed the pre-determined averages.
In addition, although salinity is generally a factor limiting the presence of the species, we could not use the variables for the analyses. We conducted an additional binary logistic regression for each of the species with salinity as an independent variable and the encoded presence or absence of each of the three species as dependent variables under a main effect model in three different models. We then plotted the presence of the three species against salinity. All analyses were conducted in SPSS v.21.0 (SPSS Inc., Chicago, USA).
The three independent multinomial logistic regressions showed that for R. amurensis (n = 92), the overall model was significant (χ 2 = 21.81, df = 4, p < 0.001) and explained 43% of the variation (Nagelkerke R 2 = 0.43), although none of the variables was significantly explaining the presence of the species in relation to the four environmental variables (Table 1). In the case of S. raddei (n = 92), the model was significant (χ 2 = 16.55, df = 4, p = 0.002) and explained 22% of the variation (Nagelkerke R 2 = 0.22), and only conductivity significantly explained the presence of the species (Table 1). A higher conductivity was associated with the presence of the species (Table 2). Finally, for D. japonicus (n = 92) the model was significant (χ 2 = 102.01, df = 4, p < 0.001) and explained 89% of the variation (Nagelkerke R 2 = 0.89). The variables pH, temperature and vegetation significantly explained the presence of the species, albeit marginally for vegetation (Table 1). A higher pH, but lower temperature and higher vegetation were associated with the presence of the species (Table 2). Finally, the three binary logistic regressions were significant for S. raddei (Wald = 5.11, df = 1, p = 0.024), but not for R. amurensis (Wald = 1.24, df = 1, p = 0.265) or D. japonicus (Wald = 0.61, df = 1, p > 0.433). Rana amurensis tadpoles were found in water with a maximum salinity of 2.87 ppt, S. raddei with a maximum salinity of 2.06 ppt and D. japonicus in a maximum salinity of 0.59 ppt (Fig. 2). A maximum of 4.14 ppt was recorded in the absence of any of the focal species.

Discussion
Our results highlighted the significant impact of pH and temperature on the presence of anuran larvae in the shallow water bodies of the septentrional eastern Palearctic. These results are in line with the general literature on amphibian tadpoles and the species of this region seem to follow the same broad patterns (Altig and McDiarmid 1999;Kuzmin et al. 2017). The relationship with pH is likely related to the acidity of the water or other chemicals for which pH is a proxy, resulting in a set of limiting factors potentially leading to the death of individuals in case of relatively weak variations in values (Clark and Lazerte 1985;Horne and Dunson 1995). Sensitivity to pH and other variations in water parameters, such as TDS (here correlated with all other variables) and salinity (here correlated with pH), may be linked to the mineral salts present in large concentration in the region and leaching into the environment, as well as pollutants of human origins (Linhoff et al. 2011;Afonina and Tashlykova 2018). In addition, greater changes are expected in the near future as global changes, locally resulting in an increase in the rate of desertification and thus linked environmental variables, are becoming increasingly prominent threats (Shnitnikov 1975;Zolotokrylin et al. 2016). The sensitivity and tolerance to variations in water chemistry are, however, likely to be regionally mitigated by genetic variation and likely adaptation (Pierce and Wooten 1992). For instance, a difference in genetic expression between populations of Dryophytes cinereus makes coastal populations resistant to salinity levels that are deadly to the same species in other regions (Albecker et al. 2021). In addition, while the increase in specific variables is beneficial to the development of tadpoles, this relationship is generally limited to a specific level. For instance, while the significant correlation with an increase in temperature is likely related to the minimum requirement for growth at a relatively rapid pace due to the short duration of the warm season at high latitudes (Herreid and Kinney 1967), an increase in temperature above a certain threshold will have a negative impact on larvae when approaching the critical thermal temperature (Deutsch et al. 2008;Simon et al. 2015).
The fact that none of the variables tested was significant for R. amurensis may relate to the species ability to withstand hypoxia and the related increase in osmotic pressure (Shekhovtsov et al. 2020) and because the species may be able to partially detoxify internal fluids. The sensitivity of S. raddei to conductivity is likely linked to local variations in the landscape as minerals leaching chemicals into the environment are unlikely to be continuously present in the environment. In addition, local adaptations to the environment can also be mitigated by the synergetic action of environmental variables. For instance, we found adult individuals in water with salinity up to 4 ppt, but only at sites where we did not find tadpoles. The significant relationship of D. japonicus to pH, temperature and vegetation may be related to the points raised above for the two first variables. In addition, the occurrence of the species may be linked to the requirements of adults for vegetation cover when breeding in this landscape (Borzée et al. 2019). It is, however, important to note the risk of artefacts for D. japonicus due to the low sample size, especially in view of the species variable vegetation requirements in relation to humidity (Maslova 2001). Despite these restrictions due to environmental limitations, some populations are able to adapt to harsh environments by physiologically adjusting. For instance, R. amurensis and S. raddei were found in unusually high salinity, a pattern mirroring the tolerance of numerous species (Smith et al. 2007;Rios-López 2008;Hua and Pierce 2013;Hopkins and Brodie Jr. 2015).
In future studies, such patterns could be related to other variables and especially variation between populations, in relation to genetic and landscape connectivity. Intra-clade plasticity may be related to adaptation to changing environments and, thus, correlated with the potential to adapt to climate change. Predictions highlight the risk of salination of landscapes and clades with coastal populations able to cope with higher salinity may be able to adapt better in the coming decades.