Ecology

Environmental and competitive filtering is most important for future tree survivalWe find that individual functional traits of each tree were most important for individual tree survival (40–87%) for all study sites, followed by forest dynamics (16–28%) and functional diversity (10–26%) (Fig. 2). Nevertheless, importance proportions substantially varied in each study site. While individual functional traits were least important in the mixed mountainous forest under reference climate (no climate warming), individual functional traits showed highest influence in the mixed temperate forest under future climate change (Fig. 2).Figure 2Relative importance of functional diversity, forest dynamics and individual traits for individual tree survival under reference climate and future climate conditions (RCP 4.5). Panels correspond to alpine needle-leaved (A), mixed mountainous (B), mixed temperate (C) and temperate broad-leaved (D) forests. Left bars in each panel illustrate reference climate (no warming) and right bars future climate (RCP4.5). Colours indicate forest dynamics (grey), functional diversity (yellow) and individual functional traits (blue). Forest dynamics include the number of locally competing trees ( > 5 m in height) and local biomass as a proxy for the successional status.Full size imageTree survival depends on a mixture of environmental (e.g. climate) and natural competitive filtering, which excludes trees with trait combinations that underperform under local conditions16. Therefore, the high importance of individual functional traits across all study sites suggests a strong environmental and competitive filtering. Under future climate, the importance of individual functional traits generally increases or remains at high levels (Fig. 2). This shows that environmental and competitive filtering through functional traits are important processes to select best performing trees for the future, although being different for each forest type.Changing forest composition and trait shifts require large functional portfolio to secure forest resistanceUnder future climate, we observe trait shifts within plant functional types and strong changes of the forest composition (Figs. 3, 4). Especially in the alpine and mixed forests, the proportion of broad-leaved trees increases to at least ~ 70% towards the end of the twenty-first century (Fig. 3). The changing climate alters environmental and competitive filtering simultaneously, whereby broad-leaved trees become more productive, survive better and increasingly outcompete needle-leaved trees. For instance, in the two mixed forests survival probabilities of broad-leaved trees (high SLA) increase by about ~ 10%, whereas the survival of needle-leaved (low SLA) trees is reduced by 10% to 30% in a warmer climate (see Supplementary Figs. S7, S8, Panel A). Locally better adapted and competitive broad-leaved trees can replace needle-leaved trees if they die and secure the forest’s overall biomass in the future. Nevertheless, our simulated forests still contain significant amounts of needle-leaved trees in the year 2099 in the two coldest study areas (Fig. 3, red and blue lines). Therefore, mixed tree communities with high functional diversity, where broad- and needle-leaved trees coexist, contain a broad range of functional niches out of which the best suitable plant strategies emerge and result in better resistance to climate change.Figure 3Forest compositions and changes in the proportion of broad-leaved trees (summergreen and evergreen plant functional type combined) under climate change (RCP 4.5) from 2000 to 2099 for each study site. The fraction of broad-leaved trees, as simulated by LPJmL-FIT for each site, increases gradually in almost every forest type reaching at least about 70% by the end of the century. Pictures depict snapshots from visualization of model output in the years 2000 and 2099, respectively. For a full animation of all sites from 2000 to 2099 please see Supplementary Video 1.Full size imageFigure 4Trait distributions of specific leaf area in year 2000 and 2099, respectively, under future climate change (RCP 4.5). Arrows indicate trait shifts within plant functional types: BL-S Broad-leaved summergreen, BL-E Broad-leaved evergreen, T-NL Temperate needle-leaved, B-NL Boreal needle-leaved. For more detailed distributions see Supplementary Figs. S3 and S4.Full size imageSimultaneously, we observe strong trait shifts in SLA within plant functional types across all study sites under climate change (Fig. 4). In general, the community of broad-leaved trees shift to lower SLA, while boreal needle-leaved trees are strongly reduced or slowly replaced by their temperate equivalent with higher SLA (Supplementary Fig. S3A, dark blue colours). In contrast, wood density distributions remain relatively broad and do not shift strongly under climate change (Supplementary Fig. S4). Throughout the century, the increasingly warmer climate filters new trait combinations leading to changes in the community composition within and across PFTs (see Supplementary Discussion B). Those trait shifts emerge from changes in the composition within PFTs and newly establishing PFTs, and could be less drastic if trait adaptation of tree individuals was considered (see “Limitations and outlook” section). The points raised above show, that trait ranges within and between PFTs should be wide to cover potential future trait shifts that secure future forest resistance.All this suggests that functionally diverse forests are more resistant to future climate changes, due to their rich portfolio of traits. Broad trait distributions both within and between PFTs form the fundament for environmental and competitive filtering to select the most productive trees, securing the forest’s overall biomass under changing conditions. But can functional diversity further strengthen forest resistance beyond portfolio effects?Functional complementarity helps young trees to surviveWe find that, in addition to port-folio effects, functional diversity increases forest resistance by supporting the survival of young trees to changing climate conditions via trait complementarity. Our results indicate, that trees benefit from functional diversity if they grow in tree communities with high FR, high FDv and low FE (Supplementary Figs. S6–S9, Panel D–F in each figure). Here, functional traits lay highly separated (FDv and FE) and span a broad range in the functional trait space (FR), enabling functional complementarity. Under these conditions the survival of trees increases up to + 16.8% (± 1.6%) depending on the study site and climate (Table 1). This effect is highest in the alpine and mountainous forests (14–17%), whereas it is less prominent or has an opposite effect in the two temperate forests (− 7% to 6%). That suggests, that complementarity effects are stronger in cold-limited and mixed forests where a marked cold winter season fosters a co-existence between broadleaved and needle-leaved trees. Both PFTs are specialized in fixing carbon during different times of the annual cycle: Due to their leaf phenology, needle-leaved trees can already be productive when broad-leaved trees are still in progress of unfolding or shedding their leaves. On the other hand, broad-leaved trees are more productive than needle-leaved trees during warmer months. If coexistence is given, these phenological differences enable complementarity and reduce competition among PFTs. That overall increases tree survival, because trees can invest more carbon in their stems and defensive structures if competition is lower. Therefore, we argue that phenological complementary can enhance tree survival and thus forest resistance. An in-depth discussion of those mechanisms is further found in Supplementary Discussion A.Table 1 Additional survival probabilities for trees in each forest site under reference climate (central column) and future climate (RCP4.5, right column) in case FR and FDv are high, while FE is low.Full size tableSurprisingly, our results show that those complementary effects are much more important for small trees ( 10 m) in right panel, respectively. Functional diversity and forest dynamics are more important for small trees compared to large trees, whereas individual functional traits matter most for large trees. This pattern was found to be consistent across all sites (see Supplementary Table S6).Full size imageFunctionally diverse understoreys unlock the synergy of filtering and complementary effectsOur findings underline the role of functionally diverse trees in the understoreys for forest resistance. On the one hand, functional diversity supports the survival of understorey trees via functional trait complementarity. On the other hand, they form the fundament for competitive and environmental filtering. Only diverse tree communities have trait pools large enough to ensure that their tree portfolio holds trait combinations best suited for changing climate conditions. Therefore, we argue that functional diversity does not only support tree survival through complementarity, but is a prerequisite for filtering resistant trees in the first place.To profit constantly from functional diversity of the understory and ensure constant adaptation, a diverse age structure is a prerequisite. Depending on the forest type, trees are distributed in a broad range of different height and age classes in our study (Supplementary Fig. S5, Supplementary Video 1). This multi-aged structure is preserved under climate change (Supplementary Fig. S5) and allows gradual changes through constant environmental and competitive filtering in the future.In this study, we simulate forests without any human interference or management. Our results are therefore to be interpreted in the context of environmental and competitive filtering as observed in natural forests. Most managed forests lack this natural filtering effect as they are less dense and diverse in their age-structure. Functionally diverse trees in the understorey could provide the fundament for climate adapted multi-aged forests, as they constantly form new better adapted tree generations with natural competition and succession allowed. Therefore, we fully underline the importance of functionally diverse understorey trees and natural competition as the fundament for future forest resistance.Management implicationsThe results of this study highlight the importance of functionally diverse understorey trees. However, browsing by game might damage new tree saplings and limit tree diversity in the understorey. In addition, invasive species like Prunus serotina or herbaceous competition might hinder forest succession and the establishment of woody native species in European forests20,21. Therefore, regulating game, limiting the spread of invasive species and controlling herbaceous vegetation should be considered in future management practices where tree diversity in the understorey seeks to be increased or maintained. Moreover, insufficient dispersal of functionally different tree species might limit the establishment of functionally diverse trees in the understorey. Future forest management may consider to artificially plant functionally different tree species if dispersal from surrounding forests cannot be guaranteed. On the other hand, forests that already contain functionally diverse trees in the understorey should be preserved.In this study, a clear trend from needle-leaved to broad-leaved trees is captured at all sites, whereas within broad-leaved PFTs a shift to lower specific leaf area and higher leaf longevities indicates that future forests might especially benefit from longer vegetation periods (earlier leaf onset, later senescence). Therefore, forests containing broad-leaved tree individuals with high phenological plasticity could be more resistant. The broad simulated wood density ranges, which persist under climate change, imply beneficial effects for forest communities entailing a range of different growing strategies, i.e. early to late successional species. Therefore, we argue that forest fragmentation should be reduced or reversed to foster some natural dispersal of early and late successional species.This study intended to explore the potentials of functionally diverse forests as a possibility to stabilize forests under climate change over a large climatic gradient. The model used in this study operates on the more general level of functional traits and their diversity rather than on species level (see “Methods” section and Supplementary Methods A). Consequently, management implications regarding suitable specific tree species are beyond the scope of this study. However, we think that our results will stimulate the discussion on the importance of functional tree traits and their diversity for species selection.Limitations and outlookThis study focussed on identifying the importance of functional diversity for future tree survival to advance our understanding on the role of biodiversity for future forest resistance using the flexible-trait Dynamic Global Vegetation Model LPJmL-FIT. The general approach of LPJmL-FIT is to simulate biogeographic dynamics purely based on environmental and competitive filtering (see “Methods” section, Supplementary Methods A). Due to missing processes in the model and the ambiguity of former human influence, drawing site-specific implications on future forest dynamics must be taken with caution (see Supplementary Discussion D).Moreover, processes not yet captured in the LPJmL-FIT model, might play a role and could lower forest resistance in the future, which is why we recommend relying on a trait space as broad as possible.Including more climate scenarios would widen the envelope of possible future pathways by considering climate model uncertainties. Insect outbreaks and pathogens might put pressure to the already drought- and temperature-stressed trees and heavily accelerate mortality especially of needle-leaved trees, although functionally diverse forests are less vulnerable to bark beetle outbreaks22,23. Multi-layered forests showed higher growth resilience to structural disturbances such as wind-throw24, likely enhancing the importance of individual tree height and reduce the survival probability of large trees25. Belowground competition and trait plasticity could favour complementarity effects further. Variable rooting strategies could further reduce competition for soil water and thereby increase individual drought resistance of trees26. Trait plasticity can contribute to tree survival by widening niche, further increasing complementarity effects. However, trait plasticity remains one of the most challenging objectives in vegetation modelling as observational data and modelling approaches are scarce27, leaving it open how far trait relationships would hold under climate change. Considering more functional traits in our analysis might increase the overall predictive power of the random forest models. Even though the explained variance increased with the number of analysed traits explaining ecosystem properties in long-term grassland experiments, such improvement is limited as abiotic factors and their interactions with plant traits might be more important for prediction28. We conclude that simulating future forest dynamics dominated by environmental and natural competitive filtering requires to integrate both, abiotic and biotic drivers on forest dynamics. Machine learning techniques are increasingly used in forest ecological research—but mainly applied in the processing of field and forest inventory data29,30. Machine learning can help to understand the complexity of interactions and provide deeper insights into the underlying ecological process in a modelling study as we have shown here using LPJmL-FIT simulation results. Random forest analyses are suitable for a variety of data and applications because they are relatively robust to different data structures. Importance analysis can help to identify the role of underlying processes in complex models and to visualize their changes in a simple way. In doing so, model development is advanced by making use of large data sets, opening the door to further theory building and deeper understanding of plant trait ecology. More

UPRLIMET is our response to a need for a consistent method for predicting the upper extent of trout in all streams across land ownerships within our region. By developing and implementing the model using LiDAR-derived flowline hydrography, we offer a standardized, spatially explicit, spatially contiguous (where LiDAR hydrography is available), and high-quality fish-distribution layer based on the probability of fish presence. UPRLIMET maps both the probability of trout and the upper limit of trout across landscapes, ownerships, and jurisdictions, and better captures the upper extent of fish in headwater reaches relative to previous approaches allowing for a cross-boundary distribution map on which decision-makers and managers can base policies and regulations.This work provides a transferable prediction modeling framework for systematically and comprehensively estimating the upper distribution limit of fish, which could be calibrated and implemented in watersheds and for fish species around the globe. Although the dependency on LiDAR-derived data here may be seen as a limitation to broader implementation of this method, the method is scalable to any resolution, and LiDAR is becoming increasingly ubiquitous in the United States through the U.S. Geological Survey 3D Elevation Program, which is funding LiDAR acquisitions across the United States. Furthermore, LiDAR data is available globally via data from GEDI and ICESAT-2 satellites that offer coarser resolution (~ 25-m) data that are still superior to either ASTER or SRTM derived-DEMs26, 27.Minimizing prediction errors for the upper limit of trout is important to decision support and management planning because it ensures that forest-harvest regulations and management prescriptions are aligned. It is important to note that the prediction error estimates from this study are derived from the NSpCV process, except for models using 20% slope thresholds or unaltered parameterization of Fransen’s model13, because it is likely that the NSpCV estimates are conservative. They tended to overestimate error, as evidenced by the fact that the Refit model (i.e. Fransen’s optimal model13 refit to our data) exhibited a larger MAE than the unchanged optimal Fransen model13. This unexpected result was likely due to applying the NSpCV routine on the Refit model, resulting in the use of many intermediate models to characterize predictive performance using randomized subsets of independent training and test data. In contrast, the optimal Fransen model13 was developed independently using the data in this study and thus error could be evaluated directly without subsampling imposed by NSpCV.The relatively low error for the two-stage model that becomes UPRLIMET suggests that it more accurately characterizes the upper limit of fish than all other models considered in this study, including the Fransen model13, which has been used for estimating upper limit of fish regionally. Although some of the models exhibited relatively small differences in error relative to the model that became UPRLIMET, small differences in predicted upper limit locations when considered in aggregate across multiple watersheds can potentially alter management decisions and expected outcomes. Differences in predictive performance and error between UPRLIMET and the optimal Fransen model13 are likely attributed to high-accuracy hydrography and hydro-topographic data (as LiDAR-derived DEMs were not available in western Washington in 2006), which allowed a finer-scale of analysis (i.e., 5-m vs 10-m reaches). Additionally, the fact that UPRLIMET was fit to data solely from western Oregon likely offers predictive performance gains when applied to western Oregon when compared to the Fransen model13 that was fit to data from western Washington.Quantifying the predicted accuracy associated with applying UPRLIMET to western Washington will require new data and is outside the intended scope of this study. However, we think it is reasonable to infer findings from UPRLIMET across regions with similar climatic and hydro-topographic conditions including northwestern California, western Oregon, western Washington, and southwestern British Columbia, especially given the broad availability of LiDAR-derived DEMs. This conclusion is supported the fact that both the Fransen13 and Refit models produced similar logistic regression coefficients (Data S5) and similar Matthews Correlation Coefficients (Data S6), suggesting that feature space of the two models is similar. This evidence is further corroborated by the high degree of overlap observed among the distributions of each of the four predictor variables for both western Oregon and Washington. We acknowledge that UPRLIMET does not contain identical predictor variables to Fransen’s model13 but maintain that they are similar enough in purpose that it is reasonable to assume that the feature space similarities are retained.When we undertook this study, we hypothesized that a prediction model based on RF would offer superior predictive performance over those based on LR, given the availability of 67 predictor variables and RF’s demonstrated superior predictive performance in ecological applications23,24,25. However, our results suggests no improvement is offered by including more than four of the 67 environmental predictors examined, and that no clear advantage is offered by employing the more complex RF model, as evidenced by the top three of the top five prediction models being four-variable LR model algorithms (Fig. 3; Data S3.) The general importance of these variables to so many models is likely due to the strong linear relationships in the response of fish or no fish in logit space given the slopes of the curves in the partial dependence profiles (Fig. 4). This finding is congruent with the fundamental premise of LR, which is to explain and predict a response with a functional relationship, whereas RF deliberately focuses only on maximizing prediction accuracy with many decision trees28. Additional advantages to prediction models based on LR include the following: relatively better extrapolation performance over RF29, the simplicity of transferring a LR model to another processing platform using the model coefficients (versus the black box of RF decisions), and the immensely reduced computational processing times associated with LR model fitting and prediction. These advantages are especially key to this work, where there may be a desire to implement the model on other landscapes without the requisite expertise in doing so using the R software30. However, there are tradeoffs, as LR is more sensitive to the influence of outliers and multi-collinearity among variables, and overfitting is an increasing concern as the number of predictor variables increase, whereas RF tends to be robust to these concerns, but is more likely to produce a high-variance, low-bias prediction model31.Although there is no single, general explanation for distribution limits of species32, the intersection of stream size, slope, and elevation together locate the upper limit of fish. Stream size corresponds to major ecosystem changes along a stream continuum including for energy sources, ecosystem metabolism, habitat characteristics, and biodiversity33, as well as the upper distribution limit of fish, as shown here. As expected, stream size accounts for the top two variables in the model suggesting that it is the major driver of the upper distribution limit of fish with the probability of trout increasing with increasing upstream stream length and upstream drainage area. Our finding proposes that downstream stream reaches are more likely to have fish. Although the underlying mechanisms have multiple influences, factors related to increasing stream size, such as increasing habitat size, habitat complexity, stability, or temperature variability34 have been shown to be important. Similarly, stream size is the most sensitive factor in intrinsic potential models for Chinook Salmon (O. tshawytscha35). Slope, the next variable of importance influencing the upper extent of fish, exerts control on physical habitats in streams, including channel morphology, hydraulics, sediment transport, substrate, and habitat36. Steep slopes drastically prevent trout from reaching areas above waterfalls or impassable chutes of over 25% slope, but trout can be found in streams channels without barriers at slopes as high as 28%7, 14, 37. Other fishes, such as Coho Salmon (O. kisutch) and steelhead (O. mykiss) are generally not found above 12% slope38. Interestingly, survival of fishes that make it upstream or are introduced above barriers may be facilitated by a geomorphic setting that is less prone to debris flows and other episodic sediment fluxes and has a greater resilience to flooding resulting from wider valley and greater floodplain connectivity39. Elevation or vertical topographic position may indirectly integrate broad influences of other landscape-scale or climate factors or also indirectly capture stream size, influencing the likelihood of fish presence. Frequently, species richness increases at lower elevations40, and we suggest that elevation also contributes to species distribution limits, as is the case for the Endangered Species Act listed Bull Trout (Salvelinus confluentus)41. The multiple factors associated with elevation correspond to the relationship found for stream size that smaller streams are less likely to have fish. Ultimately, the intersection of stream size, slope, and elevation guide us to finding the upper extent of fish in streams.Physical influences have been proposed to be more limiting to fish distributions upstream, such as near the upper extent of fish, whereas biological factors are probably more important downstream33. Although 67 environmental predictor variables representing geologic, soil, climatic, and hydro-topographic conditions at local and patch scales are evaluated (Data S1), only the hydro-topographic variables of stream size, slope, and elevation are important to predicting the upper limit of fish in UPRLIMET. In fact, the top 9 models (Fig. 3; Data S3) relied on just four to five hydro-topographic variables, most of which were patch-scale variables or elevation at 1000 m, all of which incorporate a broader extent of influence. This suggests that local scale variables that contribute to fish limits, including slope or riparian influences may need to be further explored. In addition, some of the remaining 63 variables present in UPRLIMET, such as precipitation and air temperature, are important drivers of within-network trout distributions and contribute to their connectivity. Some of these predictor variables appear in the 10th ranked 26-variable RF-O-SR1 model (Data S2; Data S4; Data S8), but the influence appears to be dubious for isolating the upper limit and explaining variation in fish occurrence because MAE of upper limit was substantially higher than the 9 models with lower MAEs (Fig. 3; Data S3), and the lower MCC of the associated RF-O sub-model (Data S6). It is likely that other combinations of the 67 predictor variables, including precipitation, may be more important when this model development and evaluation framework is applied elsewhere, especially if those areas contain fishes or are places that are vulnerable to changing water temperatures and streamflow regimes. In addition, biological factors may be a concern in other watersheds, including invasive species and fish stocking which can limit the longitudinal distribution and the upstream extent of fishes.Given the large geographic extent of this study, we expected other variables such as precipitation to be more important drivers, however due to a combination of a wet water year, a lack of precipitation gradient in the study area, coarse grain data, and location of fish in streams this was not the case. For example, 2017 was a wetter than normal water year53, and it may be that the gradient of precipitation variation in western Oregon was not strong enough to explain the variation in the spatial distribution of trout occurrence. All climate data, including the precipitation data were sourced from relatively coarse-scale (800 m) PRISM data. The inability to adequately downscale precipitation to characterize how precipitation truly varies within and between patches, especially along elevational gradients, likely confounded how the model interprets the influence of precipitation. Trout occurrence was on perennial streams, which is likely far enough downstream of locations where variation in precipitation was the dominant influence on streamflow permanence and consequently would not have been a factor.Stream network structure plays a key role in the upper limits of fish. Upper limits for fish can occur at either lateral or terminal points13 and when mapping these points, differences were seen for UPRLIMET relative to other datasets. Lateral limits end in the tributary stream just above where it connects with a mainstem stream. Terminal limits include both mid-stream terminal limits where fish drop out in the middle of a stream channel owing to a soft (i.e., transient barrier or puttering out) or hard (i.e., waterfall) edge, and confluence terminal limits where the upper limit of fish ends at the confluence. For example, when closely examining the 14 watersheds where we have overlapping information across various datasets and models, UPRLIMET and the Fransen optimal model13 exhibit substantial agreement in their lateral limits. However, the largest differences are in their terminal ends, especially terminal mid-stream limits, probably owing to hydro-topographic changes that contribute to fish occurrence at confluences, which are more pronounced than mid-stream. Accordingly, the logic in the stopping rule is likely important in identifying specific upper extent of fish distributions in reaches that end mid-stream.Differences among databases for the upper distribution limits of fish come from both the upper limit points and depiction of fish-bearing reaches, underscoring the importance of having a shared map with common coverage of the fish extent across landscapes and ownerships. Differences among mapped distributions can result from source information, relating to whether it is modeled or occurrence data. Models, such as UPRLIMET, can be applied across a broad extent based on model parameters and training data, thereby offering broad coverage for distributions (and quantifiable error) across the landscape, ownerships, and jurisdictions. However, models are limited by accuracy and fit. As such, they can incorrectly predict distributions in some areas, especially if there are prediction features not yet trained with the model data where prediction would require extrapolation of the model. This makes both the training dataset and modeled extent important considerations, as models are only as good as the data used to develop them. Updating UPRLIMET with new data as it becomes available will help to expand the prediction domain, improve accuracy, and allow the model to do more interpolation than extrapolation.Distributions based on occurrence information depend heavily on data availability, data quality, and access. Differences in data availability can lead to inconsistent coverage across landscapes and ownerships, with high coverage in some watersheds and low to no coverage in others. Inconsistent coverage can lead to errors that are difficult to quantify across landscapes, ownerships, and survey crews. Occurrence information also depends on the ability to survey watersheds and gain access across ownership types, including on private lands that do not have the same assurances of access as public lands, resulting in information asymmetry42, 43. Data quality also depends on the spatial accuracy of the points of uppermost fish, which are a function of GPS quality and error, and can drastically change the modeled results, as these points are used in the training dataset. Differences among mapped distribution limits also result from differences in field protocols on designating last fish. For example, some crews note fish distribution limits where they visually see the last fish, whereas others note it upstream of where they saw last fish, based on habitat features that would limit fish. With the advent of LiDAR-derived DEMs and associated LiDAR-derived stream hydrography, like those available in much of western Oregon, have revealed additional flowlines in watersheds compared to previous topographic maps, which adds more potential tributaries to survey for fish-distribution assessments. When these new previously unmapped tributaries are paired with a model, such as UPRLIMET, a common information set is available across landowners, managers, and agencies for the upper extent of fish. This helps policymakers determine where to apply regulations that support fisheries and forest management, based on the upper fish limit.Next steps for applying and expanding the model include addressing current data gaps. More information and observations about the upper distribution limits of fish beyond western Oregon would be needed to properly expand the spatial scope of the model. The upper extent of fish is at the detection limit of many current technologies, including global nativation satellite system (GNSS), geographic information systems (GIS), and LiDAR, especially in forested landscapes. Better precision of GNSS coordinates from observations would help greatly. From an ecological perspective, we could focus on fish distribution limits that vary seasonally or interannually to better understand which stream features and hydrologic parameters influence those endpoints. We also need information related to locations of barriers, including culverts, waterfalls, and knickpoints to understand their influence on contemporary distributions. Incorporating variables representing riparian conditions as well as leveraging higher-resolution DEMs ( More

Time-course observations on 43 flowers of Causonis japonica revealed changes in flower disc colour and sexual expression (Table 1). Temporal changes in floral features showed no difference between diploid (19 flowers) and triploid (24 flowers) individuals. For example, flowering onset times did not differ substantially between ploidy level (diploid: from 07:07 to 13:27, triploid: from 06:58 to 14:49). However, the flowering duration varied significantly from flower to flower, ranging from a minimum of one day to a maximum of six days. Regardless of the ploidy level, all flowers with damaged styles (14 flowers) exhibited brown stigmas after the male phase, then ceased floral development prior to the female phase.Table 1 Characteristics of 43 flowers of Causonis japonica.Full size tableFigure 1 shows the typical time-course changes of C. japonica flower features (flower ID 35 in Table 1) according to the RGB values (representing the activities of nectar secretion: see below). As in the other 42 examined cases shown in Table 1, the initial colour of this flower disc immediately after anthesis (male phase) was orange (Stage 1, Fig. 1a, RGB: 255, 88, 16), as reported earlier7. Immediately after the petals and stamens fell off, the flower disc colour changed to pink (Stage 2, Fig. 1b, RGB: 255, 82, 102). The styles were not yet elongated at this stage, and the flowers were asexual. In 12 cases with damaged styles and brown stigmas, the flower discs remained pink until the flowers fell off (shown as “O–P” in Table 1).Figure 1Colour change process of a flower disc in Causonis japonica (flower ID 35 in Table 1). Disc colour changed from orange (a, RGB: 255, 88, 16) to pink (b, RGB: 255, 82, 102) before recovering to orange (c, RGB: 255, 88, 16) again, then pink (d, RGB: 255, 120, 94) again. In the last stage, the flower disc turned brownish pink (e, RGB: 232, 162, 169) then fell off. The two orange colour stages were synchronised with flower sexual activity. (a) First orange stage shows stamen activity (male phase); (c) second orange stage indicates stigma maturation (female phase). Nectar secretion was active only in the orange stages and more active in the female phase; the same tendency was observed in the other cases shown in Table 1.Full size imageHowever, in the remaining 31 flowers with normally elongated styles, maturation of the styles (female phase) coincided with the flower discs again exhibiting a distinct orange colour (Stage 3, Fig. 1c, RGB: 255, 88, 16). After the female phase, the flower discs turned pink again (Stage 4, shown as “O–P–O–P” in Table 1), and brownish colouration appeared in the stigmas (Fig. 1d, RGB: 255, 120, 94). Finally, the flower discs turned to a faded pink (Fig. 1e, RGB: 232, 162, 169) just before the flowers fell off. Therefore, the above observations imply that colour-change has a strict correlation with sexual phase.The timings of the disc colour change to the second orange stage (female phase) varied depending on the onset time of each flower. Most flowers that opened before 10:00 reached the second orange stage (female phase) on the afternoon of the same day (except for two flowers, ID 1 and 2 in Table 1). Conversely, flowers that bloomed after 10:00 reached the second orange stage (female phase) at approximately noon the following day. These flowering processes were not fully synchronised in the same inflorescence; therefore, pink and orange discs often coexisted in the same inflorescence. Indeed, we can collect various stages of flowers at a time point from one population as shown in Fig. 2a.Figure 2Histology of floral discs of C. japonica. (a) Floral disc colour change observed in a triploid individual. Flowers were hand-sectioned along the longitudinal axis to show inside colouration of floral disc. Floral phase was judged from the stigma length and colour of the stigma tip; from left, initial stage with orange floral disc and short style, first pink stage with short style, second orange stage with elongated style with matured stigma, and second orange stage with elongated style. Unit of scale bar = 1 mm. (b–e) Longitudinal sections of floral discs in the initial orange stage (b, d) and pink stage (c, e). Scale bar = 500 µm. (b, c) Hand sections of living floral discs showing pigmentation of vacuoles in some scattered cells. (d, e) Resin-embedded sections of floral discs showing histology.Full size imageFigure 1 also shows the typical time-course changes of nectar activities (flower ID 35 in Table 1), which indicates active nectar secretions during both orange colour stages. That is, the flower discs secreted nectar in both male and female phases, with no visible nectar secretion in the pink stages. Moreover, nectar secretion in the female phase of this flower was higher than that in the male phase, a tendency that was also observed in other flowers; however, the total volume of nectar varied among the flowers shown in Table 1. During anthesis, we confirmed that bees, wasps, ants and other insects visited the flowers as previously described7 (Supplementary Fig. 1).Longitudinal sections of flowers in the pink-coloured and orange-coloured stages (Fig. 2a) revealed that pigmentation occurred only in a subset of parenchymatous cells in both cases (Fig. 2b, c). No structural or cytological changes were observed between the initial orange stage and the pink stages (Fig. 2d, e), suggesting that the observed oscillating colour change depends on the degradation and biosynthesis of orange pigments.To understand what pigments are involved in the dual colour change of the C. japonica flower disc, we extracted carotenoids and chlorophyll with Acetone. Anthocyanin was also extracted with Methanol-HCl. As a result, while anthocyanin content was not significantly altered throughout the stages examined, we found that carotenoid level strongly correlated with the colour change detected by naked eye. More specifically, in stages 1 and 3 the carotenoid content was high (63.8 and 65.3 µg/g dry weight, respectively), but significantly decreased in stages 2 and 4 (14.3 and 36.5 µg/g dry weight, respectively) (Table 2). Interestingly, an increase in chlorophyll content was confined to stage 4 (Table 2). Together, the observed dual colour change was ascribed to the decrease (stage 2) and the increase (stage 3) of carotenoid contents in the flower discs.Table 2 Chlorophyll and Carotenoid contents in the flower discs of C. japonica.Full size table More

Rosenzweig, M. L. Habitat selection and population interactions: the search for mechanism. Am. Nat. 137, S5–S28 (1991).Article 

Google Scholar 
Binckley, C. A. & Resetarits, W. J. Habitat selection determines abundance, richness and species composition of beetles in aquatic communities. Biol. Lett. 1, 370–374 (2005).Article 
PubMed 
PubMed Central 

Google Scholar 
Foltz, S. J. & Dodson, S. I. Aquatic Hemiptera community structure in stormwater retention ponds: A watershed land cover approach. Hydrobiologia 621, 49–62 (2009).Article 

Google Scholar 
Resetarits, W. J. Habitat selection behaviour links local and regional scales in aquatic systems: Habitat selection at multiple spatial scales. Ecol. Lett. 8, 480–486 (2005).Article 
PubMed 

Google Scholar 
Leclerc, M. et al. Quantifying consistent individual differences in habitat selection. Oecologia 180, 697–705 (2016).Article 
PubMed 

Google Scholar 
Morris, D. W. Adaptation and habitat selection in the eco-evolutionary process. Proc. R. Soc. B Biol. Sci. 278, 2401–2411 (2011).Article 

Google Scholar 
Resetarits, W. J. Colonization under threat of predation: avoidance of fish by an aquatic beetle, Tropisternus lateralis (Coleoptera: Hydrophilidae). Oecologia 129, 155–160 (2001).Article 
PubMed 

Google Scholar 
Wellborn, G. A., Skelly, D. K. & Werner, E. E. Mechanisms creating community structure across a freshwater habitat gradient. Annu. Rev. Ecol. Evol. Syst. 27, 337–363 (1996).Article 

Google Scholar 
Klečka, J. & Boukal, D. S. Who eats whom in a pool? A comparative study of prey selectivity by predatory aquatic insects. PLoS ONE 7, e37741 (2012).Article 
PubMed 
PubMed Central 

Google Scholar 
Nilsson, P. A. & Brönmark, C. Prey vulnerability to a gape-size limited predator: behavioural and morphological impacts on northern pike piscivory. Oikos 88, 539–546 (2000).Article 

Google Scholar 
Šigutová, H. et al. Specialization directs habitat selection responses to a top predator in semiaquatic but not aquatic taxa. Sci. Rep. 11, 18928 (2021).Article 
PubMed 
PubMed Central 

Google Scholar 
Pintar, M. R. & Resetarits, W. J. Match and mismatch: integrating consumptive effects of predators, prey traits, and habitat selection in colonizing aquatic insects. Ecol. Evol. 11, 1902–1917 (2021).Article 
PubMed 
PubMed Central 

Google Scholar 
Pintar, M. R. & Resetarits, W. J. Jr. Out with the old, in with the new: oviposition preference matches larval success in cope’s gray treefrog Hyla chrysoscelis. J. Herpetol. 51, 186–189 (2017).Article 

Google Scholar 
Wildermuth, H. Habitat selection and oviposition site recognition by the dragonfly Aeshna juncea (L.): an experimental approach in natural habitats (Anisoptera: Aeshnidae). Odonatologica 22, 27–44 (1993).Fortin, D., Morris, D. W. & McLoughlin, P. D. Habitat selection and the evolution of specialists in heterogeneous environments. Isr. J. Ecol. Evol. 54, 311–328 (2008).Article 

Google Scholar 
McLoughlin, P. D., Boyce, M. S., Coulson, T. & Clutton-Brock, T. Lifetime reproductive success and density-dependent, multi-variable resource selection. Proc. R. Soc. B Biol. Sci. 273, 1449–1454 (2006).Article 

Google Scholar 
Morris, D. W. Scales and costs of habitat selection in heterogeneous landscapes. Evol. Ecol. 6, 412–432 (1992).Article 

Google Scholar 
McLoughlin, P. D., Morris, D. W., Fortin, D., Wal, E. V. & Contasti, A. L. Considering ecological dynamics in resource selection functions. J. Anim. Ecol. 79, 4–12 (2010).Article 
PubMed 

Google Scholar 
Leclerc, M., Dussault, C. & St-Laurent, M.-H. Behavioural strategies towards human disturbances explain individual performance in woodland caribou. Oecologia 176, 297–306 (2014).Article 
PubMed 

Google Scholar 
Bolnick, D. I. et al. The ecology of individuals: incidence and implications of individual specialization. Am. Nat. 161, 1–28 (2003).Article 
MathSciNet 
PubMed 

Google Scholar 
Sheppard, C. E. et al. Intragroup competition predicts individual foraging specialisation in a group-living mammal. Ecol. Lett. 21, 665–673 (2018).Article 
PubMed 
PubMed Central 

Google Scholar 
Forstmeier, W. & Birkhead, T. R. Repeatability of mate choice in the zebra finch: consistency within and between females. Anim. Behav. 68, 1017–1028 (2004).Article 

Google Scholar 
Gómez-Laplaza, L. M. The influence of social status on shoaling preferences in the freshwater angelfish (Pterophyllum scalare). Behaviour 142, 827–844 (2005).Article 

Google Scholar 
Gillingham, M. P. & Parker, K. L. The importance of individual variation in defining habitat selection by moose in northern British Columbia. Alces 44, 7–20 (2008).
Google Scholar 
Lesmerises, R. & St-Laurent, M.-H. Not accounting for interindividual variability can mask habitat selection patterns: a case study on black bears. Oecologia 185, 415–425 (2017).Article 
PubMed 

Google Scholar 
van Beest, F. M. et al. Increasing density leads to generalization in both coarse-grained habitat selection and fine-grained resource selection in a large mammal. J. Anim. Ecol. 83, 147–156 (2014).Article 
PubMed 

Google Scholar 
Fretwell, S. D. & Lucas, H. L. On territorial behavior and other factors influencing habitat distribution in birds I. Theoretical development. Biotheoretica 19, 16–36 (1970).Article 

Google Scholar 
Binckley, C. A. & Resetarits, W. J. Functional equivalence of non-lethal effects: generalized fish avoidance determines distribution of gray treefrog, Hyla chrysoscelis, larvae. Oikos 102, 623–629 (2003).Article 

Google Scholar 
Kraus, J. M. & Vonesh, J. R. Feedbacks between community assembly and habitat selection shape variation in local colonization. J. Anim. Ecol. 79, 795–802 (2010).PubMed 

Google Scholar 
Pollard, C. J. et al. Removal of an exotic fish influences amphibian breeding site selection: Exotic fish removal. J. Wildl. Manag. 81, 720–727 (2017).Article 

Google Scholar 
Calenge, C., Dufour, A. B. & Maillard, D. K-select analysis: a new method to analyse habitat selection in radio-tracking studies. Ecol. Model. 186, 143–153 (2005).Article 

Google Scholar 
Freitas, C., Kovacs, K. M., Lydersen, C. & Ims, R. A. A novel method for quantifying habitat selection and predicting habitat use. J. Appl. Ecol. 45, 1213–1220 (2008).
Google Scholar 
Mitchell, L. J., Kohler, T., White, P. C. L. & Arnold, K. E. High interindividual variability in habitat selection and functional habitat relationships in European nightjars over a period of habitat change. Ecol. Evol. 10, 5932–5945 (2020).Article 
PubMed 
PubMed Central 

Google Scholar 
Richter, L. et al. So close and yet so different: the importance of considering temporal dynamics to understand habitat selection. Basic Appl. Ecol. 43, 99–109 (2020).Article 

Google Scholar 
Tyler, J. A. & Rose, K. A. Individual variability and spatial heterogeneity in fish population models. Rev. Fish Biol. Fish. 4, 91–123 (1994).Article 

Google Scholar 
Córdoba-Aguilar, A. Dragonflies and Damselflies: Model Organisms for Ecological and Evolutionary Research. (Oxford University Press, 2008).Sandall, E. L. & Fischer, B. Be a professional: attend to the insects. Am. Entomol. 65, 176–179 (2019).Article 

Google Scholar 
Blaustein, L. Oviposition site selection in response to risk of predation: evidence from aquatic habitats and consequences for population dynamics and community. in Evolutionary theory and processes: modern perspectives (ed. Wasser, S. P.) 441–456 (Kluwer, 1999).Helebrandová, J. B., Pyszko, P. & Dolný, A. Behavioural phenotypic plasticity of submerged oviposition in damselflies (Insecta: Odonata). Insects 10, 124 (2019).Article 
PubMed 
PubMed Central 

Google Scholar 
Hollis, K. & Guillette, L. What associative learning in insects tells us about the evolution of learning and fixed behavior. Int. J. Comp. Psychol. 28, 25706 (2015).Article 

Google Scholar 
Papaj, D. R. & Lewis, A. C. Insect Learning: Ecological and Evolutinary Perspectives. (Chapman & Hall, 1993).Simons, M. & Tibbetts, E. Insects as models for studying the evolution of animal cognition. Curr. Opin. Insect Sci. 34, 117–122 (2019).Article 
PubMed 

Google Scholar 
Benard, M. F. Predator-induced phenotypic plasticity in organisms with complex life histories. Annu. Rev. Ecol. Evol. Syst. 35, 651–673 (2004).Article 

Google Scholar 
Cook, W. L. & Streams, F. A. Fish predation on Notonecta (Hemiptera): relationship between prey risk and habitat utilization. Oecologia 64, 177–183 (1984).Article 
CAS 
PubMed 

Google Scholar 
Larson, D. J. The predaceous water beetles (Coleoptera: Dytiscidae) of Alberta: Systematics, natural history and distribution. Quaest. Entomol. 11, 245–498 (1985).
Google Scholar 
Svensson, B. G., Tallmark, B. & Petersson, E. Habitat heterogeneity, coexistence and habitat utilization in five backswimmer species (Notonecta spp.; Hemiptera, Notonectidae). Aquat. Insects 22, 81–98 (2000).Lock, K., Adriaens, T., Meutter, F. V. D. & Goethals, P. Effect of water quality on waterbugs (Hemiptera: Gerromorpha & Nepomorpha) in Flanders (Belgium): results from a large-scale field survey. Ann. Limnol. Int. J. Limnol. 49, 121–128 (2013).Article 

Google Scholar 
Macan, T. T. A twenty-one-year study of the water-bugs in a Moorland Fishpond. J. Anim. Ecol. 45, 913–922 (1976).Article 

Google Scholar 
Boukal, D. S. et al. Catalogue of water beetles of the Czech Republic. Klapalekiana 43 (Suppl.), 1–289 (2007).Åbjörnsson, K., Wagner, B. M. A., Axelsson, A., Bjerselius, R. & Olsén, K. H. Responses of Acilius sulcatus (Coleoptera: Dytiscidae) to chemical cues from perch (Perca fluviatilis). Oecologia 111, 166–171 (1997).Article 
PubMed 

Google Scholar 
Gioria, M., Schaffers, A., Bacaro, G. & Feehan, J. The conservation value of farmland ponds: Predicting water beetle assemblages using vascular plants as a surrogate group. Biol. Conserv. 143, 1125–1133 (2010).Article 

Google Scholar 
Bergsten, J. & Miller, K. B. Taxonomic revision of the Holarctic diving beetle genus Acilius Leach (Coleoptera: Dytiscidae): Acilius taxonomic revision. Syst. Entomol. 31, 145–197 (2005).Article 

Google Scholar 
Everard, M. Britain’s Freshwater Fishes. (Princeton University Press, 2013).Miller, K. B. & Bergsten, J. Predaceous diving beetle sexual systems. in Ecology, systematics, and the natural history of predaceous diving beetles (Coleoptera: Dytiscidae) (ed. Yee, D. A.) 199–234 (Springer Netherlands, 2014).Culler, L. E., Ohba, S. & Crumrine, P. Predator-prey interactions of dytiscids. in Ecology, Systematics, and the Natural History of Predaceous Diving Beetles (Coleoptera: Dytiscidae) (ed. Yee, D. A.) 363–379 (Springer, 2014).Baines, C. B., McCauley, S. J. & Rowe, L. Dispersal depends on body condition and predation risk in the semi-aquatic insect Notonecta undulata. Ecol. Evol. 5, 2307–2316 (2015).Article 
PubMed 
PubMed Central 

Google Scholar 
Baines, C. B., Ferzoco, I. M. & McCauley, S. J. Sex-biased dispersal is independent of sex ratio in a semiaquatic insect. Behav. Ecol. Sociobiol. 71, 119 (2017).Article 

Google Scholar 
Hungerford, H. B. The biology and ecology of aquatic and semiaquatic Hemiptera. Univ. Kans. Sci. Bull. 11, 3–334 (1919).
Google Scholar 
Streams, F. A. Intrageneric predation by Notonecta (Hemiptera: Notonectidae) in the laboratory and in nature. Ann. Entomol. Soc. Am. 85, 265–273 (1992).Article 

Google Scholar 
Halekoh, U., Højsgaard, S. & Yan, J. The R Package geepack for generalized estimating equations. J. Stat. Softw. 15, 1–11 (2006).Article 

Google Scholar 
Lenth, R. V. Least-squares means: the R package lsmeans. J. Stat. Softw. 69, 1–33 (2016).Article 

Google Scholar 
Bates, A., Maechler, M., Bolker, B. & Walker, S. Fitting linear mixed-effects models using lme4. J. Stat. Softw. 67, 1–48 (2015).Article 

Google Scholar 
Stoffel, M. A., Nakagawa, S. & Schielzeth, H. rptR: repeatability estimation and variance decomposition by generalized linear mixed-effects models. Methods Ecol. Evol. 8, 1639–1644 (2017).Article 

Google Scholar 
R Core Team. R: A language and environment for statistical computing. (2021).Harvill, M. L. The antipredatory behavior of the aquatic diving beetle, Coptotomus venustus (Say)(Coleoptera: Dytiscidae) in response to fish predation. (Texas A&M University, 1994).McCauley, S. J. & Rowe, L. Notonecta exhibit threat-sensitive, predator-induced dispersal. Biol. Lett. 6, 449–452 (2010).Article 
PubMed 
PubMed Central 

Google Scholar 
Schoeppner, N. M. & Relyea, R. A. Damage, digestion, and defence: the roles of alarm cues and kairomones for inducing prey defences. Ecol. Lett. 8, 505–512 (2005).Article 
PubMed 

Google Scholar 
Roberts, G. Why individual vigilance declines as group size increases. Anim. Behav. 51, 1077–1086 (1996).Article 

Google Scholar 
Giller, P. S. Locomotory efficiency in the predation strategies of the British Notonecta (Hempitera, Heteroptera). Oecologia 52, 273–277 (1982).Article 
CAS 
PubMed 

Google Scholar 
Gittelman, S. H. Locomotion and predatory strategy in backswimmers (Hemiptera: Notonectidae). Am. Midl. Nat. 92, 496–500 (1974).Article 

Google Scholar 
Morris, D. W. Density-dependent habitat selection: testing the theory with fitness data. Evol. Ecol. 3, 80–94 (1989).Article 

Google Scholar 
Holt, R. D. Population dynamics in two-patch environments: Some anomalous consequences of an optimal habitat distribution. Theor. Popul. Biol. 28, 181–208 (1985).Article 
MathSciNet 
MATH 

Google Scholar 
Briers, R. A. Metapopulation ecology of Notonecta in small ponds. Doctoral dissertation. (1999).Popham, E. J. The migration of aquatic bugs with special reference to the Corixidae (Hemiptera Heteroptera). Arch. Für Hydrobiol. 60, 450–496 (1964).
Google Scholar 
Doligez, B., Cadet, C., Danchin, E. & Boulinier, T. When to use public information for breeding habitat selection? The role of environmental predictability and density dependence. Anim. Behav. 66, 973–988 (2003).Article 

Google Scholar 
Pintar, M. R. & Resetarits, W. J. Aquatic beetles influence colonization of disparate taxa in small lentic systems. Ecol. Evol. 10, 12170–12182 (2020).Article 
PubMed 
PubMed Central 

Google Scholar 
Sebastián-González, E., Sánchez-Zapata, J. A., Botella, F. & Ovaskainen, O. Testing the heterospecific attraction hypothesis with time-series data on species co-occurrence. Proc. R. Soc. B Biol. Sci. 277, 2983–2990 (2010).Article 

Google Scholar 
Giller, P. S. & McNeill, S. Predation strategies, resource partitioning and habitat selection in Notonecta (Hemiptera/Heteroptera). J. Anim. Ecol. 50, 789–808 (1981).Article 

Google Scholar 
Buxton, V. L., Enos, J. K., Sperry, J. H. & Ward, M. P. A review of conspecific attraction for habitat selection across taxa. Ecol. Evol. 10, 12690–12699 (2020).Article 
PubMed 
PubMed Central 

Google Scholar 
Ferzoco, I. M. C., Baines, C. B. & McCauley, S. J. Co-occurring Notonecta (Hemiptera: Heteroptera: Notonectidae) species differ in their behavioral response to cues of Belostoma (Hemiptera: Heteroptera: Belostomatidae) predation risk. Ann. Entomol. Soc. Am. 112, 402–408 (2019).Article 

Google Scholar 
Roughgarden, J. Evolution of niche width. Am. Nat. 106, 683–718 (1972).Article 

Google Scholar 
Ruckstuhl, K. E. Sexual segregation in vertebrates: proximate and ultimate causes. Integr. Comp. Biol. 47, 245–257 (2007).Article 
CAS 
PubMed 

Google Scholar 
Hochkirch, A., Gröning, J. & Krause, S. Intersexual niche segregation in Cepero’s ground-hopper Tetrix ceperoi. Evol. Ecol. 21, 727–738 (2007).Article 

Google Scholar 
Romey, W. L. & Wallace, A. C. Sex and the selfish herd: sexual segregation within nonmating whirligig groups. Behav. Ecol. 18, 910–915 (2007).Article 

Google Scholar 
Main, M. B., Weckerly, F. W. & Bleich, V. C. Sexual segregation in ungulates: new directions for research. J. Mammal. 77, 449–461 (1996).Article 

Google Scholar 
Trivers, R. L. Parental investment and sexual selection. in Sexual Selection and the Descent of Man 1871–1971 (ed. Campbell, B.) (Aldine Publishing Company, 1972).Bonduriansky, R. The evolution of male mate choice in insects: a synthesis of ideas and evidence. Biol. Rev. Camb. Philos. Soc. 76, 305–339 (2001).Article 
CAS 
PubMed 

Google Scholar 
Foster, S. E. & Soluk, D. A. Protecting more than the wetland: The importance of biased sex ratios and habitat segregation for conservation of the Hine’s emerald dragonfly Somatochlora hineana Williamson. Biol. Conserv. 127, 158–166 (2006).Article 

Google Scholar 
Miller, K. B. The phylogeny of diving beetles (Coleoptera: Dytiscidae) and the evolution of sexual conflict. Biol. J. Linn. Soc. 79, 359–388 (2003).Article 

Google Scholar 
Watson, P. J., Stallmann, R. R. & Arnqvist, G. Sexual conflict and the energetic costs of mating and mate choice in water striders. Am. Nat. 151, 46–58 (1998).Article 
CAS 
PubMed 

Google Scholar 
Rowe, L., Krupa, J. J. & Sih, A. An experimental test of condition-dependent mating behavior and habitat choice by water striders in the wild. Behav. Ecol. 7, 474–479 (1996).Article 

Google Scholar 
McLain, D. K. & Pratt, A. E. The cost of sexual coercion and heterospecific sexual harassment on the fecundity of a host-specific, seed-eating insect (Neacoryphus bicrucis). Behav. Ecol. Sociobiol. 46, 164–170 (1999).Article 

Google Scholar 
Stone, G. N. Female foraging responses to sexual harassment in the solitary bee Anthophora plumipes. Anim. Behav. 50, 405–412 (1995).Article 

Google Scholar 
Martens, A. & Rehfeldt, G. Female aggregation in Platycypha caligata (Odonata: Chlorocyphidae): A tactic to evade male interference during oviposition. Anim. Behav. 38, 369–374 (1989).Article 

Google Scholar 
Kolar, V. & Boukal, D. S. Habitat preferences of the endangered diving beetle Graphoderus bilineatus: implications for conservation management. Insect Conserv. Divers. 13, 480–494 (2020).Article 

Google Scholar 
Wilcox, C. Habitat size and isolation affect colonization of seasonal wetlands by predatory aquatic insects. Isr. J. Zool. 47, 459–475 (2001).Article 

Google Scholar 
Baines, C. B., Ferzoco, I. M. C. & McCauley, S. J. Phenotype-by-environment interactions influence dispersal. J. Anim. Ecol. 88, 1263–1274 (2019).Article 
PubMed 

Google Scholar 
Liao, W., Venn, S. & Niemelä, J. Diving beetle (Coleoptera: Dytiscidae) community dissimilarity reveals how low landscape connectivity restricts the ecological value of urban ponds. Landsc. Ecol. 37, 1049–1058 (2022).Article 

Google Scholar  More

Stocks, G., Seales, L., Paniagua, F., Maehr, E. & Bruna, E. M. The geographical and institutional distribution of ecological research in the tropics. Biotropica 40, 397–404 (2008).Article 

Google Scholar 
Bernstein, J. M., Murphy, J. C., Voris, H. K., Brown, R. M. & Ruane, S. Phylogenetics of mud snakes (Squamata: Serpentes: Homalopsidae): A paradox of both undescribed diversity and taxonomic inflation. Mol. Phylogenet. Evol. 160, 107109 (2021).Article 
PubMed 

Google Scholar 
Murphy, J. C., Voris, H. K. & Karns, D. R. The dog-faced water snakes, a revision of the genus Cerberus Cuvier, (Squamata, Serpentes, Homalopsidae), with the description of a new species. Zootaxa 3484, 1–34 (2012).Article 

Google Scholar 
Stuart, B. L. The harvest and trade of reptiles at U Minh Thuong National Park, southern Viet Nam. Traffic Bull. 20, 25–34 (2004).
Google Scholar 
Brooks, S. E., Allison, E. H. & Reynolds, J. D. Vulnerability of Cambodian water snakes: Initial assessment of the impact of hunting at Tonle Sap Lake. Biol. Conserv. 139, 401–414 (2007).Article 

Google Scholar 
Murphy, J. C. Homalopsid Snakes (Evolution in the Mud (Krieger Publishing, Malabar, 2007).
Google Scholar 
Karns, D. R., Murphy, J. C. & Voris, H. K. Semi-aquatic snake communities of the central plain region of Thailand. Trop. Nat. Hist. 10, 1–25 (2010).
Google Scholar 
Jayne, B. C., Voris, H. K. & Heang, K. B. Diet, feeding behavior, growth, and numbers of a population of Cerberus rynchops (Serpentes: Homalopsinae) in Malaysia: a contribution in celebration of the distinguished scholarship of Robert F. Inger on the occasion of his sixty-fifth birthday. Fieldiana Zoology, Series 50 (Field Museum of Natural History, Chicago, IL, 1988).Chim, C. K. & Diong, C. H. A mark-recapture study of a dog-faced water snake Cerberus schneiderii (Colubridae: Homalopsidae) population in Sungei Buloh Wetland Reserve Singapore. Raffles Bull. Zool. 61, 811–825 (2013).
Google Scholar 
Shine, R., Ambariyanto, Harlow, P. S. & Mumpuni. Ecological attributes of two commercially-harvested python species in northern Sumatra. J. Herpet. 33, 249–257 (1999).Natusch, D. J., Lyons, J. A., Riyanto, A., Khadiejah, S. & Shine, R. Detailed biological data are informative, but robust trends are needed for informing sustainability of wildlife harvesting: A case study of reptile offtake in Southeast Asia. Biol. Conserv. 233, 83–92 (2019).Article 

Google Scholar 
Natusch, D. J., Lyons, J. A., Riyanto, A. & Shine, R. Harvest effects on blood pythons in North Sumatra. J. Wildl. Manage. 84, 249–255 (2020).Article 

Google Scholar 
Shine, R., Harlow, P. S. & Keogh, J. S. The influence of sex and body size on food habits of a giant tropical snake, Python reticulatus. Funct. Ecol. 12, 248–258 (1988).Article 

Google Scholar 
Shine, R., Harlow, P. S. & Keogh, J. S. The allometry of life-history traits: Insights from a study of giant snakes (Python reticulatus). J. Zool. 244, 405–414 (1998).Article 

Google Scholar 
Shine, R. & Harlow, P. S. Reticulated pythons in Sumatra: biology, harvesting and sustainability. Biol. Conserv. 87, 349–357 (1999).Article 

Google Scholar 
Hoesel, J. K. P. Ophidia Javanica (Museum Zoologicum Bogoriense, Kebun Raya, Indonesia, 1959).Voris, H. K. & Murphy, J. C. The prey and predators of homalopsine snakes. J. Nat. Hist. 36, 1621–1632 (2002).Article 

Google Scholar 
Wall, F. A popular treatise on the common Indian Snakes. Part 26. J. Bombay Nat. Hist. Soc. 26, 89–97 (1918).Gorman, G. C., Licht, P. & McCollum, F. Annual reproductive patterns in three species of marine snakes from the central Philippines. J. Herpetol. 15, 335–354 (1981).Article 

Google Scholar 
Auffenberg, W. The herpetofauna of Komodo, with notes on adjacent areas. Bull. Florida State Mus. Biol. Sci. 25, 39–156 (1980).
Google Scholar 
Alcala, A. C. Guide to Philippine Flora and Fauna. Vol. X. Amphibians and Reptiles (Natural Resource Management Center, Ministry of Natural Resources and the University of the Philippines, Manila, Philippines, 1986).Harlow, P. S. & Taylor, J. E. Reproductive ecology of the jacky dragon (Amphibolurus muricatus): An agamid lizard with temperature-dependent sex determination. Austral. Ecol. 25, 640–652 (2000).Article 

Google Scholar 
Saint Girons, H. & Pfeffer, P. Notes sur l’ecologie des serpents du Cambodge. Zool. Mededelingen 47, 65–87 (1972).Kusrini, M. D. et al. Abundance, demography, and harvesting of water snakes from agricultural landscapes in West Java, Indonesia. Wildl. Res. In review (2022).Shine, R. Sexual differences in morphology and niche utilization in an aquatic snake Acrochordus arafurae. Oecologia 69, 260–267 (1986).Article 
PubMed 

Google Scholar 
Houston, D. & Shine, R. Sexual dimorphism and niche divergence: Feeding habits of the Arafura filesnake. J. Anim. Ecol. 62, 737–748 (1993).Article 

Google Scholar 
Shine, R., Reed, R., Shetty, S. & Cogger, H. Relationships between sexual dimorphism and niche partitioning within a clade of sea-snakes (Laticaudinae). Oecologia 133, 45–53 (2002).Article 
CAS 
PubMed 

Google Scholar 
Vincent, S. E., Herrel, A. & Irschick, D. J. Sexual dimorphism in head shape and diet in the cottonmouth snake (Agkistrodon piscivorus). J. Zool. 264, 53–59 (2004).Article 

Google Scholar 
Perkins, M. W., Cloyed, C. S. & Eason, P. K. Intraspecific dietary variation in niche partitioning within a community of ecologically similar snakes. Evol. Ecol. 34, 1017–1035 (2020).Article 

Google Scholar 
Shine, R. & Goiran, C. Sexual dimorphism in size and shape of the head in the sea snake Emydocephalus annulatus (Hydrophiinae, Elapidae). Sci. Rep. 11, 20026 (2021).Article 
CAS 
PubMed 
PubMed Central 

Google Scholar 
Shine, R. Intersexual dietary divergence and the evolution of sexual dimorphism in snakes. Am. Nat. 138, 103–122 (1991).Article 

Google Scholar 
Bonnet, X., Shine, R., Naulleau, G. & Vacher-Vallas, M. Sexual dimorphism in snakes: Different reproductive roles favour different body plans. Proc. R. Soc. B 265, 179–183 (1998).Article 
PubMed Central 

Google Scholar 
Shine, R., Olsson, M. M., Moore, I. T., LeMaster, M. P. & Mason, R. T. Why do male snakes have longer tails than females?. Proc. R. Soc. B 266, 2147–2151 (1999).Article 
PubMed Central 

Google Scholar  More

Pablo Ignacio Plaza and Sergio Agustín Lambertucci from the National University of Comahue and the Argentine Research Council in Argentina quantified the contribution of vultures to reducing greenhouse gas emissions by developing two contrasting scenarios. The first assumes that all the dead animals that the vultures can consume are disposed of, whereas in the second scenario, the dead animals are left to decompose in the environment without scavengers. The results show that the current vulture population can reduce emissions by up to 60.7 teragrams CO2 equivalent per year. A decline in vulture populations decreases their mitigation capacity by 30%. The study highlights that vultures are essential to keep our climate cool. More

Rearing O. biroi pupae in social isolation and collecting pupal fluidIn O. biroi colonies, larvae and pupae develop in discrete and synchronized cohorts26. Ten days after the first larvae had entered pupation in a large stock colony, the entire colony was anaesthetized using a CO2 pad, and white pupae were separated using a paintbrush. Pupae were individually placed in 0.2 ml PCR tubes with open lid. These tubes were then placed inside 1.5 ml Eppendorf tubes with 5 µl sterile water at the bottom to provide 100% relative humidity. The outer tubes were closed and kept in a climate room at 25 °C. The inner tube in this design prevents the pupa from drowning in the water reservoir. The outer tubes were kept closed throughout the experiment, except for once a day when the tubes were opened to remove pupal social fluid. Pulled glass capillaries were prepared as described elsewhere29, and used to remove and/or collect secretion droplets. We were careful to leave no remains of the secretion behind on the pupae or the inside of the tubes. To ensure that all secretion had been removed, pupae were taken out of the tube after fluid collection and briefly placed on a tissue paper to absorb any excess liquid. The inner tubes were replaced if needed—for example, if fluid traces were visible on the old tube after collection. Each pupa was checked daily for secretion (absent or present), onset of melanization and eclosion, and whether the pupa was alive (responding to touch). Control groups of 30 pupae and 30 adult ants from the same stock colony and cohort as the isolated pupae were placed in Petri dishes with a plaster of Paris floor, and the same parameters as for the isolated pupae were scored daily. Experiments ended when all pupae had either eclosed or died. Newly eclosed (callow) workers moved freely inside the tube and showed no abnormalities when put in a colony. A pupa was declared dead if it did not shed its pupal skin and did not respond to touch three days after all pupae in the control group had eclosed.To calculate the average secretion volume per secreting pupa (Fig. 1d), the total volume collected daily from a group of isolated pupae (142–166 pupae) was divided by the number of pupae from which fluid had been collected that day. The total volume was determined by multiplying the height of the fluid’s meniscus in the capillary by πr², where r is the inner radius of the capillary (0.29 mm). While pupae were secreting, pupal whole-body wash samples were collected daily. The pupae were removed from colonies with adults and washed promptly with 1500 µl LC–MS grade water. Whole-body wash samples were lyophilized and reconstituted in 15 µl LC–MS grade water.Collecting additional ant species and honeybees, rearing pupae in social isolation, and collecting pupal fluidsColonies of the ants N. flavipes, T. sessile, P. pennsylvanica and Lasius neoniger were collected in NY state, USA (Central Park, Manhattan; Pelham Bay Park, Bronx; Prospect Park, Brooklyn; and Woodstock). Solenopsis invicta colonies were collected in Athens, GA, USA. M. mexicanus colonies were collected in Piñon Hills, CA, USA. Colonies comprised of queens, workers and brood were maintained in the laboratory in airtight acrylic boxes with plaster of Paris floors. Colonies were fed a diet of insects (flies, crickets and mealworms). White pupae were socially isolated, cocoons were removed in the case of P. pennsylvanica, and secretion droplets were collected from melanized pupae as described for O. biroi. A. mellifera pupae of unknown age were socially isolated from hive fragments (A&Z Apiaries, USA) and reared as described for O biroi, except that the rearing temperature was set to 32 °C. Relative humidity was set to either 100% to replicate conditions used for the different ant species, or to 75% as recommended in the literature30.Injecting dye and tracking pupal fluidInjection needles were prepared as in previous studies31. Injections were performed using an Eppendorf Femtojet with a Narishige micromanipulator. The Femtojet was set to Pi 1000 hPa and Pc 60 hPa. Needles were broken by gently touching the capillary tip to the side of a glass slide. To inject, melanized pupae were placed on ‘Sticky note’ tape (Post-it), with the abdomen tip forward and the ventral side upward. Pupae were injected with blue food colouring (McCormick) into the exuvium for 1–2 s by gently piercing the pupal case at the abdominal tip with the needle. During successful injections, no fluid was discharged from the pupa when the needle was removed, and the moulting fluid inside the exuvium was immediately stained. Pupae were washed in water three times to remove any excess dye. Following injections, 10 pupae were reared in social isolation to confirm the secretion of dyed droplets. For experiments, injected pupae were transferred to colonies with adult ants (Figs. 1f and  4c) or to colonies with adult ants and larvae (Figs. 3b and  4c) to track the distribution of the pupal social fluid.After spending 24 h with dye-injected pupae, adults were taken out of the colony, briefly immersed in 95% ethanol, and transferred to PBS. Digestive systems were dissected in cold PBS and mounted in DAKO mounting medium. Crop and stomach images (Fig. 1f, inset and Fig. 4c, inset) were acquired with a Revolve microscope (Echo). Larvae are translucent, and the presence of dye in the digestive system can be assayed without dissection. Whole-body images of larvae were acquired with a Leica Z16 APO microscope equipped with a Leica DFC450 camera and Leica Application Suite version 4.12.0 (Leica Microsystems). In the experiment on larval growth (Fig. 3c), larval length was measured from images using ImageJ32.Occluding pupaeTen pupae were placed on double-sided tape on a glass coverslip with the ventral side up. The area between the pupae was covered with laser-cut filter paper to prevent adults from sticking to the tape. The pupae were then placed in a 5 cm diameter Petri dish with a moist plaster of Paris floor. To block pupal secretion, the tip of the gaster was occluded with a drop of oil-paint (Uni Paint Markers PX-20), which has no discernible toxic effect7. Secreting pupae received a drop of the same paint on their head to control for putative differences resulting from the paint. Pupae were left in isolation for one day before adults were added to the assay chamber.Behavioural tracking of adult preference assayVideos were recorded using BFS-U3-50S5C-C: 5.0 MP, 35 FPS, Sony IMX264, Colour cameras (FLIR) and the Motif Video Recording System (Loopbio). To assess adult preference (Fig. 1g), physical contact of adults with pupae was manually annotated for the first 10 min after the first adult had encountered (physically contacted) a pupa.Protein profilingWe extracted 30 µl of pupal social fluid and whole-body wash samples with 75:25:0.2 acetonitrile: methanol: formic acid. Extracts were vortexed for 10 min, centrifuged at 16,000g and 4 °C for 10 min, dried in a SpeedVac, and stored at −80 °C until they were analysed by LC–MS/MS.Protein pellets were dissolved in 8 M urea, 50 mM ammonium bicarbonate, and 10 mM dithiothreitol, and disulfide bonds were reduced for 1 h at room temperature. Alkylation was performed by adding iodoacetamide to a final concentration of 20 mM and incubating for 1 h at room temperature in the dark. Samples were diluted using 50 mM ammonium bicarbonate until the concentration of urea had reached 3.5 M, and proteins were digested with endopeptidase LysC overnight at room temperature. Samples were further diluted to bring the urea concentration to 1.5 M before sequencing-grade modified trypsin was added. Digestion proceeded for 6 h at room temperature before being halted by acidification with TFA and samples were purified using in-house constructed C18 micropurification tips.LC–MS/MS analysis was performed using a Dionex3000 nanoflow HPLC and a Q-Exactive HF mass spectrometer (both Thermo Scientific). Solvent A was 0.1% formic acid in water and solvent B was 80% acetonitrile, 0.1% formic acid in water. Peptides were separated on a 90-minute linear gradient at 300 nl min−1 across a 75 µm × 100 mm fused-silica column packed with 3 µm Reprosil C18 material (Dr. Maisch). The mass spectrometer operated in positive ion Top20 DDA mode at resolution 60 k/30 k (MS1/MS2) and AGC targets were 3 × 106/2 × 105 (MS1/MS2).Raw files were searched through Proteome Discoverer v.1.4 (Thermo Scientific) and spectra were queried against the O. biroi proteome using MASCOT with a 1% FDR applied. Oxidation of M and acetylation of protein N termini were applied as a variable modification and carbamidomethylation of C was applied as a static modification. The average area of the three most abundant peptides for a matched protein33 was used to gauge protein amounts within and between samples.Functional annotation and gene ontology enrichmentTo supplement the current functional annotation of the O. biroi genome34, the full proteome for canonical transcripts was retrieved from UniProtKB (UniProt release 2020_04) in FASTA format. We then applied the EggNog-Mapper tool35,36 (http://eggnog-mapper.embl.de, emapper version 1.0.3-35-g63c274b, EggNogDB version 2) using standard parameters (m diamond -d none –tax_scope auto –go_evidence non-electronic –target_orthologs all –seed_ortholog_evalue 0.001 –seed_ortholog_score 60 –query-cover 20 –subject-cover 0) to produce an expanded annotation for all GO trees (Molecular Function, Biological Process, Cellular Components). The list of proteins identified in the pupal fluid was evaluated for functional enrichment in these GO terms, P-values were adjusted with an FDR cut-off of 0.05, and the network plots were visualized using the clusterProfiler package37.Metabolite profilingFor bulk polar metabolite profiling, we used 10 µl aliquots of pupal social fluid and whole-body wash (pooled samples). For the time-series metabolite profiling, 1 µl of pupal social fluid and whole-body wash was used. Samples were extracted in 180 µl cold LC–MS grade methanol containing 1 μM of uniformly labelled 15N- and 13C-amino acid internal standards (MSK-A2-1.2, Cambridge Isotope Laboratories) and consecutive addition of 390 µl LC–MS grade chloroform followed by 120 µl of LC–MS grade water.The samples were vortexed vigorously for 10 min followed by centrifugation (10 min at 16,000g and 4 °C). The upper polar metabolite-containing layer was collected, flash frozen and SpeedVac-dried. Dried extracts were stored at −80 °C until LC–MS analysis.LC–MS was conducted on a Q-Exactive benchtop Orbitrap mass spectrometer equipped with an Ion Max source and a HESI II probe, which was coupled to a Vanquish UPLC system (Thermo Fisher Scientific). External mass calibration was performed using the standard calibration mixture every three days.Dried polar samples were resuspended in 60 µl 50% acetonitrile, and 5 µl were injected into a ZIC-pHILIC 150 × 2.1 mm (5 µm particle size) column (EMD Millipore). Chromatographic separation was achieved using the following conditions: buffer A was 20 mM ammonium carbonate, 0.1% (v/v) ammonium hydroxide (adjusted to pH 9.3); buffer B was acetonitrile. The column oven and autosampler tray were held at 40 °C and 4 °C, respectively. The chromatographic gradient was run at a flow rate of 0.150 ml min−1 as follows: 0–22 min: linear gradient from 90% to 40% B; 22–24 min: held at 40% B; 24–24.1 min: returned to 90% B; 24.1 −30 min: held at 90% B. The mass spectrometer was operated in full-scan, polarity switching mode with the spray voltage set to 3.0 kV, the heated capillary held at 275 °C, and the HESI probe held at 250 °C. The sheath gas flow was set to 40 units, the auxiliary gas flow was set to 15 units. The MS data acquisition was performed in a range of 55–825 m/z, with the resolution set at 70,000, the AGC target at 10 × 106, and the maximum injection time at 80 ms. Relative quantification of metabolite abundances was performed using Skyline Daily v 20.1 (MacCoss Lab) with a 2 ppm mass tolerance and a pooled library of metabolite standards to confirm metabolite identity (via data-dependent acquisition). Metabolite levels were normalized by the mean signal of 8 heavy 13C,15N-labelled amino acid internal standards (technical normalization).The raw data were searched for a targeted list of ~230 polar metabolites and the corresponding peaks were integrated manually using Skyline Daily software. We were able to assign peaks to 107 compounds based on high mass accuracy ( More

arising from Y. Wang et al. Nature https://doi.org/10.1038/s41586-021-04016-x (2021)A unique challenge for environmental DNA (eDNA)-based palaeoecological reconstructions and extinction estimates is that organisms can contribute DNA to sediments long after their death. Recently, Wang et al.1 discovered mammoth eDNA in sediments that are between approximately 4.6 and 7 thousand years (kyr) younger than the most recent mammoth fossils in North America and Eurasia, which they interpreted as mammoths surviving on both continents into the Middle Holocene epoch. Here we present an alternative explanation for these offsets: the slow decomposition of mammoth tissues on cold Arctic landscapes is responsible for the release of DNA into sediments for thousands of years after mammoths went extinct. eDNA records are important palaeobiological archives, but the mixing of undatable DNA from long-dead organisms into younger sediments complicates the interpretation of eDNA, particularly from cold and high-latitude systems.All animal tissues, including faeces, contribute DNA to eDNA records2, but the durations across which tissues can contribute genetic information must vary depending on tissue type and local rates of destruction and decomposition. On high-latitude landscapes, soft tissues and skeletal remains of large mammals may persist, unburied, for millennia3,4,5. For example, unburied antlers of caribou (Rangifer tarandus) from Svalbard (Norway) and Ellesmere Island (Canada) have been dated3,4 to between 1 and 2 cal kyr bp (calibrated kyr before present). Elephant seal (Mirounga leonina) remains near the Antarctic coastline5,6 can persist for more than 5,000 years. This is in contrast to bones in warmer settings, which persist for only centuries or decades7,8. Because bones are particularly resistant to decay, quantifying how their persistence changes across environments enables us to constrain the durations that dead individuals generally contribute to eDNA archives. To do this, we consolidated data on the oldest radiocarbon-dated surface-collected bones from different ecosystems. We included bones that we are reasonably confident persisted without being completely buried (‘never buried’), and bones for which exhumation cannot be confidently excluded (‘potentially never buried’). Pairing bone persistence with mean annual temperatures (MAT) from their sample localities, we find a strong link between the local temperature and the logged duration of bone persistence (Fig. 1, never buried bones: R2 = 0.94, P  More