Study species and environments
We conducted this study within Bangcheon stream in Jeju Island, South Korea (N33.4722, E126.5432). Jeju Island, a volcanic island, is largely covered by basalt layers that quickly drain surface water underground. As a consequence, in most streams, including our study site, water flows only during rain. In these ephemeral streams, oviposition sites were available in the form of pools of varying sizes. Along the stream, the pools were patchily distributed, with distances varying from less than a metre between each pool and its nearest neighbour (Fig 1A). Each pool remained isolated until rainfall occurred. Heavy rains usually swept most, if not all, organisms in the pools. Rains did not alter the location of each pool, although some pools temporarily merged into larger ones just after heavy rains. There were no fish in our study sites.
We chose a specific study site occupying approximately 730 m2 area. This area spanned a longitudinal distance of approximately 120 metres along the stream and contained 321 different pools when estimated just after rain when all available pools can be identified. In the study site, three frogs and one salamander species mainly occurred and reproduced during our study period from 8th March to 24th August, 2021: Rana uenoi, Hynobius quelpaertensis, Bombina orientalis, and Dryphytes japonicus. While Kaloula borealis also laid some eggs in our study site, their numbers were only a few and we did not include K. borealis in our analysis. The timing of reproduction largely overlapped between R. ueonoi and H. quelpaertensis, and between B. orientalis and D. japonicus (Fig. 1B). For R. uenoi, although reproduction started before our study period, we only analysed the eggs laid during our study period.
Field surveys to study oviposition preference in relation to abiotic and biotic factors
To monitor the oviposition choices of females, we visited all available pools within our study site between 0800 and 1300 on every non-rainy day of our study period. We refrained from visiting the stream on rainy days for the safety of the experimenters. Each pool was marked with a unique number on a rock surface adjacent to it, using non-toxic silicone.
We recorded several biotic and abiotic characteristics of each pool daily (see Table 1): (i) whether water was present or absent (Dry; a binary variable), (ii) the amount of leaf litter (Leaves; an ordinal variable), (iii) the quantity of stones as potential refuges (Stones; ordinal), (iv) the turbidity of the pool water (Turbidity; ordinal), (v) the presence of aquatic predators (Predators; binary), and (vi) the presence of mosquito larvae (MosqLarvae; binary). Variables concerning the presence of amphibian juveniles were coded for each species as follows: (vii) the presence of newly spawned eggs (NewEggs; binary), (viii) the presence of conspecific eggs that were not newly spawned (ConEggs; binary), (ix) the presence of conspecific tadpoles (ConTadpoles; binary), (x) the presence of heterospecific eggs that were not newly spawned (HeteroEggs; binary), and (xi) the presence of heterospecific tadpoles (HeteroTadpoles; binary).
Although we initially started by counting the number of con-/heterospecific eggs and tadpoles, we decided to use these variables in binary form due to the heterogeneity of errors among pools; pools with larger size, more turbid water, or a higher number of stones were more likely to have inaccuracies in their counts. Leaves, Stones, and Turbidity (all ordinal variables) were subjectively but consistently assessed and rated by a single experimenter (DO) throughout the survey. The criteria can be found in Table 1. Mosquito larvae have primarily been described as competitors, but also as predators or prey of tadpoles (Blaustein and Margalit 1994, 1996; Mokany and Shine 2003b). While we were unsure of the ecological relationship between our studied species and mosquito larvae, we nonetheless analysed correlational patterns between their presence and female oviposition choice. Invertebrates other than mosquito larvae, such as mayfly or Chironomidae larvae, were also occasionally present, but their occurrence was very low (less than 0.3% of pools) and thus were not included in our analysis.
Occasionally, a pool was transiently divided into multiple pools (under consecutive non-rainy days), or two pools merged together (this rarely happened just after rain). When a pool divided, we recorded the characteristics of each separate pool and used the averaged value for analysis, retaining the original pool ID for these transiently separated pools. When pools merged, we noted the characteristics of the merged pool and applied them to both original pools. During the breeding seasons of R. uenoi and H. quelpaertensis, from March to May, all 321 pools were surveyed. However, for the breeding seasons of B. orientalis and D. japonicus, spanning May to August, we limited the survey to 192 pools in the lower parts of our study site, where more oviposition activities were observed. This reduced scale was necessary to complete the surveys in a timely manner, as the increased oviposition activities resulted in longer times required to survey each pool.
Egg clutches from different females were clearly distinguishable for R. uenoi and H. quelpaertensis, in which eggs were either densely grouped together or found in a sac. However, distinguishing eggs from different mothers in the same pool posed a challenge for B. orientalis and D. japonicus. This difficulty arose because not only was the number of eggs highly variable, but also the boundaries between different egg clutches were not clearly defined. Therefore, instead of counting the number of egg clutches, we employed a binary variable to denote whether newly spawned eggs were found in each pool each day. This approach, however, resulted in an inability to detect instances where more than one female laid eggs in the same pool on the same day.
In addition to the pool characteristics that we surveyed daily, we also measured (i) pool volume (PoolVolume) and (ii) canopy coverage (CanopyCover). Unlike the aforementioned variables, these two characteristics were measured only once during our study period. While PoolVolume could change daily (decreasing gradually due to evaporation but increasing during rain), all pools were expected to be simultaneously affected by these changes. Therefore, we measured PoolVolume only once and used this as a size index for each pool. PoolVolume was estimated just after a rainy day when all pools were fully filled with water. The three-dimensional shape of each pool was highly irregular; thus, we crudely estimated PoolVolume by multiplying the surface area by the deepest depth of each pool. The surface area was measured by aerial photography using a drone (MAVIC PRO, DJI, Guangdong, China), with a ruler placed next to each pool for area calibration. CanopyCover was estimated using a mobile phone camera (Galaxy Note 8, Samsung, Seoul, South Korea). This measurement was taken on a summer day when foliage was dense. We positioned the camera just above the suspected centroid of the water surface, facing upwards, and took a photograph. The percentage of plant cover in the image was then calculated, using this percentage as an index of CanopyCover with ImageJ 1.80 (National Institute of Health, Maryland, USA). The camera’s angle of view was 77°, as specified by the manufacturer.
Abiotic factor analysis
We considered Dry, Leaves, Stones, Turbidity, PoolVolume, and CanopyCover as abiotic factors that could potentially affect female oviposition choice. All analyses were conducted separately for each species. Firstly, we calculated a representative index for each variable for each pool, generating a pool-level dataset. For Dry, we calculated the proportion of the days that each pool found to be dry during the breeding season of each species. For Leaves, Stones, and Turbidity, we first examined whether the intra-pool variation of each variable across the study period was smaller than the inter-pool variation. For this, we employed the intraclass correlation coefficient (ICC) test in the ‘irr’ package in R, using a one-way consistency model (Koo and Li 2016). Within-pool variability was lower than among-pool variability for all three variables (all P < 0.001, estimated ICC ranged from 0.39 to 0.77). Then, for each pool, we averaged each variable across the breeding season and used these averaged values to characterise each pool. PoolVolume was log-transformed to reduce distributional skewness. After all these processes, we had one representative value for each abiotic factor for each pool. We excluded transiently appearing pools in the abiotic factor analysis (i.e., when either a pool was divided into multiple pools during consecutive non-rainy days due to the exposure of bottom surfaces, or multiple pools were connected to each other just after rain) because averaging across the entire breeding season was impossible for those temporary pools.
As an index of each pool's oviposition frequency, we employed two different methods depending on the species. For R. uenoi and H. quelpaertensis, we used a binary variable to indicate whether eggs were found in each pool at least once during the breeding season. This approach was chosen because (i) the total number of egg clutches per pool was generally one or a few, and (ii) a small fraction of pools (14% for R. uenoi and 7% for H. quelpaertensis) were used for oviposition. Conversely, for B. orientalis and D. japonicus, a higher number of pools were used for oviposition (56% and 36%, respectively), with each pool often being used multiple times on different dates. Therefore, we counted the number of days new eggs were found in each pool throughout the breeding season and used this count variable for our analysis.
We fitted generalised linear models (GLZs) to examine the effect of abiotic variables on female oviposition choice. Firstly, we examined the presence of multicollinearity among predictors and found that Leaves and Turbidity were highly correlated (Pearson’s r > 0.78 for all species data). To resolve the multicollinearity issue, we removed the Turbidity variable from further analysis, leaving only five abiotic factors (Dry, Leaves, Stones, Pool Volume, and Canopy Cover). For R. uenoi and H. quelpaertensis, we fitted binomial GLZs using the binary response of whether each pool was used for oviposition at least once and the five abiotic factors as predictors. For B. orientalis and D. japonicus, we used the number of days new eggs were found during the breeding season for each pool as a response variable and the abiotic factors as predictors. Overdispersion was detected when we fitted Poisson GLZs, so we instead fitted negative binomial GLZs (Lindén and Mäntyniemi 2011). Visual inspections of predictors and our response variables suggested putative polynomial relationships for Leaves and PoolVolume for B. orientalis and D. japonicus; thus, we additionally included quadratic terms of Leaves and PoolVolume for these two species. We did not include interaction terms in our predictors because (i) we did not specifically predict interactive effects among the predictors, and (ii) with five main predictors, including unplanned interaction terms leads to the addition of 10 extra predictors (more when polynomial terms are simultaneously considered) for each species, which makes interpretations challenging in our multi-species datasets.
To estimate the parameters of variables and predict their effects, we applied multimodel inference and performed model averaging using the ‘MuMIn’ package (Grueber et al. 2011). We averaged the parameters of models that appeared within the 0.95 cumulative sum of AIC weights using conditional averaging.
Biotic factor analysis
We considered Predators, MosqLarvae, ConEggs, ConTadpoles, HeteroEggs, HeteroTadpoles as biotic factors that could potentially affect female oviposition site choice. For the main results, we did not distinguish between eggs and tadpoles, thus using the merged categories of ConEggs/Tadpoles (where either or both conspecific eggs or tadpoles were present) and HeteroEggs/Tadpoles as variables, in order to keep the main results from being overly complex. However, we nonetheless performed additional analyses, considering eggs and tadpoles separately, and included these results in the Supporting Information, which are also discussed.
Unlike abiotic factors, which did not vary much throughout the breeding season within each pool, biotic conditions changed frequently; thus, generating and analysing averaged pool-level data was inappropriate. The most ecologically relevant approach may be analysing each day separately, as the female's decision was made by comparing the conditions of currently available pools. However, this was statistically challenging in our data because only a limited number of pools were found with new eggs on most days, making GLZ analysis for each day impossible due to the excessive number of pools with no new eggs each day. Instead, we pooled the data across our survey period for each species, generated frequency tables, and examined whether the frequency of eggs laid in pools with specific biotic conditions deviated from what would be expected if females randomly laid eggs without considering the biotic condition. Pooling the data across all survey dates was necessary due to the low number of new egg clutches found in each individual survey (often fewer than three), which precluded us from estimating the expected frequencies for each specific date. In the biotic factor analysis, we were able to include temporary pools in the analysis, but their numbers were minor (less than 1% of the total pools).
Firstly, we excluded the pools that had never been used during the study period from the analysis. We considered that these pools were not used for oviposition, regardless of the biotic conditions, likely because the abiotic conditions were not favourable to frogs. For example, pools that were substantially smaller than other pools were never used for oviposition in all four species. Thus, we compared the biotic conditions among the pools that had been used for oviposition at least once for each species.
Then, for each biotic factor, we calculated the frequency of available pools and the pools containing a specific biotic factor. For instance, to investigate the potential influence of the presence of predators on female oviposition choice, we determined the frequency of available pools and the pools containing predators during each survey. Specifically, we used the biotic conditions from the day before newly spawned eggs were found because female oviposition mostly, if not always, occurred during the night. Therefore, the conditions experienced by females were not those in which the new eggs were found, but rather the conditions of the previous day. These frequencies were subsequently pooled across all survey dates. Utilising this information, we estimated the expected frequencies of pools with and without egg clutches if female frogs randomly laid their eggs irrespective of the presence of predators. Subsequently, we constructed a contingency table using the frequencies (expected and observed frequencies of pools with egg clutches in both predator-free and predator-present pools) and conducted a chi-square test. We performed the above procedures for all biotic factors.