Philippa Kaur

An experiment in secondary forests in the Democratic Republic of the Congo finds that calcium, an overlooked soil nutrient, is scarcer than phosphorus, and represents a potentially greater limitation on tropical forest growth.Ecology can reveal distributional patterns and dynamics in nature. One approach used is studying the elemental composition of plants, which has been linked to ecological processes such as growth, diversity or water use efficiency. More recently, elemental composition has been detected as a cofactor in governing the carbon sink capacity of plants, and therefore climate change mitigation1,2,3. This discovery has added an extra layer of urgency to the field, which now aims to better understand and predict global change. The study of nitrogen and/or phosphorus has until now received most of the attention of plant ecologists: nitrogen is the most abundant element in dry leaves after hydrogen and carbon, forms the main structure of proteins and is strongly linked to photosynthesis4. Phosphorus represents around one-tenth of nitrogen abundance in leaves and is key in energy storage and nucleic acids. However, although these represent only two of the many chemical elements that are in flux throughout ecosystems, whether others may have a dominant role in ecosystem dynamics is an open question. Writing in Nature Ecology & Evolution, Bauters et al.5 share some evidence to motivate broadening out from the dominant focus on nitrogen and phosphorus in terrestrial ecology, by revealing a crucial limiting role of calcium in the dynamics of tropical forest succession. More

Virtual hydrostatic model parametersVarious morphological characteristics were held constant in order to isolate and manipulate the variable of conch shape. A CT-scanned Nautilus pompilius conch was essentially morphed into ammonoid-like conch shapes, populating the Westermann morphospace22 while holding constant septal morphology, septal spacing, and shell/septal thicknesses (Fig. 9). Furthermore, body chamber proportions were determined by iteratively computing soft body volumes that yield Nautilus-like chamber liquid (~ 12% of the phragmocone volume retained)67,68. Septal spacing was measured as the angle from the ventral attachment of the current and previous septa, and the spiraling axis of the conch. Because septal spacing differs in early ontogeny (Fig. S11), only measurements from the 7th to 33rd (terminal) septum were considered. The average angle of 23.46° ± 3.32° (standard deviation) was rounded to 23° and held constant throughout the ontogeny of the hydrostatic models.Figure 9Hydrostatic models of theoretical planispiral cephalopods. These models were constructed by morphing a Nautilus pompilius conch into ammonoid shapes (see “Methods”): (a) oxycone, (b) serpenticone, (c) sphaerocone, and (d) morphospace center. The centers of buoyancy and mass are denoted by the tips of the blue (upper) and red (lower) cones. Prime symbols (′) refer to transparent, transverse views of each respective conch shape. (e) Westermann morphospace22 showing relative positions of these conch shapes. All models were rendered in MeshLab76.Full size imageShell and septal thicknesses were measured with digital calipers from a physical specimen of Nautilus pompilius (Table S13). These measurements were recorded as a ratio of inner whorl height (measured from the ventral point on the current whorl to the ventral point on the previous whorl). These ratios were used in the theoretical models to define shell and septum thicknesses (3.1% of inner whorl height for shell thickness and 2.1% of inner whorl height for septal thickness; Table S13).Hydrostatic model constructionThe near-endmember models were constructed from representative ammonoid specimens (Sphenodiscus lobatus and S. lenticularis—oxycone; Dactylioceras commune—serpenticone; Goniatites crenistria—sphaerocone). Lateral and transverse views were measured from figured specimens for the oxycone (Fig. 5 of Kennedy et al.69), serpenticone (Fig. 2 of Kutygin and Knyazev70), and sphaerocone (Figs. 17 and 20 of Korn and Ebbighausen71). These models were constructed with array algorithms similar to earlier hydrostatic models9,35,72, which were used in a piecewise manner to account for allometric changes in coiling throughout ontogeny (Table S14). These arrays replicated the adult whorl section backwards and translated, rotated, and scaled each successive one. These whorl sections were bridged together to create a single tessellated surface representing the outer interface of the shell. Shell thickness was defined by shrinking the original whorl section so that the thickness between the two was equal to 3.1% of the inner whorl height (Table S13), then using the same array to build the internal interface of the shell. The morphospace center was constructed from previously used conch measurements18 and averaging the whorl section shape in blender (Fig. S12). The corresponding Westermann morphospace parameters (Fig. S13) for each morphology are reported in Table S15.Virtual models of the septa were derived from the CT-scan of Nautilus pompilius (Fig. S14). A single septum was isolated from the adult portion of the phragmocone then smoothed to delete the siphuncular foramen. This septum was placed within the whorl section of each theoretical model and stretched in the lateral directions until it approximately fit. The “magnetize” tool in Meshmixer (Autodesk Inc.) was used to attach the septal margin to the new whorl section so that the Nautilus suture was transferred to the new whorl section. The septum was then smoothed to reconcile the first order curves with the new location of the septal margin. The respective septum for each theoretical model was then replicated with the same array instructions used to build the shell. Because each replicated object was rotated one degree (Table S14), 22 septa were deleted in between every two so that the septal spacing was equal to 23° (Fig. S11).For each theoretical model, the septa were unified with the model of the shell using Boolean operations in Netfabb (Autodesk Inc.). To perform hydrostatic calculations, virtual models must be created for each material of unique density. The virtual model of the shell constrains the shape of the soft body (within the body chamber) and chamber volumes (within the phragmocone). These internal interfaces were isolated from the model of the shell, then their faces inverted for proper, outward-facing orientations of their normals. A conservative soft body estimate was created, aligning with previously published reconstructions64,65,73. The profile shape of this soft body was scaled and maintained between each model. External interfaces of the shell and soft body were also isolated to create a model of the water displaced by each theoretical cephalopod. Each of these models are necessary for hydrostatic calculations (buoyancy and the distribution of organismal mass).Each hydrostatic model is stored in an online repository (Dataset S1; https://doi.org/10.5281/zenodo.5684906). The hydrostatic centers of each virtual model and their volumes and masses are listed in Tables S16 and S17.Hydrostatic calculationsEach theoretical model was scaled to have equal volume (near one kilogram; 0.982 kg–a result of arbitrarily scaling the sphaerocone model to 15 cm in conch diameter). An object is neutrally buoyant when the sum of organismal mass is equal to the mass of water displaced (the principle of Archimedes). The percentage of chamber liquid can be computed to satisfy this condition.$${Phi } = frac{{left( {frac{{{text{V}}_{{{text{wd}}}} {uprho }_{{{text{wd}}}} – {text{V}}_{{{text{sb}}}} {uprho }_{{{text{sb}}}} – {text{V}}_{{{text{sh}}}} {uprho }_{{{text{sh}}}} }}{{{text{V}}_{{{text{ct}}}} }}} right) – left( {{uprho }_{{{text{cl}}}} } right)}}{{left( {{uprho }_{{{text{cg}}}} – {uprho }_{{{text{cl}}}} } right)}}$$
(1)
where Vwd and ρwd are the volume and density of the water displaced, Vsb and ρsb are the volume and density of the soft body, Vsh and ρsh are the volume and density of the shell, ρcl is the density of cameral liquid, ρcg is the density of cameral gas, and Vct is the total volume of all chambers. A soft body density of 1.049 g/cm3 is used based on bulk density calculations of Nautilus-like tissues74, a seawater-filled mantle cavity, and thin calcitic mouthparts21. A shell density of 2.54 g/cm374, cameral liquid density of 1.025 g/cm375, and cameral gas density of 0.001 g/cm3 are adopted from recent hydrostatic studies.Other hydrostatic properties depend on the relative positions of the centers of buoyancy and mass. The center of buoyancy is equal to the center of volume of water displaced. This center and the centers of each virtual model of unique density were computed in the program MeshLab76. The individual centers for each organismal model (soft body, shell, cameral liquid and cameral gas) were used to compute the total center of mass, with an average weighted by material density:$$M = frac{{sum left( {L*m_{o} } right)}}{{sum m_{o} }}$$
(2)
where M is the total center of mass in a principal direction, L is the center of mass of a single object measured with respect to an arbitrary datum in each principal direction, and (m_{o}) is the mass of each object with unique density. Equation 2 was used in the x, y, and z directions to compute the 3D coordinate position of the center of mass. The centers of mass for the chamber contents (liquid and gas) were set equal to the center of volume of all chambers, a minor assumption considering the capillary retention of liquid around the septal margins in the living animals62.The hydrostatic stability index (St) is computed from the relative location of the centers of buoyancy (B) and mass (M), normalized by the cube root of volume (V) for a dimensionless metric that is independent of scale:$$S_{t} = frac{{ sqrt {left( {B_{x} – M_{x} } right)^{2} + left( {B_{y} – M_{y} } right)^{2} + left( {B_{z} – M_{z} } right)^{2} } }}{{sqrt[3]{V}}}$$
(3)
where the subscripts correspond to the x, y, and z components of each hydrostatic center.Apertural orientations were measured in blender after orienting each model so that the center of buoyancy was vertically aligned above the center of mass. Apertural angles of 0° correspond to a horizontally facing soft body, while angles of + 90° and − 90° correspond to upward- and downward-facing orientations, respectively.Thrust angles were measured from the hyponome location (ventral edge of the aperture) to the midpoint of the hydrostatic centers, with respect to the horizontal. Thrust angles of 0° infer idealized horizontal backward transmission of energy into movement, while thrust angles of + 90° and − 90° infer more efficient transmission of energy into downward and upward vertical movement, respectively.Biomimetic robot constructionTo isolate the variable of shell shape on swimming capabilities, only the external shape, and static orientation of each virtual hydrostatic model were used to build physical, 3D printed robots. That is, each model has artificially high hydrostatic stability (Tables S3) to nullify the effect of the thrust angle (the angle at which thrust energy passes through the hydrostatic centers and most efficiently transmits energy into movement; Table S4). Less stable morphotypes (e.g., serpenticones and sphaerocones) are more sensitive to the constraints imposed by this hydrostatic property.Space constraints inside each model were determined by first constructing a propulsion system and electronic components that operate the motor. The models use impeller-based water pumps (Figs. 1d and 10a) driven by a brushed DC motor. This system creates a partial vacuum by centrifugal acceleration, drawing water from a “mantle cavity” and expelling it out of a “hyponome”. This system was optimized by iteratively designing models in Blender77, then testing 3D-printed, stand-alone water pumps. After three iterations, a four-blade impeller and gently tapering hyponome (inner diameter at distal end = 6.7 mm) were chosen. The electronic components used to drive the motor consist of an Arduino Pro Micro microcontroller, a motor driver, and two batteries (Fig. 10). A 3.7 V battery operates the microcontroller, and a larger 7.4 V battery supplies power to the motor. Communication is achieved via infrared, allowing specification of the jet pulse duration, number of pulses, and the power level of the motor (using pulse-width modulation; PWM). Each of these electronic components fold into a compact cartridge capable of being plugged into 3D-printed models of each investigated shell shape (Figs. 2 and 10). Each model was designed with brackets to hold the electronics cartridge in place. The sphaerocone had the most severe space constraints, with low conch diameter to volume ratio. After determining the space required for the electronics (Fig. 10) this model was scaled to 15 cm, and all other models were scaled to have similar volumes (with subtle volume differences due to minor differences in soft body shape compared to the hydrostatic models).Figure 10Biomimetic cephalopod robot components. (a) Ventral view of the sphaerocone biomimetic robot (before covering the pump and mantle cavities) with assembled electronics cartridge to the right. (b) View of electronic components that fit into the cartridge. (c) Electronics cartridge placed in robot. These two halves are fit together with wax to create a water-tight seal. Each model component is denoted by letters in circles: A = Arduino microcontroller, B = microcontroller charger / voltage regulator, C = motor driver, D = infrared sensor, E = indicator LED, F = microcontroller battery (3.7 V), G = motor battery (7.4 V), H = brushed motor, I = impeller and water pump cavity, J = electronics cartridge. The colors of annotations correspond to components depicted in Figs. 1 and 2.Full size imageIn addition to having a propulsion system, biomimetic cephalopod robots must also be capable of neutral buoyancy, while assuming the proper orientation in the water. These robots, and their once-living counterparts, each have differing material densities and associated mass distributions for each component. To reconcile these differences, the total mass and total centers of mass for each model were manipulated by controlling the volume and 3D distribution of the 3D-printed PETG (polyethylene terephthalate glycol) thermoplastic. That is, the shape of this material holds each model component in place while correcting for these differences in hydrostatics. The PETG mass required for neutral buoyancy was found by subtracting the mass of every other model component from the mass of the water displaced by the model (i.e., electronics cartridge, bismuth counterweight, liquid, motor, batteries, electronic components, and self-healing rubber; Table S1). This model configuration also allows buoyancy to be fine-tuned in water, compensating for potential density differences between the virtual water and the actual water in the experimental settings. That is, each virtual model accounts for ~ 9 g of internal liquid, but the actual volume of this liquid can be adjusted in the physical robot with a syringe through a self-healing rubber valve (Table S1; Fig. 1).The 3D position of the total center of mass was manipulated by accounting for the local centers of mass of each material of unique density. Materials like the batteries, motor, and electronic components were each assigned bulk density values because they are made up of composite materials. While this is an approximation, their contributions to the total center of mass are low because they account for small fractions of the total model mass (Tables S1 and S2). These components, like all others, were digitally modeled in Blender77 and their volumes and centers of mass were computed in the program MeshLab76. A dense, bismuth counterweight was also modeled, and positioned to artificially stabilize each model (pulling the z component of the total center of mass downward, while maintaining the horizontal components). The virtual model of this counterweight was used to make a 3D-printed mold, allowing a high heat silicone mold to be casted. The bismuth counterweight was cast from this silicone mold and filed/sanded to the dimensions of its virtual counterpart. Hyponomes were oriented horizontally, to yield movement in this direction. To maintain the same static orientation as the virtual model (same x and y center of mass components), the PETG center of mass was computed with the following equation:$$D_{PETG} = frac{{Mmathop sum nolimits_{i = 1}^{n} m_{i} – mathop sum nolimits_{i = 1}^{n} (D_{i} m_{i} )}}{{left( {m_{PETG} } right)}}$$
(4)
where DPETG is the location of the PETG center of mass from an arbitrary datum in each principal direction. M is the total center of mass in a particular principal direction, mi is the mass of each model component, Di is the local center of mass of each model component in a particular principal direction and mPETG is the mass of the PETG required for a neutrally buoyant condition. See Tables S1 and S2 for a list of model components and measurements.Each model was 3D printed with an Ultimaker S5 3D printer using clear (natural) PETG in separate parts, allowing the internal components to be implanted (i.e., brushed DC motors and bismuth counterweights). Each model part was chemically welded together with 100% dichloromethane, with minor amounts of cyanoacrylate glue used to fill seams (e.g., the water pump lid; Fig. 10a). Each final model consists of the main body (housing the water pump, motor, and counterweight), and a “lid” with brackets that house the electronics cartridge (Figs. 2 and 10). The main body and lid were fused together before each experiment by placing wax (paraffin-beeswax blend) along a tongue and groove seam, heating it with a hairdryer, then vigorously squeezing each part together. Surplus wax extruded from the seam was removed and smoothed, producing a water-tight seal.Thrust calibrationEven though each model was designed to have equal mantle cavity and pump cavity volumes, they produced slightly different thrusts. These differences were likely due to variable degrees of friction between the impellers and the surrounding water pumps. To correct for these differences, the thrust produced by each model was measured with a Vernier Dual-Range Force Sensor (0.01 N resolution). Each robot was attached at the hyponome location, through a series of pulleys, and to the sensor with fishing line (Fig. S1; similar to the methods used for living cephalopods78). Force was recorded for 30-s intervals at a sample rate of 0.05 s. During this time, each model was recorded jetting with a 6-s pulse for 15 trials (Fig. S2A). Each trial had initial noise from setting up the model, then peaked randomly when the fishing line became taught, then stabilized after some period of oscillation. Only the stabilized portion of the thrust profile was used to record thrust at 100% voltage for each model (Fig. S2B). The true zero datum was also subtracted from each of these trials. The lowest thrust from each of the models was used as a baseline (serpenticone and oxycone). Each model was recorded again for 15 trials by lowering the motor voltage in increments of 5% until they yielded similar thrusts (0.3 N) to the original serpenticone and oxycone trials (Fig. S2C). The final power levels were then determined for each model and adjusted with pulse-width modulation (PWM) through the microcontroller: serpenticone (100%), oxycone (100%), sphaerocone (95%), and morphospace center (85%).The peak thrust measured for 1 kg extant Nautilus is around 2 N16. The time-averaged thrust during each pulse is around 23% of this value (0.46 N16). This computed value slightly overpredicts observed maximum velocities for this animal (33 cm/s instead of 25 cm/s), so the appropriate time-averaged thrust is probably slightly lower. The motor in the robots quickly reaches its maximum thrust (~ 0.3 N) once initiated then quickly declines after shutting off (Fig. S2). Therefore, the thrust produced by the robots can be treated as a conservative Nautilus-like jet thrust close to the behavior of escape jetting. One-second pulse and refill intervals are also on par with values reported for extant Nautilus16.Robot buoyancyEach of the models were made near neutrally buoyant by adjusting the allotted ~ 9 g of internal liquid with a syringe through a self-healing rubber valve. The single-pulse experiments were performed in an external pool (ranging ~ 23.5 to 26.5 °C). The three-pulse and maneuverability experiments were performed in an internal pool (the Crimson Lagoon at the University of Utah). This internal pool had slightly higher temperatures (~ 28 °C), yielding lower ambient water densities than the virtual water. These conditions required slightly less internal liquid (~ 2–5 g). These differences in internal liquid masses produced negligibly small shifts in mass distributions because they are very small proportions of total robot masses (Table S1).Perfect neutral buoyancy cannot be practically achieved, but this condition can be closely approached. Each of the biomimetic robots experience subtle upward or downward movements of the course of their 5–15 s long trials due to slightly positive or negative buoyancies. Because these differences in buoyancy influence the vertical component of movement, only the horizontal components are considered for discussion. However, a comparison of velocities computed from full, 3D movement (Eq. 5) and restricted 2D components (Eq. 6) reveals that these differences are minor (Figs. S7 and S8). These comparisons demonstrate that model buoyancy did not substantially influence kinematics other than gross trajectories (Figs. 4 and S9).3D motion trackingAfter adjusting buoyancy, each model was positioned underwater with a grabber tool. This tool was fitted with a bundle of fiber-optic cable (Fig. S4) attached to an infrared remote control. Arduino code (Dataset S2) was uploaded to the microcontroller in the robot allowing jet pulse duration, number of pulses, and power to be adjusted with this remote control. After an infrared pulse is received, the motor activates, and activity is indicated by a green LED that illuminates the model from the inside. This light is used to determine time-zero for each trial of motion tracking.After sending an infrared signal, the movement of each model was recorded with a submersible camera rig fitted with two waterproof cameras (Fig. 3). Each of the four models were monitored during a single, one-second jet for at least 9 trials each. Additionally, the laterally compressed morphotypes (serpenticone and oxycone) were monitored during three, one-second pulses for 10 trials each. The inflated morphotypes (sphaerocone and morphospace center) were not able to be monitored over longer distances because they had the tendency to rotate about the vertical axis, obscuring views of the tracking points. In addition to horizontal movement, turning efficiency (maneuverability about the vertical axis) was monitored by directing the cameras with a top-down view of each model. A 90° elbow attachment for the hyponome was fit to each model to investigate the ease or difficulty of rotation. Each model was designed to spin counter-clockwise when viewed from above so that the influence of the motor’s angular momentum was consistent between models.Footage was recorded with two GoPro Hero 8 Black cameras at 4K resolution and 24 (23.975) frames per second, with linear fields of view. Motion tracking was performed with the software DLTdv879 to record the pixel locations of each tracking point (Figs. 1c and S4). These coordinates were transformed into 3D coordinates in meters using the program easyWand580. The tracking points on each model were used for wand calibration because the distances between these sets of points were fixed. Standard deviations of the reproduced tracking point distances of less than 1 cm were considered suitable.The 3D position datasets allowed velocity, acceleration, rocking, to be computed for each experiment. Additionally angular displacement and angular velocity was of interest for the rotation experiments about the vertical axis. Velocity was computed under two scenarios: (1) using the 3D movement direction between each timestep (Eq. 5), and (2) only considering the horizontal movement direction between each time step (Eq. 6). The latter scenario was preferred to nullify the influences of model buoyancies, which were not perfectly neutral and caused some degree of vertical movement.$$V_{i} = frac{{sqrt {left( {x_{i} – x_{i – 1} } right)^{2} + left( {y_{i} – y_{i – 1} } right)^{2} + left( {z_{i} – z_{i – 1} } right)^{2} } }}{{left( {t_{i} – t_{i – 1} } right)}}$$
(5)
$$V_{i} = frac{{sqrt {left( {x_{i} – x_{i – 1} } right)^{2} + left( {y_{i} – y_{i – 1} } right)^{2} } }}{{left( {t_{i} – t_{i – 1} } right)}}$$
(6)
where V and t are velocity and time, and the subscripts i and i −1 refer to the current and previous time steps, respectively. Coordinate components are denoted by x, y, and z at each timestep. The averaged 3D location of both tracking points was used for each model (i.e., midpoints). Note that Eq. (5) uses the 3D form of the Theorem of Pythagoras, whereas Eq. (6) uses the 2D version. Time zero for each trial was defined as the frame where the robot was illuminated by the internal LED, indicating motor activity. Acceleration was modeled by fitting a linear equation to the datapoints during the one-second pulse interval(s) using the curve fitting toolbox in MATLAB R2020A.The artificially high hydrostatic stability of each model was designed to nullify rocking during movement. This behavior was computed for each model during the one-pulse experiments with the following equation:$$theta_{dv} = cos^{ – 1} left( {frac{{left( {z_{2} – z_{1} } right)}}{{sqrt {left( {x_{2} – x_{1} } right)^{2} + left( {y_{2} – y_{1} } right)^{2} + left( {z_{2} – z_{1} } right)^{2} } }}} right) – theta_{tp}$$
(7)
where (theta_{dv}) is the angle deviated from true vertical and (theta_{tp}) is the angle of the tracking points measured from the vertical in a static setting. The subscripts 1 and 2 of the x, y, and z coordinates refer to the anterior and posterior tracking points, respectively.Maneuverability about the vertical axis was determined by computing the angle between the horizontal components of each tracking point. The net angle from the starting angle for each trial was tabulated. Angular velocity was determined by dividing the change in angle between each frame by the frame duration (1/23.975 fps).Links to example motion tracking footage, and robotic models are deposited in an online repository60,61,63 (Dataset S2; https://doi.org/10.5281/zenodo.6180801). More

Sample collection and hydrographySampling was conducted aboard the University of Victoria’s MSV John Strickland either weekly, biweekly or monthly between 11 March 2010 and 15 November 2017 in Saanich Inlet at 48.59°N, 123.50°W (Fig. 1). To standardize measurements and due to biological significance, seawater was collected from the euphotic zone. Sampling depths corresponded to approximately 100, 50, 15, and 1% of the photosynthetically active radiation (PAR) at the surface (Io). These “light” depths were either determined using a CTD-mounted PAR sensor or a Secchi disk. CTD profiles were performed prior to each seawater cast to measure depth, temperature and conductivity of the water column, and PAR, fluorescence, and dissolved oxygen (when available).Seawater from each light depth was collected using Niskin or GO-FLO bottles on either a rosette sampler or an oceanographic wire. When possible, individual samples were collected directly from the Niskin or GO-FLO bottles. When time was not sufficient to allow direct sampling, bulk samples of seawater from each depth were collected into acid-washed polyethylene carboys, kept cold in the dark, and homogenized before sub-sampling for the individual measurements.Dissolved nutrientsFor the measurements of nitrite (NO2−), nitrate and nitrite (NO3− + NO2−), phosphate (PO43−) and Si(OH)4, seawater samples from each light depth were syringe-filtered through a combusted 0.7 µm (nominal porosity) glass fibre filter into acid-washed 30-mL polypropylene bottles and immediately frozen. All nutrient samples were stored at −20 °C until analysis. Concentrations of NO2−, NO3− + NO2−, PO43−, and Si(OH)4 were determined using an Astoria Nutrient Autoanalyzer (Astoria-Pacific, OR, USA) following the methodology of Barwell-Clarke and Whitney22. During 2014 and 2015, samples for the measurement of Si(OH)4 were collected separately from those for the other nutrients, filtered with a 0.6 µm polycarbonate membrane filter and stored at 4 °C. During this period, Si(OH)4 concentrations were determined manually using the molybdate blue colorimetric methodology23. Replicate (2 or 3) nutrient samples were taken at each depth; average data are presented in the published dataset and the figures (Fig. 2).Fig. 2Dissolved macronutrient concentrations in the euphotic zone of Saanich Inlet from March 11, 2010 to November 15, 2017. Left panels show depth-integrated concentrations (black bars on top) and time-series profiles (filled contour/scatter plots on bottom) for (A) nitrate plus nitrite (NO3− + NO2−), (B) phosphate (PO43−) and (C) silicic acid (Si(OH)4). In the time-series profiles, 2012–2013 data are not interpolated due to single-depth sampling. Grey shaded regions in top panels indicate the phytoplankton growing seasons considered for this study (March 1st – October 30th). Right panels show monthly-averaged depth profiles for the entire 8-year period, illustrating euphotic zone seasonality for each nutrient. The color scale bars on the far right apply to both the time-series vertical profiles and the 8-year seasonal plots. Sampling depths are indicated by round symbols. The year labels are positioned under the tick marks corresponding to January.Full size imageSuspended particulate matterTotal chlorophyll-aChlorophyll-a (Chl-a) was used as a proxy for phytoplankton biomass (Fig. 3A). For total Chl-a analysis, seawater samples (0.25–1 L) were gently vacuum filtered onto 0.7 µm (nominal porosity) glass fiber filters, which were then stored at −20 °C until analysis. Chl-a concentrations were determined using the acetone extraction and acidification method24,25. Acidification of samples decreased the likelihood of overestimation of Chl-a concentrations due to the presence of chlorophyll degradation products26. Filters were submerged in 10 mL of 90% acetone, sonicated for 10 minutes in an ice bath, and left to extract at −20 °C for 22 h. Following the extraction period, samples were allowed to equilibrate to room temperature (~2 h). Fluorescence of the acetone solution containing the extracted Chl-a was measured before and after acidification with 1.2 N hydrochloric acid using a Turner 10-AU fluorometer. The final concentrations of total Chl-a were calculated from measurements made before (Fo) and after (Fa) acidification using Eq. (1)25. The coefficient (τ) of Eq. (1), adapted from Strickland and Parsons25, was derived from a calibration of the Turner 10-AU fluorometer with known pure chlorophyll standards (Table 2).$${rm{C}}{rm{h}}{rm{l}} mbox{-} {rm{a}}(mu g,{L}^{-1})=frac{tau }{tau -1}ast ({rm{F}}{rm{o}}-{rm{F}}{rm{a}})ast 0.814ast left(frac{{rm{V}}{rm{o}}{rm{l}}.{rm{A}}{rm{c}}{rm{e}}{rm{t}}{rm{o}}{rm{n}}{rm{e}},{rm{e}}{rm{x}}{rm{t}}{rm{r}}{rm{a}}{rm{c}}{rm{t}}{rm{e}}{rm{d}}}{{rm{V}}{rm{o}}{rm{l}}.{rm{S}}{rm{e}}{rm{a}}{rm{w}}{rm{a}}{rm{t}}{rm{e}}{rm{r}},{rm{f}}{rm{i}}{rm{l}}{rm{t}}{rm{e}}{rm{r}}{rm{e}}{rm{d}}}right)$$
(1)
Fig. 3Biological particulate concentrations in the euphotic zone of Saanich Inlet from March 11, 2010 – November 15, 2017. Left panels show time-series profiles (filled contour/scatter plots) of (A) total chlorophyll-a (Total Chl-a), (B) particulate carbon (PC), (C) particulate nitrogen (PN), and (D) particulate biogenic silica (bSiO2). The 2012–2013 data are not interpolated due to single-depth sampling. In A, the bar plot in the top panel shows percent contribution of different size fractions to total Chl-a. In (B–D), black bars in top panels show depth-integrated concentrations. Grey shaded regions in bar plots indicate phytoplankton growing seasons considered for this study (March 1st – October 30th). Right panels show monthly-averaged depth profiles for the entire 8-year period, illustrating euphotic zone seasonality for each particulate. The color scale bars on the far right apply to both the time-series vertical profiles and the 8-year seasonal plots. Sampling depths are indicated by round symbols. The year labels are positioned under the tick marks corresponding to January.Full size imageSize fractionated chlorophyll-aTo determine the percent contributions of “pico” (0.7–2 µm), “small nano” (2–5 µm), “large nano” (5–20 µm) and “micro” ( >20 µm) phytoplankton to total Chl-a, seawater samples (0.25–1 L) separate from those used for total Chl-a) were consecutively filtered through 20, 5 and 2 µm polycarbonate membrane filters and 0.7 µm (nominal porosity) glass fiber filters. Between 2013 and 2017, the “pico” and “small nano” size classes were collected as one fraction (0.7–5 µm). Analysis of Chl-a concentrations for each size fraction followed the same procedure outlined for total Chl-a.Particulate carbon and nitrogenParticulate C and N measurements were obtained from seawater samples incubated for carbon (ρC) and nitrate uptake (ρNO3) rates (see section on “Uptake rates of carbon and nitrate” for methodology) (Fig. 3B,C). PC and PN measurements presented in this dataset were taken at the end of ρC and ρNO3 incubations; however, original (‘ambient’) values can be back calculated by subtracting the amount of C and N taken up during the incubation period from the final PC and PN values. The differences between after-incubation PC and PN data and back-calculated ambient values were not significantly different than the measurement error.Particulate biogenic silicaParticulate biogenic silica was used as a proxy for siliceous phytoplankton biomass (Fig. 3D). Seawater samples (0.5–1 L) from each depth were gently vacuum filtered through 0.6 µm polycarbonate membrane filters. Filters were folded and placed in polypropylene centrifuge tubes, dried for 48 h at 60 °C, and then stored in a desiccator at room temperature until analysis. Filters were digested with 4 mL of 0.2 M NaOH for 30–45 min in a water bath at 95 °C27. After digestion, samples were neutralized with 0.1 N HCl and cooled rapidly in an ice bath. Samples were centrifuged to separate out the undissolved lithogenic silica, and colorimetric analysis was performed on the supernatant. The transmittance of the samples, standards, and reverse-order reagent blanks were read at 820 nm using a Beckman DU 530 ultraviolet-visible (UV/Vis) spectrophotometer27,28.Uptake rates of carbon and nitrateSeawater samples (~0.5–1 L) were gently collected into clear polycarbonate bottles. One additional sample was collected from the 100% light depth, into a dark polycarbonate bottle, which did not allow light penetration. After the addition of the isotopic tracers (see below), bottles were placed into an acrylic incubator with constant seawater flow to maintain surface seawater temperature. Three acrylic tubes wrapped in colored and neutral density photo-film (to obtain 50, 15, and 1% of surface PAR) were used to incubate sample bottles under the same in-situ light conditions from which samples were collected. Samples from the 100% light level were placed inside the same acrylic incubator, but outside of the film-covered tubes. A LI-COR® LI-190 Quantum sensor was installed next to the incubator and continuously recorded incoming PAR for the entire incubation period. During sampling in 2010 and 2011, all experiments were performed using a shipboard incubator. For sampling from 2012 onwards, all experiments were done using an incubator on land (University of Victoria Aquatic Facility), which was connected to a seawater system maintained at local surface seawater temperature (approximately 9–12 °C depending on the time of year).Rates of C (ρC) and NO3 (ρNO3) uptake were determined using a stable isotope tracer-technique29,30 (Fig. 4). A single seawater sample from each light depth received a dual spike, with NaH13CO3 (99% 13C purity, Cambridge Isotope Laboratories) for the determination of ρC and Na15NO3 (98 + % 15N purity, Cambridge Isotopes Laboratories) for the determination of ρNO3. Isotope additions were made at approximately 10% of ambient dissolved inorganic carbon (DIC) and NO3− concentrations.Fig. 4Carbon and nitrate uptake rates in the euphotic zone of Saanich Inlet from March 11, 2010 – November 15, 2017. Left panels show depth-integrated rates (black-bars on top) and time-series profiles (filled contour/scatter plots below) of (A) carbon (ρC) and (B) nitrate (ρNO3) uptake rates. In the time-series profiles, 2012–2013 data are not interpolated due to single-depth sampling. Grey shaded regions in depth-integrated plots indicate phytoplankton growing seasons for this study (March 1st – October 30th). Right panels show monthly-averaged depth profiles for the entire 8-year period, illustrating euphotic zone seasonality for carbon and nitrate uptake. The color scale bars on the far right apply to both the total time-series vertical profiles and the 8-year seasonal plots. Sampling depths are indicated by round symbols. The year labels are positioned under the tick marks corresponding to January.Full size imageSpiked seawater samples were incubated for 24 h, except from 2010 to 2013 when the incubation period was 4 to 6 hr. After incubation, the entire sample was gently vacuum filtered onto a combusted 0.7 µm (nominal porosity) glass fibre filter. Filters were dried for 48 h at 60 °C and kept in a desiccator at room temperature until analysis. Filters were packed into pellets and sent to the Stable Isotope Facility at the University of California (UC) Davis for analysis of 13C and 15N enrichment, and total C and N content by continuous flow isotope ratio mass spectrometry and elemental analysis, respectively. For these measurements, UC Davis uses either an Elementar Vario EL Cube or Micro Cube elemental analyzer (Elementar Analysensysteme GmbH, Hanau, Germany) interfaced to either a PDZ Europa 20–20 isotope ratio mass spectrometer (Sercon Ltd., Cheshire, UK) or an Isoprime VisION IRMS (Elementar UK Ltd, Cheadle, UK).Carbon and NO3− uptake rates were calculated using Eq. 3 of Hama et al.29, and Eq. 3 and 6 of Dugdale and Wilkerson30, respectively.For samples incubated for less than 24 h, the daily C or NO3− uptake rates (ρX) were calculated using a PAR extrapolation method shown in Eq. (2):$$rho X(mu mol,{L}^{-1}da{y}^{-1})=left(rho X(mu mol,{L}^{-1}h{r}^{-1})divleft(frac{PAR,during,incubation}{Total,Daily,PAR}right)right)ast 24$$
(2)
Additionally, to account for NO3− uptake occurring under no light, ρNO3− was measured in dark bottles and this rate was added to the ρNO3− of each sample incubated for less than 24 h. The ρNO3− DARK was calculated following Eq. (3):$$rho N{O}_{3}DARK(mu mol,{L}^{-1}da{y}^{-1})=left(rho N{O}_{3}DARK(mu mol,{L}^{-1}h{r}^{-1})divleft(frac{Total,Daily,PAR-PAR,during,incubation}{Total,Daily,PAR}right)right)ast 24$$
(3)
PAR data used in Eqs. (2) and (3) came from the LI-COR® LI-190 Quantum sensor that was mounted beside the incubator. The seawater DIC value for each sample was calculated using a regression equation relating water density to DIC for Saanich Inlet31. Ambient NO3− concentrations were measured as described above. More

Losos, E. & Leigh, E. G. Tropical Forest Diversity and Dynamism: Findings from a Large-Scale Plot Network (Univ. Chicago Press, 2004).Pan, Y. et al. A large and persistent carbon sink in the world’s forests. Science 333, 988–993 (2011).CAS 
PubMed 

Google Scholar 
Hansen, M. C. C. et al. High-resolution global maps of 21st-century forest cover change. Science 342, 850–854 (2013).CAS 
PubMed 

Google Scholar 
Chazdon, R. L. Beyond deforestation: restoring degraded lands. Science 1458, 1458–1460 (2008).
Google Scholar 
Global Forest Resources Assessment 2010 (FAO, 2010).Rozendaal, D. M. A. & Chazdon, R. L. Demographic drivers of tree biomass change during secondary succession in northeastern Costa Rica. Ecol. Appl. 25, 506–516 (2015).PubMed 

Google Scholar 
Poorter, L. et al. Biomass resilience of Neotropical secondary forests. Nature 530, 211–214 (2016).CAS 
PubMed 

Google Scholar 
Chazdon, R. L., Broadbent, E. N., Rozendaal, D. M. A., Bongers, F. & Al, E. Carbon sequestration potential of second-growth forest regeneration in the Latin American tropics. Sci. Adv. 2, e1501639 (2016).PubMed 
PubMed Central 

Google Scholar 
Lohbeck, M. et al. Functional diversity changes during tropical forest succession. Perspect. Plant Ecol. Evol. Syst. 14, 89–96 (2012).
Google Scholar 
Poorter, L. et al. Wet and dry tropical forests show opposite successional pathways in wood density but converge over time. Nat. Ecol. Evol. 3, 928–934 (2019).PubMed 

Google Scholar 
Townsend, A. R., Cleveland, C. C., Houlton, B. Z., Alden, C. B. & White, J. W. Multi-element regulation of the tropical forest carbon cycle. Front. Ecol. Environ. 9, 9–17 (2011).
Google Scholar 
Medvigy, D. et al. Observed variation in soil properties can drive large variation in modelled forest functioning and composition during tropical forest secondary succession. New Phytol. 223, 1820–1833 (2019).Powers, J. S., Mar, E. & Marín-Spiotta, E. Ecosystem processes and biogeochemical cycles during secondary tropical forest succession. Annu. Rev. Ecol. Evol. Syst. 48, 497–519 (2017).
Google Scholar 
Davidson, E. A. et al. Recuperation of nitrogen cycling in Amazonian forests following agricultural abandonment. Nature 447, 995–998 (2007).CAS 
PubMed 

Google Scholar 
Davidson, E. A. & Martinelli, L. A. in Amazonia and Global Change (eds Keller, M. et al.) 299–309 (American Geophysical Union, 2013).Vitousek, P. M. & Howarth, R. W. Nitrogen limitation on land and in the sea: how can it occur? Biogeochemistry 13, 87–115 (1991).
Google Scholar 
Batterman, S. A. et al. Key role of symbiotic dinitrogen fixation in tropical forest secondary succession. Nature 502, 224–227 (2013).CAS 
PubMed 

Google Scholar 
Bauters, M., Mapenzi, N., Kearsley, E., Vanlauwe, B. & Boeckx, P. Facultative nitrogen fixation by legumes in the central Congo basin is downregulated during late successional stages. Biotropica 48, 281–284 (2016).
Google Scholar 
Van Langenhove, L. et al. Regulation of nitrogen fixation from free-living organisms in soil and leaf litter of two tropical forests of the Guiana shield. Plant Soil 450, 93–110 (2020).PubMed 

Google Scholar 
Vitousek, P. M. Litterfall, nutrient cycling, and nutrient limitation in tropical forests. Ecology 65, 285–298 (1984).CAS 

Google Scholar 
Kaspari, M. et al. Multiple nutrients limit litterfall and decomposition in a tropical forest. Ecol. Lett. 11, 35–43 (2008).PubMed 

Google Scholar 
Cleveland, C. C. et al. Relationships among net primary productivity, nutrients and climate in tropical rain forest: a pan-tropical analysis. Ecol. Lett. 14, 939–947 (2011).PubMed 

Google Scholar 
Chadwick, O. A., Derry, L. A., Vitousek, P. M., Huebert, B. J. & Hedin, L. O. Changing sources of nutrients during four million years of ecosystem development. Nature 397, 491–497 (1999).CAS 

Google Scholar 
Hedin, L. O. et al. Nutrient losses over four million years of tropical forest development. Ecology 84, 2231–2255 (2003).
Google Scholar 
Sanchez, P. A., Villachica, J. H. & Bandy, D. E. Soil fertility dynamics after clearing a tropical rainforest in Peru. Soil Sci. Soc. Am. J. 47, 1171 (1983).CAS 

Google Scholar 
Davidson, E. A. et al. Nitrogen and phosphorus limitation of biomass growth in a tropical secondary forest. Ecol. Appl. 14, 150–163 (2004).
Google Scholar 
Wardle, D. A., Walker, L. R. & Bardgett, R. D. Ecosystem properties and forest decline in contrasting long-term chronosequences. Science 305, 509–513 (2004).CAS 
PubMed 

Google Scholar 
Wassen, M. J., Venterink, H. O., Lapshina, E. D. & Tanneberger, F. Endangered plants persist under phosphorus limitation. Nature 437, 547–550 (2005).CAS 
PubMed 

Google Scholar 
Waring, B. G., Becknell, J. M. & Powers, J. S. Nitrogen, phosphorus, and cation use efficiency in stands of regenerating tropical dry forest. Oecologia 178, 887–897 (2015).PubMed 

Google Scholar 
De longe, M., D’odorico, P. & Lawrence, D. Feedbacks between phosphorus deposition and canopy cover: the emergence of multiple stable states in tropical dry forests. Glob. Change Biol. 14, 154–160 (2008).
Google Scholar 
Bauters, M. et al. Fire-derived phosphorus fertilization of African Tropical Forests. Nat. Commun. 12, 5129 (2021).Vitousek, P. M. & Reiners, W. A. Ecosystem succession and nutrient retention: a hypothesis. Bioscience 25, 376–381 (1975).CAS 

Google Scholar 
Gallarotti, N. et al. In-depth analysis of N 2O fluxes in tropical forest soils of the Congo Basin combining isotope and functional gene analysis. ISME J. 15, 3357–3374 (2021).Gorham, E., Vitousek, P. M. & Reiners, W. A. The regulation of chemical budgets over the course of terrestrial ecosystem succession. Annu. Rev. Ecol. Syst. 10, 53–84 (1979).CAS 

Google Scholar 
Markewitz, D., Davidson, E., Moutinho, P. & Nepstad, D. Nutrient loss and redistribution after forest clearing on a highly weathered soil in Amazonia. Ecol. Appl. 14, 177–199 (2004).
Google Scholar 
Lawrence, D. et al. Ecological feedbacks following deforestation create the potential for a catastrophic ecosystem shift in tropical dry forest. Proc. Natl Acad. Sci. USA 104, 20696–20701 (2007).CAS 
PubMed 
PubMed Central 

Google Scholar 
Veldkamp, E., Schmidt, M., Powers, J. S. & Corre, M. D. Deforestation and reforestation impacts on soils in the tropics. Nat. Rev. Earth Environ. 1, 590–605 (2020).
Google Scholar 
Sanchez, P. A. Properties and Management of Soils in the Tropics (John Wiley and Sons, 1976).Turner, B. L. & Engelbrecht, B. M. J. Soil organic phosphorus in lowland tropical rain forests. Biogeochemistry 103, 297–315 (2011).Sullivan, B. W. et al. Biogeochemical recuperation of lowland tropical forest during succession. Ecology 100, e02641 (2019).Sardans, J. et al. Empirical support for the biogeochemical niche hypothesis in forest trees. Nat. Ecol. Evol. 13, 184–194 (2021).White, P. J. & Broadley, M. R. Calcium in plants. Ann. Bot. 92, 487–511 (2003).CAS 
PubMed 
PubMed Central 

Google Scholar 
Vitousek, P. M., Porder, S., Houlton, B. Z. & Chadwick, O. A. Terrestrial phosphorus limitation: mechanisms, implications, and nitrogen–phosphorus interactions. Ecol. Appl. 20, 5–15 (2010).PubMed 

Google Scholar 
Huggett, B. A., Schaberg, P. G., Hawley, G. J. & Eagar, C. Long-term calcium addition increases growth release, wound closure, and health of sugar maple (Acer saccharum) trees at the Hubbard Brook Experimental Forest. Can. J. For. Res. 37, 1692–1700 (2007).CAS 

Google Scholar 
Marschner, P. Marschner’s Mineral Nutrition of Higher Plants 3rd edn (Elsevier/Academic Press 2002).Walker, L. R., Wardle, D. A., Bardgett, R. D. & Clarkson, B. D. The use of chronosequences in studies of ecological succession and soil development. J. Ecol. 98, 725–736 (2010).
Google Scholar 
Bauters, M. et al. Soil nutrient depletion and tree functional composition shift following repeated clearing in secondary forests of the Congo Basin. Ecosystems 24, 1422–1435 (2021).Turner, B. L., Brenes-arguedas, T. & Condit, R. Pervasive phosphorus limitation of tree species but not communities in tropical forests. Nature 555, 367–370 (2018).CAS 
PubMed 

Google Scholar 
Wright, S. J. Plant responses to nutrient addition experiments conducted in tropical forests. Ecol. Monogr. 89, e01382 (2019).Lugli, L. F. et al. Rapid responses of root traits and productivity to phosphorus and cation additions in a tropical lowland forest in Amazonia. New Phytol. 230, 116–128 (2021).Vitousek, P. M. M. & Sanford, R. L. Nutrient cycling in moist tropical forest. Annu. Rev. Ecol. Syst. 17, 137–167 (1986).
Google Scholar 
Kaspari, M. & Powers, J. S. Biogeochemistry and geographical ecology: embracing all twenty-five elements required to build organisms. Am. Nat. 188, S62–S73 (2016).PubMed 

Google Scholar 
Nykvist, N. in Soils of Tropical Forest Ecosystems (eds Schulte, A. & Ruhiyat, D.) 87–91 (Springer, 1998).Bunyavejchewin, S., Sinbumroong, A., Turner, B. L. & Davies, S. J. Natural disturbance and soils drive diversity and dynamics of seasonal dipterocarp forest in Southern Thailand. J. Trop. Ecol. 35, 95–107 (2019).
Google Scholar 
Quesada, C. A. et al. Variations in chemical and physical properties of Amazon forest soils in relation to their genesis. Biogeosciences 7, 1515–1541 (2010).CAS 

Google Scholar 
Gerland, P. et al. World population stabilization unlikely this century. Science 346, 234–237 (2014).Makelele, I. A. et al. Afrotropical secondary forests exhibit fast diversity and functional recovery, but slow compositional and carbon recovery after shifting cultivation. J. Veg. Sci. 32, e13071 (2021).Van Langenhove, L. et al. Atmospheric deposition of elements and its relevance for nutrient budgets of tropical forests. Biogeochemistry 149, 175–193 (2020).
Google Scholar 
Staelens, J. et al. Calculating dry deposition and canopy exchange with the canopy budget model: review of assumptions and application to two deciduous forests. Water Air Soil Pollut. 191, 149–169 (2008).CAS 

Google Scholar 
Hofhansl, F. et al. Topography strongly affects atmospheric deposition and canopy exchange processes in different types of wet lowland rainforest, southwest Costa Rica. Biogeochemistry 106, 371–396 (2011).
Google Scholar 
Schrijver, A. De, Nachtergale, L. & Staelens, J. Comparison of throughfall and soil solution chemistry between a high-density Corsican pine stand and a naturally regenerated silver birch stand. Environ Pollut. 131, 93–105 (2004).Eriksson, E. & Khunakasem, V. Chloride concentration in groundwater, recharge rate and rate of deposition of chloride in the Israel coastal plain. J. Hydrol. 7, 178–197 (1969).
Google Scholar 
Malhi, Y. et al. An international network to monitor the structure, composition and dynamics of Amazonian forests (RAINFOR). J. Veg. Sci. 13, 439 (2002).
Google Scholar 
Réjou-Méchain, M., Tanguy, A., Piponiot, C., Chave, J. & Hérault, B. biomass: an R package for estimating above-ground biomass and its uncertainty in tropical forests. Methods Ecol. Evol. 8, 1163–1167 (2017).
Google Scholar 
Chave, J. et al. Improved allometric models to estimate the aboveground biomass of tropical trees. Glob. Change Biol. 20, 3177–3190 (2014).
Google Scholar 
Malhi, Y. et al. The Global Ecosystems Monitoring network: monitoring ecosystem productivity and carbon cycling across the tropics. Biol. Conserv. 253, 108889 (2021).D’Angelo, E., Crutchfield, J. & Vandiviere, M. Rapid, sensitive, microscale determination of phosphate in water and soil. J. Environ. Qual. 30, 2206–2209 (2001).Rowland, A. P. & Haygarth, P. M. Determination of total dissolved phosphorus in soil solutions. J. Environ. Qual. 26, 410–415 (1997).CAS 

Google Scholar 
Vance, E. D., Brookes, P. C. & Jenkinson, D. S. An extraction method for measuring soil microbial biomass C. Soil Biol. Biochem. 19, 703–707 (1987).CAS 

Google Scholar 
Brookes, P. C., Powlson, D. S. & Jenkinson, D. S. Measurement of microbial biomass phosphorus in soil. Soil Biol. Biochem. 14, 319–329 (1982).CAS 

Google Scholar 
Kaiser, C. et al. Belowground carbon allocation by trees drives seasonal patterns of extracellular enzyme activities by altering microbial community composition in a beech forest soil. New Phytol. 187, 843–858 (2010).CAS 
PubMed 
PubMed Central 

Google Scholar 
Pérez-Harguindeguy, N. et al. New handbook for standardised measurement of plant functional traits worldwide. Aust. J. Bot. 61, 167–234 (2013).Poorter, L. et al. Multidimensional tropical forest recovery. Science 374, 1370–1376 (2021). More

Gabor transformThe Gabor transform has rarely been used as a feature in a landscape classification OBIA approach but has been used in other OBIA processes such as fingerprint enhancement and human iris detection and for data dimensionality reduction24,29,30,31,32,33,34,35. Gabor filters are a bandpass filter applied to an image to identify texture. The different Gabor bandpass filters mathematically model the visual cortical cells of mammalian brains and thus is expected to improve segmentation and classification accuracy when compared to a human delineated and classified image26,27.Samiappan et al.36 compared Gabor filters to other texture features (grey-level co-occurrence matrix, segmentation-based fractal texture analysis, and wavelet texture analysis) within the GEOBIA process, of a wetland, using sub-meter resolution multispectral imagery. These Gabor filters performed comparably, in overall classification accuracy and Kappa coefficients, with other texture features. However, they were still outperformed by all other texture features. This study did not use any other data for analysis for determining the performance of Gabor filters when paired with data sources such as spectral, NDVI, or LiDAR36,37. Wang et al.38 paired a Gabor transformation with a fast Fourier transformation for edge detection on an urban landscape image that contained uniform textures with promising results. Su30 used the textural attributes derived from Gabor filters for classification but had similar results to Samiappan et al.36 where they found that Gabor features were one of the least useful/influential that contributed to the classification of a mostly agricultural landscape.Gabor filters are a Fourier influenced wavelet transformation, or bandpass filter, that identifies texture as intervals in a 2-D Gaussian modulated sinusoidal wave. This modulation differentiates the Gabor transform from the Fourier transform23,26. These Gabor transformed wavelets are parameterized by the angle at which they alter the image and the frequency of the wavelet. Rather than smoothing an image at the cost of losing detail through Fourier transforms or median filters, Gabor transformed images identify the repeated pattern of localized pixels and gives them similar values if they are a part of the same repeated sequence. Gabor features can closely emulate the visual cortex of mammalian brains that utilize texture to identify objects26,27. This is based on the evaluation of neurons associated with the cortical vertex that respond to different images or light profiles39. Marcelja27 identified that cortical cells responded to signals that are localized frequencies of light like what is represented by the Gabor transformations. Within the frequency domain, the Gabor transform can be defined by Eq. (1):$$Gleft(u, v;f, theta right)= {e}^{-frac{{pi }^{2}}{{f}^{2}} ({gamma }^{2}({u}^{{prime}}-f{)}^{2}+{n}^{2}{v}^{{{prime}}2})}$$
(1)

where (f) is the user-determined frequency (or wavelength); (theta) is the user-determined orientation at which the wavelet is applied to the image; (gamma) and (n) are the standard deviations of the Gaussian function in either direction23,38. These parameters define the shape of the band pass filter and determines its effect on one-dimensional signals. Daugman26, created a 2-D application of this filter in Eq. (2);$$gleft(u,vright)= {e}^{-{pi }^{2}/{f}^{2}[{gamma }^{2}{left({u}^{{prime}}-fright)}^{2}+{n}^{2}{{v}^{{prime}}}^{2}]}$$
(2)

where u’ = ucos − vsin θ θ and v’ = usin − vcos θ.In order to implement Gabor filters on multi-band spectral images, we used Matlab’s Gabor feature on the University of Iowa’s Neon high performance computer (HPC)40 which has up to 512 GB of RAM, which was necessary for processing these images. The first implementation of Gabor filters was performed on a 1610 × 687 single band pixel array (a small subset of the study area), a filter bank of 4 orientations and 8 wavelengths, on a 32 GB RAM computer, and took approximately 8 h to complete. Filter banks are a set of Gabor filters with different parameters that is applied to the spectral image and are required to identify different textures with different orientations and frequencies. By lowering the number of wavelengths from 8 to 4 on an 8128 × 8128 single band pixel array on the same machine 32 GB RAM, the processing was reduced to an hour. Using the HPC, this was further reduced to approximately 90 s using the same filter bank. Before implementing on the HPC, the original spectral image was divided into manageable subsets with overlap in order to prevent ‘edge-effect.’ These images were converted to greyscale by averaging values across all three bands33. When wavelengths become too long, they no longer attribute the textural information desired from the image and therefore add unnecessary computing time. The wavelengths that were used for the filter bank were selected as increasing powers of two starting from 2.82842712475 ((24/sqrt{2})) up to the pixel length of the hypotenuse of the input image. From this, we used only 2.82842712475, 7.0710678, 17.6776695, and 44.19417382. The directional orientation was selected as 45° intervals, from 0 to 180: 0, 45, 90, 135. These parameters were based on the reasoning outlined within Jain and Farrokhina25. More directional orientations could have been included but four were used for computational efficiency. The radial frequencies were selected so that they could capture the different texture in the landscape represented by consistent changes in pixels values within each landcover class. When frequencies are too wide or fine of a width they no longer represent the textures of the different landcover classes and thus are not included. This selection of filter bank parameters are similar or the same as other studies that look into the use of Gabor features for OBIA25,30,31.From the different combinations of parameters (four directions and four frequencies) in the Gabor Transform filter bank, sixteen magnitude response images were created from the converted greyscale three band average image. To limit high local variance within the output Gabor texture images, a Gaussian filter was applied. The magnitude response values were normalized across the 16 different bands so that a Principal Component Analysis (PCA) could be applied. The first principal component of the PCA, from these Gabor transformed images, was used for this study since it limits the computation time to process 16 separate Gabor features, in addition to the other data sources, while still retaining the most amount of information from the different Gabor response features. The Gabor band that was used for this study can be viewed in Fig. 2.Figure 2Gabor transformation. Gabor transformed image of study area derived from original image using the first principal component of all gabor outputs using the filter bank parameters. Software: ArcMap (10.x).Full size imageSegmentationFor this study, we used the watershed algorithm for the segmentation of GEOBIA, implemented by ENVI version 5.0 Feature Extraction tool, due to its ubiquitous use within GEOBIA, its ability to create a hierarchy of segmented objects, and support within the literature as a reliable algorithm37,41,39,43. The watershed algorithm can either use a gradient image or intensity image for segmentation. Based on the observed results, this study used the intensity method. The intensity method averages the value of pixels across bands. Scale, a user-defined parameter, is selected to identify the threshold that decides if a given intensity value within the gradient image can be a boundary. This allows the user to decide the size of the objects created. A secondary, user-defined, parameter defines how similar, adjacent, objects need to be before they are combined or merged. The user arbitrarily selects the parameter value based on how it reduces both under and over segmentation. The parameters selected for this study were visually chosen based on a compromise between over and under segmentation relative to the hand demarcated objects.The merging of two separate objects was based on the full lambda schedule where the user selects a merging threshold ({t}_{i, j}) which is defined by Eq. (3):$${t}_{i, j}= frac{frac{left|{O}_{i}right|cdot left|{O}_{j}right|}{left|{O}_{i}right|+ left|{O}_{j}right|}cdot {Vert {u}_{i}-{u}_{j}Vert }^{2}}{mathrm{length}(mathrm{vartheta }left({O}_{i},{O}_{j}right))}$$
(3)

where ({O}_{i}) is the object of the image, (left|{O}_{i}right|) is the area of (i), ({u}_{i}) is the average of object (i), ({u}_{j}) is the average of object (j), (Vert {u}_{i}-{u}_{j}Vert) is the Euclidean distance between the average values of the pixel values in regions (i) and (j), and (mathrm{length}left(mathrm{vartheta }left({O}_{i},{O}_{j}right)right)) is the length of the shared boundary of ({O}_{i}) and ({O}_{j}).To compare the segmentation of a riparian landscape, with and without Gabor features, we conducted segmentation on two separate sets of data. One dataset was a normalized stacked layer of NDVI and CHM (see Fig. 3) with the original multispectral image used as ancillary data; the other dataset differed only by the inclusion of the Gabor feature. For both instances, the bands were converted to an intensity image by averaging across bands rather than being converted into a gradient image for segmentation. The dataset that included the Gabor features had a scale parameter set at 30 with merge settings at 95 and 95.7 for the sub and super-objects, respectively. The dataset that did not include the Gabor features had a scale parameter of 10 with merge settings at 95.6 and 98.5 for the sub and super-objects, respectively. This resulted in the creation of 87,198 and 62,905 segments for the sub and super objects, respectively, that were created when the Gabor feature was included. 191,050 and 51,664 segments were created for the sub and super objects when the Gabor features, respectively, were not included within the segmentation process. As you will see in the next section, these segments also represent the number of training data that will be included within the supervised classification.Figure 3CHM and NDVI. LiDAR derived canopy height model (top) and normalized difference vegetation index derived from original spectral image. Software: ArcMap (10.x).Full size imageTo create a hierarchy of land cover classes, two sets of segmentation parameters needed to be selected for each dataset. One set of parameters would be used for the sub-objects within the hierarchy and the other set would be used to create super-objects. All parameters used the intensity and full lambda schedule algorithms for the watershed method. The only setting that changed between the sub and super-objects, for either dataset, was the merge parameter which helped maintain similar boundaries as much as possible. Despite this, boundaries could moderately change due to the Euclidean distance, between the pixel values of (i) and (j), changing from the merging of objects; causing ({t}_{i, j}) to cross the threshold which results in a new boundary being drawn. A representation of these results can be viewed and visually compared to the hand demarcated objects in Fig. 4.Figure 4Automated and manual segmented comparison. Juxtaposition of hand delineated, sub-objects, and super-objects for segments generated using the Gabor features. Software: ArcMap (10.x).Full size imageTraining dataThe training data, used for this study, is the transfer of class attributes from hand demarcated and classified segments to automatically segmented objects based on the majority overlap of the hand demarcated segments. Experts identified them using two different classification schemes referenced from the General Wetland Vegetation Classification System44. The 7-class scheme within this system identified objects of either being forest, marsh, agriculture, developed, open water, grass/forbs, or sand/mud. The 13-class scheme identified objects of either being agriculture, developed, grass/forbs, open water, road/levee, sand/mud, scrub-shrub, shallow marsh, submerged aquatic vegetation, upland forest, wet forest, wet meadow, and wet shrub. Not every class from the 7-class scheme will have a sub-class (i.e. developed, open water) but some do for example wet and upland forest are sub-objects of the forest class and wet meadow and shallow marsh are sub-objects of marsh. Figure 5 visually illustrates both classification schemes across the study area.Figure 5Hand delineated objects of both scales. Software: ArcMap (10.x).Full size imageENVI’s feature extraction tool calculates several landscape, spectral, and textural metrics. These attributes were used for each random forest classifier. The Gabor and Hierarchical features will be included selectively to be able to compare their contributions to the (out-of-bag) OOB classification errors. When Gabor features are included within the classification, they are computed the same way as the other image bands.Random forestThe random forest classifier was implemented in R using the random forest module45. The number of trees, that were randomly generated, was large enough (n = 250) to where the Strong law of large numbers would take effect as indicated by the decrease in the change of accuracy. The default number of variables randomly sampled as candidates at each split variable (mtry parameter) was the total number of variables divided by 3 for each dataset. R also generates two separate variable indices: mean decrease in accuracy and mean decrease Gini. Mean decrease in accuracy refers to the accuracy change in the random forest when a single variable is left out. This is a practical metric to determine the usefulness of a variable. The Gini index measures the purity change within a dataset when it is split based upon a given variable within a decision tree.The random forest classification accuracy will be based on the OOB error. The random forest algorithm trains numerous decision trees on random subsets of the training set leaving out a number of training samples when training each decision tree. The samples that are left out of each decision tree are then classified by the decision tree that they were not included within during the training step. The OOB error is the average error of each predicted bootstrapped sample across the ensemble of decision trees within the random forest algorithm.Figure 6 illustrates how the Gabor and hierarchal features were included within the classification of the super and sub-objects.Figure 6Classification procedure. Schematic flow chart illustrating how the Gabor and hierarchal features were included within the classification of the super and sub-objects. OOB classification error included in parenthesis.Full size imageHierarchical schemeTo attribute the hierarchical structure to the sub-objects, we first classified the larger segments that were created with and without the Gabor features using the broader 7-class scheme. These classified super objects were then converted to raster to calculate the majority overlap with the smaller sub-objects. This gave the sub-objects an attribute, the broader 7-class scheme, that could be used to contribute to the classification of the sub-objects with the finer 13-class scheme. This builds the hierarchical relationship between the two class schemes into the supervised classification of the sub-objects. Figure 6 illustrates how the hierarchal structure was included within two of the four sub-object’s list of features used within classification. This methodological approach aligns with O’Neill et al.21 landscape ecology principle that a super-object’s class could be a useful property in defining or predicting a sub-object. This is also different than the more common rule-based approach of iteratively classifying the landscape into smaller and smaller sub-classes22.Segmentation assessmentMost studies rely upon the accuracy assessment of their classifiers to provide support for their analysis results. However, this does not provide evidence whether a new data fusion technique improves the ability to delineate objects of interest within an image. To assess the performance of our segmented polygons, this study evaluated the segments created with and without the Gabor feature using a method highlighted in Xiao et al.37.Our segmentation results were evaluated using an empirical discrepancy measure, used frequently in image segmentation evaluation37,46,47. Discrepancy measures utilize ground truth images that represent the “correct” delineated/classified image to compare the semi-automated image results. In our study, the objects that were delineated and classified by experts from the U.S. Fish and Wildlife Service, were used as training data for our random forest classifier and as ground truth for the discrepancy measure. The discrepancy measure used the percentage of right segmented pixels (PR) in the whole image. To calculate PR, we converted the classified segmented and ground truth polygons to raster and measured the ratio of incorrect pixels to total amount of pixels which was converted to a percentage.Additionally, landscape metrics were calculated using FRAGSTATS48, an open source program commonly used for calculating landscape metrics. FRAGSTATS computed these metrics from thematic raster maps that represent the land cover types of interest. These thematic classes, used for analysis, were the classified objects at both the super and sub-object level. Since we are not attempting to compare the segmentation results for any specific class or area, we calculated metrics on a landscape level. Landscape metrics will represent the segmentation patterns for the entire study area.FRAGSTATS can calculate various metrics representing different aspects of the landscape. The metrics for analysis attempts to understand object geometry. The metrics calculated, for these analyses, were the average and standard deviation for the area (AREA), the fractal dimension index (FRAC), and the perimeter area ratio (PARA). The number of patches (NP) was also included in each result. To take a more landscape centric approach, the area weighted mean was chosen over a simple average. More

Huxley, J. Evolution. The Modern Synthesis (Allen & Unwin, 1942).Bay, R. A. & Palumbi, S. R. Rapid acclimation ability mediated by transcriptome changes in reef-building corals. Genome Biol. Evol. 7, 1602–1612 (2015).CAS 
PubMed 
PubMed Central 

Google Scholar 
Palumbi, S. R., Barshis, D. J., Traylor-Knowles, N. & Bay, R. A. Mechanisms of reef coral resistance to future climate change. Science 344, 895–898 (2014).CAS 
PubMed 

Google Scholar 
Bang, C. et al. Metaorganisms in extreme environments: do microbes play a role in organismal adaptation? Zoology 127, 1–19 (2018).PubMed 

Google Scholar 
Fraune, S., Forêt, S. & Reitzel, A. M. Using Nematostella vectensis to study the interactions between genome, epigenome, and bacteria in a changing environment. Front. Mar. Sci. 3, 1–8 (2016).
Google Scholar 
Kolodny, O. & Schulenburg, H. Opinion piece Microbiome-mediated plasticity directs host evolution along several distinct time scales. Phil. Trans. R. Soc. B 375, 20190589 (2020).Reshef, L., Koren, O., Loya, Y., Zilber-Rosenberg, I. & Rosenberg, E. The coral probiotic hypothesis. Environ. Microbiol. 8, 2068–2073 (2006).CAS 
PubMed 

Google Scholar 
Webster, N. S. & Reusch, T. B. H. Microbial contributions to the persistence of coral reefs. ISME J. 11, 2167–2174 (2017).PubMed 
PubMed Central 

Google Scholar 
Totton, A. K. The British sea anemones. Nature 135, 977–978 (1935).
Google Scholar 
Hand, C. & Uhlinger, K. R. The unique, widely distributed, estuarine sea anemone, Nematostella vectensis Stephenson: a review, new facts, and questions. Estuaries 17, 501–501 (1994).
Google Scholar 
Darling, J. A., Reitzel, A. M. & Finnerty, J. R. Regional population structure of a widely introduced estuarine invertebrate: Nematostella vectensis Stephenson in New England. Mol. Ecol. 13, 2969–2981 (2004).CAS 
PubMed 

Google Scholar 
Darling, J. A. et al. Rising starlet: the starlet sea anemone, Nematostella vectensis. BioEssays 27, 211–221 (2005).CAS 
PubMed 

Google Scholar 
Hand, C. & Uhlinger, K. R. The culture, sexual and asexual reproduction, and growth of the sea anemone Nematostella vectensis. Biol. Bull. 182, 169–176 (1992).CAS 
PubMed 

Google Scholar 
Pearson, C. V. M., Rogers, A. D. & Sheader, M. The genetic structure of the rare lagoonal sea anemone, Nematostella vectensis Stephenson (Cnidaria; Anthozoa) in the United Kingdom based on RAPD analysis. Mol. Ecol. 11, 2285–2293 (2002).CAS 
PubMed 

Google Scholar 
Reitzel, A. M., Darling, J. A., Sullivan, J. C. & Finnerty, J. R. Global population genetic structure of the starlet anemone Nematostella vectensis: multiple introductions and implications for conservation policy. Biol. Invasions 10, 1197–1213 (2008).
Google Scholar 
Stefanik, D. J., Friedman, L. E. & Finnerty, J. R. Collecting, rearing, spawning and inducing regeneration of the starlet sea anemone, Nematostella vectensis. Nat. Protoc. 8, 916–923 (2013).PubMed 

Google Scholar 
Fritzenwanker, J. H. & Technau, U. Induction of gametogenesis in the basal cnidarian Nematostella vectensis (Anthozoa). Dev. Genes Evol. 212, 99–103 (2002).PubMed 

Google Scholar 
Mortzfeld, B. M. et al. Response of bacterial colonization in Nematostella vectensis to development, environment, and biogeography. Environ. Microbiol. 18, 1764–1781 (2016).PubMed 

Google Scholar 
Baldassarre, L. et al. Contribution of maternal and paternal transmission to bacterial colonization in Nematostella vectensis. Front. Microbiol. 12, 2892 (2021).
Google Scholar 
Domin, H. et al. Predicted bacterial interactions affect in vivo microbial colonization dynamics in Nematostella. Front. Microbiol. 9, 728 (2018).Guest, J. J. R. et al. Contrasting patterns of coral bleaching susceptibility in 2010 suggest an adaptive response to thermal stress. PLoS ONE 7, e33353–e33353 (2012).CAS 
PubMed 
PubMed Central 

Google Scholar 
Puisay, A., Pilon, R., Goiran, C. & Hédouin, L. Thermal resistances and acclimation potential during coral larval ontogeny in Acropora pulchra. Mar. Environ. Res. 135, 1–10 (2018).CAS 
PubMed 

Google Scholar 
Van Oppen, M. J. H., Oliver, J. K., Putnam, H. M. & Gates, R. D. Building coral reef resilience through assisted evolution. Proc. Natl Acad. Sci. USA 112, 2313 (2015).
Google Scholar 
Torda, G. et al. Rapid adaptive responses to climate change in corals. Nat. Clim. Change 7, 627–636 (2017).
Google Scholar 
Yu, Xiaopeng et al. Thermal acclimation increases heat tolerance of the scleractinian coral Acropora pruinosa,. Sci. Total Environ. 733, 139319–139319 (2020).CAS 
PubMed 

Google Scholar 
Jury, C. P. & Toonen, R. J. Adaptive responses and local stressor mitigation drive coral resilience in warmer, more acidic oceans. Proc. R. Soc. B Biol. Sci. 286, 20190614–20190614 (2019).
Google Scholar 
Sully, S., Burkepile, D. E., Donovan, M. K., Hodgson, G. & van Woesik, R. A global analysis of coral bleaching over the past two decades. Nat. Commun. 10, 5 (2019).
Google Scholar 
Thomas, L. et al. Mechanisms of thermal tolerance in reef-building corals across a fine-grained environmental mosaic: lessons from Ofu,. Am. Samoa. Front. Mar. Sci. 4, 434 (2018).
Google Scholar 
Oliver, T. A. & Palumbi, S. R. Many corals host thermally resistant symbionts in high-temperature habitat. Coral Reefs 30, 241–250 (2011).
Google Scholar 
Kenkel, C. D. & Matz, M. V. Gene expression plasticity as a mechanism of coral adaptation to a variable environment. Nat. Ecol. Evol. 1, 14 (2017).Barker, V. Exceptional thermal tolerance of coral reefs in American Samoa a review. Curr. Clim. Change Rep. 4, 427 (2018).
Google Scholar 
Bourne, D., Iida, Y., Uthicke, S. & Smith-Keune, C. Changes in coral-associated microbial communities during a bleaching event. ISME J. 2, 350–63 (2008).CAS 
PubMed 

Google Scholar 
Carrier, T. J. & Reitzel, A. M. The hologenome across environments and the implications of a host-associated microbial repertoire. Front. Microbiol. 8, 802 (2017).Koren, O. & Rosenberg, E. Bacteria associated with mucus and tissues of the coral Oculina patagonica in summer and winter. Appl. Environ. Microbiol. 72, 5254–5259 (2006).CAS 
PubMed 
PubMed Central 

Google Scholar 
Littman, R., Willis, B. L. & Bourne, D. G. Metagenomic analysis of the coral holobiont during a natural bleaching event on the Great Barrier Reef. Environ. Microbiol. Rep. 3, 651–60 (2011).CAS 
PubMed 

Google Scholar 
Ziegler, M., Seneca, F. O., Yum, L. K., Palumbi, S. R. & Voolstra, C. R. Bacterial community dynamics are linked to patterns of coral heat tolerance. Nat. Commun. 8, 14213–14213 (2017).CAS 
PubMed 
PubMed Central 

Google Scholar 
Thurber, R. V. et al. Metagenomic analysis of stressed coral holobionts. Environ. Microbiol. 11, 2148–2163 (2009).CAS 

Google Scholar 
van Oppen, M. J. H. & Blackall, L. L. Coral microbiome dynamics, functions and design in a changing world. Nat. Rev. Microbiol. 17, 557–567 (2019).PubMed 

Google Scholar 
Moran, N. A. & Yun, Y. Experimental replacement of an obligate insect symbiont. Proc. Natl Acad. Sci. USA 112, 2093–2096 (2015).CAS 
PubMed 
PubMed Central 

Google Scholar 
Ainsworth, T. D. T. et al. The coral core microbiome identifies rare bacterial taxa as ubiquitous endosymbionts. ISME J. 9, 2261–2274 (2015).CAS 

Google Scholar 
Hester, E. R., Barott, K. L., Nulton, J., Vermeij, M. J. A. & Rohwer, F. L. Stable and sporadic symbiotic communities of coral and algal holobionts. ISME J. 10, 1157–1169 (2016).CAS 
PubMed 

Google Scholar 
Bourne, D. G., Morrow, K. M. & Webster, N. S. Insights into the coral microbiome: underpinning the health and resilience of reef ecosystems. Annu. Rev. Microbiol. 70, 340 (2016).
Google Scholar 
Pollock, F. J. et al. Reduced diversity and stability of coral-associated bacterial communities and suppressed immune function precedes disease onset in corals. R. Soc. Open Sci. 6, 31312497 (2019).Zilber-Rosenberg, I. & Rosenberg, E. Role of microorganisms in the evolution of animals and plants: the hologenome theory of evolution. FEMS Microbiol. Rev. 32, 723–735 (2008).CAS 
PubMed 

Google Scholar 
Elena, S. F. & Lenski, R. E. Evolution experiments with microorganisms: the dynamics and genetic bases of adaptation. Nat. Rev. Genet. 4, 457–469 (2003).CAS 
PubMed 

Google Scholar 
Hehemann, J. H. et al. Transfer of carbohydrate-active enzymes from marine bacteria to Japanese gut microbiota. Nature 464, 908–912 (2010).CAS 
PubMed 

Google Scholar 
Bourne, D. G. Microbiological assessment of a disease outbreak on corals from Magnetic Island (Great Barrier Reef, Australia). Coral Reefs 24, 304–312 (2005).
Google Scholar 
Leach, W. B., Carrier, T. J. & Reitzel, A. M. Diel patterning in the bacterial community associated with the sea anemone Nematostella vectensis. Ecol. Evol. 9, 9935–9947 (2019).PubMed 
PubMed Central 

Google Scholar 
Pootakham, W. et al. Heat-induced shift in coral microbiome reveals several members of the Rhodobacteraceae family as indicator species for thermal stress in Porites lutea. MicrobiologyOpen 8, e935 (2019).Webster, N. Host-associated coral reef microbes respond to the cumulative pressures of ocean warming and ocean acidification. Sci. Rep. 6, 19324 (2016).Van, K. L., Ae, A., Schupp, P. & Slattery, M. The distribution of dimethylsulfoniopropionate in tropical Pacific coral reef invertebrates. Coral Reefs 25, 321–327 (2006).
Google Scholar 
Rypien, K. L., Ward, J. R. & Azam, F. Antagonistic interactions among coral-associated bacteria. Environ. Microbiol. 12, 28–39 (2010).CAS 
PubMed 

Google Scholar 
Blazejak, A., Erséus, C., Amann, R. & Dubilier, N. Coexistence of bacterial sulfide oxidizers, sulfate reducers, and spirochetes in a gutless worm (oligochaeta) from the Peru margin. Appl. Environ. Microbiol. 71, 1553–1561 (2005).CAS 
PubMed 
PubMed Central 

Google Scholar 
Dubilier, N. et al. Phylogenetic diversity of bacterial endosymbionts in the gutless marine oligochete Olavius loisae (Annelida). Mar. Ecol. Prog. Ser. 178, 271–280 (1999).
Google Scholar 
Rincón-Rosales, R., Lloret, L., Ponce, E. & Martínez-Romero, E. Erratum: Rhizobia with different symbiotic efficiencies nodulate Acaciella angustissima in Mexico, including Sinorhizobium chiapanecum sp. nov. which has common symbiotic genes with Sinorhizobium mexicanum (FEMS Microbiology Ecology (2009) 67 (103-117)). FEMS Microbiol. Ecol. 68, 255–255 (2009).
Google Scholar 
Rosenberg, E. & DeLong, E. F., Stackebrandt, E., Lory, S., Thompson, F. The Prokaryotes—Prokaryotic Biology and Symbiotic Associations. (Springer, 2013).Kimura, H., Higashide, Y. & Naganuma, T. Endosymbiotic microflora of the Vestimentiferan Tubeworm (Lamellibrachia sp.) from a Bathyal Cold Seep. Mar. Biotechnol. 5, 593–603 (2003).CAS 

Google Scholar 
Melillo, A. A., Bakshi, C. S. & Melendez, J. A. Francisella tularensis antioxidants harness reactive oxygen species to restrict macrophage signaling and cytokine production. J. Biol. Chem. 285, 27553–27560 (2010).CAS 
PubMed 
PubMed Central 

Google Scholar 
Rabadi, S. M. et al. Antioxidant defenses of Francisella tularensis modulate macrophage function and production of proinflammatory cytokines. J. Biol. Chem. 291, 5009–5021 (2016).CAS 
PubMed 

Google Scholar 
McBride, M. J. in The Prokaryotes: Other Major Lineages of Bacteria and The Archaea. Vol. 9783642389542, 643–676 (Springer-Verlag Berlin Heidelberg, 2014).Augustin, R., Fraune, S. & Bosch, T. C. G. How Hydra senses and destroys microbes. Semin. Immunol. 22, 54–58 (2010).CAS 
PubMed 

Google Scholar 
Augustin, R. et al. A secreted antibacterial neuropeptide shapes the microbiome of Hydra. Nat. Commun. 8, 698 (2017).Franzenburg, S. et al. Distinct antimicrobial peptide expression determines host species-specific bacterial associations. Proc. Natl Acad. Sci. USA 110, E3730–E3738 (2013).CAS 
PubMed 
PubMed Central 

Google Scholar 
Fraune, S., Abe, Y. & Bosch, T. C. G. G. Disturbing epithelial homeostasis in the metazoan Hydra leads to drastic changes in associated microbiota. Environ. Microbiol. 11, 2361–9 (2009).CAS 
PubMed 

Google Scholar 
Brennan, J. J. et al. Sea anemone model has a single Toll-like receptor that can function in pathogen detection, NF-κB signal transduction, and development. Proc. Natl Acad. Sci. USA 114, E10122–E10131 (2017).CAS 
PubMed 
PubMed Central 

Google Scholar 
Sullivan, J. C. et al. Two alleles of NF-κB in the sea anemone Nematostella vectensis are widely dispersed in nature and encode proteins with distinct activities. PLoS ONE 4, e7311 (2009).Wolenski, F. S. et al. Characterization of the core elements of the NF-B signaling pathway of the sea anemone Nematostella vectensis. Mol. Cell. Biol. 31, 1076–1087 (2011).CAS 
PubMed 

Google Scholar 
Gáliková, M., Klepsatel, P., Senti, G. & Flatt, T. Steroid hormone regulation of C. elegans and Drosophila aging and life history. Exp. Gerontol. 46, 141–147 (2011).PubMed 

Google Scholar 
Taubenheim, J., Kortmann, C. & Fraune, S. Function and evolution of nuclear receptors in environmental-dependent postembryonic development. Front. Cell Dev. Biol. 9, 653792 (2021).PubMed 
PubMed Central 

Google Scholar 
Becker, P. B. & Workman, J. L. Nucleosome remodeling and epigenetics. Cold Spring Harb. Perspect. Biol. 5, a017905–a017905 (2013).PubMed 
PubMed Central 

Google Scholar 
Barno, A. R., Villela, H. D. M., Aranda, M., Thomas, T. & Peixoto, R. S. Host under epigenetic control: a novel perspective on the interaction between microorganisms and corals. BioEssays 43, 2100068.Reitzel, A. M. et al. Physiological and developmental responses to temperature by the sea anemone Nematostella vectensis. Mar. Ecol. Prog. Ser. 484, 115–130 (2013).
Google Scholar 
Chua, C. M., Leggat, W., Moya, A. & Baird, A. H. Temperature affects the early life history stages of corals more than near future ocean acidification. Mar. Ecol. Prog. Ser. 475, 85–92 (2013).
Google Scholar 
Ericson, J. A. et al. Combined effects of two ocean change stressors, warming and acidification, on fertilization and early development of the Antarctic echinoid Sterechinus neumayeri. Polar Biol. 35, 1027–1034 (2012).
Google Scholar 
Sheppard Brennand, H., Soars, N., Dworjanyn, S. A., Davis, A. R. & Byrne, M. Impact of ocean warming and ocean acidification on larval development and calcification in the sea urchin Tripneustes gratilla. PLoS ONE 5, e11372 (2010).Bernal, M. A. et al. Phenotypic and molecular consequences of stepwise temperature increase across generations in a coral reef fish. Mol. Ecol. 27, 4516–4528 (2018).CAS 
PubMed 

Google Scholar 
Clark, M. S. et al. Molecular mechanisms underpinning transgenerational plasticity in the green sea urchin Psammechinus miliaris. Sci. Rep. 9, 1–12 (2019).
Google Scholar 
Donelson, J. et al. Rapid transgenerational acclimation of a tropical reef fish to climate change. Nat. Clim. Change 2, 30–32 (2012).
Google Scholar 
Miller, G. M., Watson, S. A., Donelson, J. M., McCormick, M. I. & Munday, P. L. Parental environment mediates impacts of increased carbon dioxide on a coral reef fish. Nat. Clim. Change 2, 858–861 (2012).CAS 

Google Scholar 
Munday, P. L. Transgenerational acclimation of fishes to climate change and ocean acidification. F1000Prime Rep. 6, 99–99 (2014).PubMed 
PubMed Central 

Google Scholar 
Ryu, T. et al. An epigenetic signature for within-generational plasticity of a reef fish to ocean warming. Front. Mar. Sci. 7, 284 (2020).Veilleux, H. et al. Molecular processes of transgenerational acclimation to a warming ocean. Nat. Clim. Change 5, 1074–1078 (2015).CAS 

Google Scholar 
Zhao, C. et al. Transgenerational effects of ocean warming on the sea urchin Strongylocentrotus intermedius. Ecotoxicol. Environ. Saf. 151, 212–219 (2018).CAS 
PubMed 

Google Scholar 
Eirin-Lopez, J. M. & Putnam, H. M. Marine Environmental Epigenetics. Annu. Rev. Mar. Sci. 11, 335–368 (2019).
Google Scholar 
Fallet, M., Luquet, E., David, P. & Cosseau, C. Epigenetic inheritance and intergenerational effects in mollusks. Gene 729, 144166–144166 (2020).CAS 
PubMed 

Google Scholar 
Putnam, H. M. & Gates, R. D. Preconditioning in the reef-building coral Pocillopora damicornis and the potential for trans-generational acclimatization in coral larvae under future climate change conditions. J. Exp. Biol. 218, 2365–2372 (2015).PubMed 

Google Scholar 
Daxinger, L. & Whitelaw, E. Transgenerational epigenetic inheritance: more questions than answers. Genome Res. 20, 1623–1628 (2010).CAS 
PubMed 
PubMed Central 

Google Scholar 
Ptashne, M. Epigenetics: core misconcept. Proc. Natl Acad. Sci. USA 110, 7101–7103 (2013).CAS 
PubMed 
PubMed Central 

Google Scholar 
Rivera, H. E., Chen, C.-Y., Gibson, M. C. & Tarrant, A. M. Plasticity in parental effects confers rapid larval thermal tolerance in the estuarine anemone Nematostella vectensis. J. Exp. Biol. 224, jeb236745 (2021).Hirose, E. & Fukuda, T. Vertical transmission of photosymbionts in the colonial ascidian Didemnum molle: The larval tunic prevents symbionts from attaching to the anterior part of larvae. Zool. Sci. 23, 669–674 (2006).
Google Scholar 
Padilla-Gamiño, J. L., Pochon, X., Bird, C., Concepcion, G. T. & Gates, R. D. From parent to gamete: vertical transmission of Symbiodinium (Dinophyceae) ITS2 sequence assemblages in the reef building coral Montipora capitata. PLoS ONE 7, e38440–e38440 (2012).PubMed 
PubMed Central 

Google Scholar 
Sharp, K. H., Eam, B., John Faulkner, D. & Haygood, M. G. Vertical transmission of diverse microbes in the tropical sponge Corticium sp. Appl. Environ. Microbiol. 73, 622–629 (2007).CAS 
PubMed 

Google Scholar 
Sipkema, D. et al. Similar sponge-associated bacteria can be acquired via both vertical and horizontal transmission. Environ. Microbiol. 17, 3807–3821 (2015).CAS 
PubMed 

Google Scholar 
Apprill, A., Marlow, H. Q., Martindale, M. Q. & Rappé, M. S. The onset of microbial associations in the coral Pocillopora meandrina. ISME J. 3, 685–699 (2009).PubMed 

Google Scholar 
Sharp, K. H., Distel, D. & Paul, V. J. Diversity and dynamics of bacterial communities in early life stages of the Caribbean coral Porites astreoides. ISME J. 6, 790–801 (2012).CAS 
PubMed 

Google Scholar 
Lesser, M. P., Stat, M. & Gates, R. D. The endosymbiotic dinoflagellates (Symbiodinium sp.) of corals are parasites and mutualists. Coral Reefs 32, 603–611 (2013).
Google Scholar 
Ceh, J., Raina, J. B., Soo, R. M., van Keulen, M. & Bourne, D. G. Coral-bacterial communities before and after a coral mass spawning event on Ningaloo Reef. PLoS ONE 7, e36920 (2012).Ricardo, G. F., Jones, R. J., Negri, A. P. & Stocker, R. That sinking feeling: suspended sediments can prevent the ascent of coral egg bundles. Sci. Rep. 6, 21567 (2016).Leite, D. C. A. D. et al. Broadcast spawning coral Mussismilia Hispida can vertically transfer its associated bacterial core. Front. Microbiol. 8, 176–176 (2017).PubMed 
PubMed Central 

Google Scholar 
Epstein, H. E. et al. Microbiome engineering: enhancing climate resilience in corals. Front. Ecol. Environ. 17, 108 (2019).
Google Scholar 
Peixoto, R. S. et al. Beneficial microorganisms for corals (BMC) Proposed mechanisms for coral health and resilience. Front. Microbiol. 8, 341 (2017).PubMed 
PubMed Central 

Google Scholar 
Chakravarti, L. J., Beltran, V. H. & van Oppen, M. J. H. Rapid thermal adaptation in photosymbionts of reef-building corals. Glob. Change Biol. 23, 4675–4688 (2017).
Google Scholar 
Damjanovic, K., Blackall, L. L., Webster, N. S. & van Oppen, M. J. H. H. The contribution of microbial biotechnology to mitigating coral reef degradation. Microb. Biotechnol. 10, 1236–1243 (2017).PubMed 
PubMed Central 

Google Scholar 
Damjanovic, K., Van Oppen, M. J. H., Menéndez, P. & Blackall, L. L. Experimental inoculation of coral recruits with marine bacteria indicates scope for microbiome manipulation in Acropora tenuis and Platygyra daedalea. Front. Microbiol. 10, 1702 (2019).Rosado, P. M. et al. Marine probiotics: increasing coral resistance to bleaching through microbiome manipulation. ISME J. 13, 921–936 (2019).CAS 
PubMed 

Google Scholar 
Fraune, S. et al. Bacteria-bacteria interactions within the microbiota of the ancestral metazoan Hydra contribute to fungal resistance. ISME J. 9, 1543–1556 (2015).CAS 
PubMed 

Google Scholar 
Fadrosh, D. W. et al. An improved dual-indexing approach for multiplexed 16 S rRNA gene sequencing on the Illumina MiSeq platform. Microbiome 2, 6 (2014).PubMed 
PubMed Central 

Google Scholar 
Rausch, P. et al. Analysis of factors contributing to variation in the C57BL/6 J fecal microbiota across German animal facilities. Int. J. Med. Microbiol. 306, 343–355 (2016).PubMed 

Google Scholar 
Caporaso, J. G. et al. QIIME allows analysis of high-throughput community sequencing data. Nat. Methods 7, 335–336 (2010).CAS 
PubMed 
PubMed Central 

Google Scholar 
Faith, J. J. et al. The long-term stability of the human gut microbiota. Science 341, 1237439–1237439 (2013).PubMed 
PubMed Central 

Google Scholar 
Segata, N. et al. Metagenomic biomarker discovery and explanation. Genome Biol. 12, R60–R60 (2011).PubMed 
PubMed Central 

Google Scholar 
Andrews, S. FastQC: A Quality Control Tool for High Throughput Sequence Data [Online]. Available online at: http://www.bioinformatics.babraham.ac.uk/projects/fastqc/ (2010).Bolger, A. M., Lohse, M. & Usadel, B. Trimmomatic: a flexible trimmer for Illumina sequence data. Bioinformatics 30, 2114–2120 (2014).CAS 
PubMed 
PubMed Central 

Google Scholar 
Kim, D., Langmead, B. & Salzberg, S. L. HISAT: a fast spliced aligner with low memory requirements. Nat. Methods 12, 357–360 (2015).CAS 
PubMed 
PubMed Central 

Google Scholar 
Pertea, M. et al. StringTie enables improved reconstruction of a transcriptome from RNA-seq reads. Nat. Biotechnol. 33, 290–295 (2015).CAS 
PubMed 
PubMed Central 

Google Scholar 
Shao, M. & Kingsford, C. accurate assembly of transcripts through phase-preserving graph decomposition. Nat. Biotechnol. 35, 1167–1169 (2017).CAS 
PubMed 
PubMed Central 

Google Scholar 
Niknafs, Y. S., Pandian, B., Iyer, H. K., Chinnaiyan, A. M. & Iyer, M. K. TACO produces robust multisample transcriptome assemblies from RNA-seq. Nat. Methods 14, 68–70 (2016).PubMed 
PubMed Central 

Google Scholar 
Pertea, M. & Pertea, G. GFF Utilities: GffRead and GffCompare. F1000Research 9, 304–304 (2020).
Google Scholar 
Manni, M., Berkeley, M. R., Seppey, M., Simão, F. A. & Zdobnov, E. M. BUSCO update: novel and streamlined workflows along with broader and deeper phylogenetic coverage for scoring of eukaryotic, prokaryotic, and viral genomes. Mol. Biol. Evol. 38, 4647–4654 (2021).PubMed 
PubMed Central 

Google Scholar 
Liao, Y., Smyth, G. K. & Shi, W. FeatureCounts: an efficient general purpose program for assigning sequence reads to genomic features. Bioinformatics 30, 923–930 (2014).CAS 
PubMed 

Google Scholar 
Love, M. I., Huber, W. & Anders, S. Moderated estimation of fold change and dispersion for RNA-seq data with DESeq2. Genome Biol. 15, 550–550 (2014).PubMed 
PubMed Central 

Google Scholar 
Law, C. W., Chen, Y., Shi, W. & Smyth, G. K. Voom: Precision weights unlock linear model analysis tools for RNA-seq read counts. Genome Biol. 15, R29–R29 (2014).PubMed 
PubMed Central 

Google Scholar  More

Reich, P. B. et al. The evolution of plant functional variation: traits, spectra, and strategies. Int. J. Plant Sci. 164, S143–S164 (2003).
Google Scholar 
Cornelissen, J. H. C. et al. A handbook of protocols for standardised and easy measurement of plant functional traits worldwide. Aust. J. Bot. 51, 335–380 (2003).
Google Scholar 
Liu, X. J. & Ma, K. P. Plant functional traits concepts, applications and future directions. Sci. Sin. Vitae 45, 325–339 (2015).
Google Scholar 
Diaz, S., Cabido, M. & Casanoves, F. Plant functional traits and environmental filters at a regional scale. J. Veg. Sci. 9, 113–122 (1998).
Google Scholar 
Kraft, N. J. B., Godoy, O. & Levine, J. M. Plant functional traits and the multidimensional nature of species coexistence. Proc. Natl Acad. Sci. USA 112, 797–802 (2015).ADS 
CAS 
PubMed 
PubMed Central 

Google Scholar 
Barton, K. E. Tougher and thornier: general patterns in the induction of physical defence traits. Func. Ecol. 30, 181–187 (2016).
Google Scholar 
Adler, P. B., Fajardo, A., Kleinhesselink, A. R. & Kraft, N. J. B. Trait-based tests of coexistence mechanisms. Ecol. Lett. 16, 1294–1306 (2013).PubMed 

Google Scholar 
Wright, S. J. et al. Functional traits and the growth–mortality trade-off in tropical trees. Ecology 91, 3664–3674 (2010).PubMed 

Google Scholar 
Wright, I. J. et al. The worldwide leaf economics spectrum. Nature 428, 821–827 (2004).ADS 
CAS 
PubMed 

Google Scholar 
Ruiz-Jaen, M. C. & Potvin, C. Can we predict carbon stocks in tropical ecosystems from tree diversity? Comparing species and functional diversity in a plantation and a natural forest. New Phytol. 189, 978–987 (2011).PubMed 

Google Scholar 
Grubb, P. J. A positive distrust in simplicity-lessons from plant defences and from competition among plants and among animals. J. Ecol. 80, 585–610 (1992).
Google Scholar 
Hanley, M. E., Lamont, B. B., Fairbanks, M. M. & Rafferty, C. M. Plant structural traits and their role in anti-herbivore defence. Perspect. Plant Ecol. 8, 157–178 (2007).
Google Scholar 
Burns, K. C. Spinescence in the New Zealand flora: parallels with Australia. N. Z. J. Bot. 54, 273–289 (2016).
Google Scholar 
Goheen, J. R., Young, T. P., Keesing, F. & Palmer, T. M. Consequences of herbivory by native ungulates for the reproduction of a savanna tree. J. Ecol. 95, 129–138 (2007).
Google Scholar 
Goldel, B., Kissling, W. D. & Svenning, J.-C. Geographical variation and environmental correlates of functional trait distributions in palms (Arecaceae) across the New World. Bot. J. Linn. Soc. 179, 602–617 (2015).
Google Scholar 
Alves-Silva, E. & Del-Claro, K. Herbivory causes increases in leaf spinescence and fluctuating asymmetry as a mechanism of delayed induced resistance in a tropical savanna tree. Plant Ecol. Evol. 149, 73–80 (2016).
Google Scholar 
Cooper, S. M. & Ginnett, T. F. Spines protect plants against browsing by small climbing mammals. Oecologia 113, 219–221 (1998).ADS 
PubMed 

Google Scholar 
Charles-Dominique, T. et al. Spiny plants, mammal browsers, and the origin of African savannas. Proc. Natl Acad. Sci. USA 113, E5572–E5579 (2016).CAS 
PubMed 
PubMed Central 

Google Scholar 
Ratnam, J., Tomlinson, K. W., Rasquinha, D. N. & Sankaran, M. Savannahs of Asia: antiquity, biogeography, and an uncertain future. Philos. Trans. R. Soc. B. 371, 20150305 (2016).
Google Scholar 
Scholes, R. & Archer, S. Tree-grass interactions in savannas. Annu. Rev. Ecol. Syst. 28, 517–544 (1997).
Google Scholar 
Cerling, T. E. Development of grasslands and savannas in East Africa during the Neogene. Palaeogeogr. Palaeoclimatol. Palaeoecol. 97, 241–247 (1992).
Google Scholar 
Brown, R. W. Additions to the flora of the Green River formation. U. S. Geol. Surv. Prof. Paper, U. S. Gov. Print. Off. 154-J, 279–292 (1929).Manchester, S. Oligocene fossil plants of the John Day Formation, Oregon. Or. Geol. 49, 115d–127d (1987).
Google Scholar 
Meyer, H. W. & Manchester, S. R. Oligocene Bridge Creek flora of the John Day Formation, Oregon (Univ. California Press, 1997).Lancucka-Srodoniowa, M. Tortonian flora from the “Gdow Bay” in the south of Poland. Acta Palaeobot. 7, 1–134 (1966).
Google Scholar 
Yuan, J. et al. Rapid drift of the Tethyan Himalaya terrane before two-stage India-Asia collision. Natl Sci. Rev. 8, nwaa173 (2021).PubMed 

Google Scholar 
Spicer, R. A. et al. Why the ‘Uplift of the Tibetan Plateau’is a myth. Natl Sci. Rev. 8, nwaa091 (2021).PubMed 

Google Scholar 
Spicer, R. A. Tibet, the Himalaya, Asian monsoons and biodiversity–In what ways are they related? Plant Divers. 39, 233–244 (2017).PubMed 
PubMed Central 

Google Scholar 
DeCelles, P. G., Kapp, P., Gehrels, G. E. & Ding, L. Paleocene-Eocene foreland basin evolution in the Himalaya of southern Tibet and Nepal: implications for the age of initial India-Asia collision. Tectonics 33, 824–849 (2014).ADS 

Google Scholar 
Royden, L. H., Burchfiel, B. C. & van der Hilst, R. D. The geological evolution of the Tibetan Plateau. Science 321, 1054–1058 (2008).ADS 
CAS 
PubMed 

Google Scholar 
Deng, T., Wu, F. X., Zhou, Z. K. & Su, T. Tibetan Plateau: an evolutionary junction for the history of modern biodiversity. Sci. China Earth Sci. 63, 172–187 (2020).ADS 

Google Scholar 
Favre, A. et al. The role of the uplift of the Qinghai‐Tibetan Plateau for the evolution of Tibetan biotas. Biol. Rev. 90, 236–253 (2015).PubMed 

Google Scholar 
Su, T. et al. A Middle Eocene lowland humid subtropical “Shangri-La” ecosystem in central Tibet. Proc. Natl Acad. Sci. USA 117, 32989–32995 (2020).ADS 
CAS 
PubMed 
PubMed Central 

Google Scholar 
Scientific Expedition Team to the Qinghai-Xizang Plateau. Vegetation of Xizang (Tibet) (Sci. Press, 1988).Liu. X. H. Paleoelevation History and Evolution of the Cenozoic Lunpola basin, Central Tibet. Doctoral thesis (Institute of Tibetan Plateau Research, Chinese Academy of Sciences, 2018).Xiong, Z. Y. et al. The rise and demise of the Paleogene Central Tibetan Valley. Sci. Adv. 8, eabj0944 (2022).CAS 
PubMed 
PubMed Central 

Google Scholar 
Reichgelt, T., West, C. K. & Greenwood, D. R. The relation between global palm distribution and climate. Sci. Rep. 8, 4721 (2018).ADS 
PubMed 
PubMed Central 

Google Scholar 
Farnsworth, A. et al. Paleoclimate model-derived thermal lapse rates: towards increasing precision in paleoaltimetry studies. Earth Planet. Sci. Lett. 564, 116903 (2021).CAS 

Google Scholar 
Spicer, R. A. et al. Why do foliar physiognomic climate estimates sometimes differ from those observed? Insights from taphonomic information loss and a CLAMP case study from the Ganges Delta. Palaeogeogr. Palaeoclimatol. Palaeoecol. 302, 381–395 (2011).
Google Scholar 
Walter, H. Vegetation of the Earth and Ecological Systems of the Geo-biosphere (Springer Berlin Heidelb., 1973).Burley, J. Encyclopedia of Forest Sciences (Acad. Press, 2004).Deng, T. et al. A mammalian fossil from the Dingqing Formation in the Lunpola Basin, northern Tibet, and its relevance to age and paleo-altimetry. Sci. Bull. 57, 261–269 (2012).CAS 

Google Scholar 
Ma, P. F. et al. Late Oligocene-early Miocene evolution of the Lunpola Basin, central Tibetan Plateau, evidences from successive lacustrine records. Gondwana Res. 48, 224–236 (2017).ADS 

Google Scholar 
Hempson, G. P., Archibald, S. & Bond, W. J. A continent-wide assessment of the form and intensity of large mammal herbivory in Africa. Science 350, 1056–1061 (2015).ADS 
CAS 
PubMed 

Google Scholar 
Spicer, R. A. The formation and interpretation of plant fossil assemblages. Adv. Bot. Res. 16, 95–191 (1989).
Google Scholar 
Gibson, D. J. Grasses and Grassland Ecology (Oxford Univ. Press, 2009).Eltringham, S. K. The Hippos: Natural History and Conservation (Princeton Univ. Press, 1999).Jiang, H. et al. Oligocene Koelreuteria (Sapindaceae) from the Lunpola Basin in central Tibet and its implication for early diversification of the genus. J. Asian Earth Sci. 175, 99–108 (2019).ADS 

Google Scholar 
Liu, J. et al. Biotic interchange through lowlands of Tibetan Plateau suture zones during Paleogene. Palaeogeogr. Palaeoclimatol. Palaeoecol. 524, 33–40 (2019).
Google Scholar 
Jia, L. B. et al. First fossil record of Cedrelospermum (Ulmaceae) from the Qinghai-Tibetan Plateau: implications for morphological evolution and biogeography. J. Syst. Evol. 57, 94–104 (2019).
Google Scholar 
Su, T. et al. No high Tibetan Plateau until the Neogene. Sci. Adv. 5, eaav2189 (2019).ADS 
CAS 
PubMed 
PubMed Central 

Google Scholar 
Zhang, Y. L., Li, B. Y. & Zheng, D. A discussion on the boundary and area of the Tibetan Plateau in China. Geol. Res. 21, 1–8 (2002).
Google Scholar 
Yao, T. D. et al. From Tibetan Plateau to Third Pole and Pan-Third Pole. Bull. Chin. Acad. Sci. 32, 924–931 (2017).
Google Scholar 
Spicer, R. A., Farnsworth, A. & Su, T. Cenozoic topography, monsoons and biodiversity conservation within the Tibetan Region: an evolving story. Plant Divers. 42, 229–254 (2020).PubMed 
PubMed Central 

Google Scholar 
Liu, X. H., Xu, Q. & Ding, L. Differential surface uplift: Cenozoic paleoelevation history of the Tibetan Plateau. Sci. China Earth Sci. 59, 2105–2120 (2016).ADS 
CAS 

Google Scholar 
Ding, L., Li, Z. Y. & Song, P. P. Core fragments of Tibetan Plateau from Gondwanaland united in Northern Hemisphere. Bull. Chin. Acad. Sci. 32, 945–950 (2017).
Google Scholar 
Deng, T. & Ding, L. Paleoaltimetry reconstructions of the Tibetan Plateau: progress and contradictions. Natl Sci. Rev. 2, 417–437 (2015).CAS 

Google Scholar 
Li, S. F. et al. Orographic evolution of northern Tibet shaped vegetation and plant diversity in eastern Asia. Sci. Adv. 7, eabc7741 (2021).ADS 
PubMed 
PubMed Central 

Google Scholar 
Ding, L. et al. The Andean-type Gangdese Mountains: Paleoelevation record from the Paleocene–Eocene Linzhou Basin. Earth Planet. Sci. Lett. 392, 250–264 (2014).ADS 
CAS 

Google Scholar 
Deng, T. et al. Review: implications of vertebrate fossils for paleo-elevations of the Tibetan Plateau. Glob. Planet. Change 174, 58–69 (2019).ADS 

Google Scholar 
Westerhold, T. et al. An astronomically dated record of Earth’s climate and its predictability over the last 66 million years. Science 369, 1383–1387 (2020).ADS 
CAS 
PubMed 

Google Scholar 
Nemani, R. R. et al. Climate-driven increases in global terrestrial net primary production from 1982 to 1999. Science 300, 1560–1563 (2003).ADS 
CAS 
PubMed 

Google Scholar 
Hopkins, W. G. Introduction to Plant Physiology (John Wiley & Sons, 1999).Sun, J. M., Liu, W. G., Liu, Z. H. & Fu, B. H. Effects of the uplift of the Tibetan Plateau and retreat of Neotethys ocean on the stepwise aridification of mid-latitude Asian interior. Bull. Chin. Acad. Sci. 32, 951–958 (2017).
Google Scholar 
Zong, G. F. Cenezoic Mammals and Environment of Hengduan Mountains Region (China Ocean Press, 1996).Deng, T. et al. An Oligocene giant rhino provides insights into Paraceratherium evolution. Commun. Biol. 4, 639 (2021).PubMed 
PubMed Central 

Google Scholar 
Young, T. P., Stanton, M. L. & Christian, C. E. Effects of natural and simulated herbivory on spine lengths of Acacia drepanolobium in Kenya. Oikos 101, 171–179 (2003).
Google Scholar 
Karban, R. & Myers, J. H. Induced plant responses to herbivory. Annu. Rev. Ecol. Syst. 20, 331–348 (1989).
Google Scholar 
Huntly, N. Herbivores and the dynamics of communities and ecosystems. Annu. Rev. Ecol. Syst. 22, 477–503 (1991).
Google Scholar 
Asner, G. P. et al. Large-scale impacts of herbivores on the structural diversity of African savannas. Proc. Natl Acad. Sci. USA 106, 4947–4952 (2009).ADS 
CAS 
PubMed 
PubMed Central 

Google Scholar 
Sankaran, M., Augustine, D. J. & Ratnam, J. Native ungulates of diverse body sizes collectively regulate long‐term woody plant demography and structure of a semi‐arid savanna. J. Ecol. 101, 1389–1399 (2013).
Google Scholar 
Staver, A. C. & Bond, W. J. Is there a ‘browse trap’? Dynamics of herbivore impacts on trees and grasses in an African savanna. J. Ecol. 102, 595–602 (2014).
Google Scholar 
Bakker, E. S. et al. Combining paleo-data and modern exclosure experiments to assess the impact of megafauna extinctions on woody vegetation. Proc. Natl Acad. Sci. USA 113, 847–855 (2016).ADS 
CAS 
PubMed 

Google Scholar 
Spicer, R. A. et al. The topographic evolution of the Tibetan Region as revealed by palaeontology. Palaeobio. Palaeoenv. 101, 213–243 (2021).
Google Scholar 
Rowley, D. B. & Currie, B. S. Palaeo-altimetry of the late Eocene to Miocene Lunpola basin, central Tibet. Nature 439, 677–681 (2006).ADS 
CAS 
PubMed 

Google Scholar 
Sun, J. M. et al. Palynological evidence for the latest Oligocene-early Miocene paleoelevation estimate in the Lunpola Basin, central Tibet. Palaeogeogr. Palaeoclimatol. Palaeoecol. 399, 21–30 (2014).
Google Scholar 
DeCelles, P. G., Kapp, P., Ding, L. & Gehrels, G. E. Late Cretaceous to middle Tertiary basin evolution in the central Tibetan Plateau: Changing environments in response to tectonic partitioning, aridification, and regional elevation gain. Geol. Soc. Am. Bull. 119, 654–680 (2007).ADS 

Google Scholar 
Tang, H. et al. Extinct genus Lagokarpos reveals a biogeographic connection between Tibet and other regions in the Northern Hemisphere during the Paleogene. J. Syst. Evol. 57, 670–677 (2019).
Google Scholar 
Wang, T. X. et al. Fossil fruits of Illigera (Hernandiaceae) from the Eocene of central Tibetan Plateau. J. Syst. Evol. 59, 1276–1286 (2021).
Google Scholar 
Del Rio, C. et al. Asclepiadospermum gen. nov., the earliest fossil record of Asclepiadoideae (Apocynaceae) from the early Eocene of central Qinghai-Tibetan Plateau, and its biogeographic implications. Am. J. Bot. 107, 126–138 (2020).PubMed 

Google Scholar 
Xu, Z. Y. The Tertiary and its petroleum potential in the Lunpola Basin, Tibet. Oil Gas. Geol. 1, 153–158 (1980).
Google Scholar 
Zhang, K. X. et al. Paleogene-Neogene stratigraphic realm and sedimentary sequence of the Qinghai-Tibet Plateau and their response to uplift of the plateau. Sci. China Earth Sci. 53, 1271–1294 (2010).ADS 

Google Scholar 
Wu, Y. F. & Chen, Y. Y. Fossil cyprinid fishes from the late Tertiary of north Xizang, China. Vertebrata Palasiat. 18, 15–20 (1980).
Google Scholar 
Wu, F. X., Miao, D. S., Chang, M. M., Shi, G. L. & Wang, N. Fossil climbing perch and associated plant megafossils indicate a warm and wet central Tibet during the late Oligocene. Sci. Rep. 7, 878 (2017).ADS 
PubMed 
PubMed Central 

Google Scholar 
Cai, C. Y., Huang, D. Y., Wu, F. X., Zhao, M. & Wang, N. Tertiary water striders (Hemiptera, Gerromorpha, Gerridae) from the central Tibetan Plateau and their palaeobiogeographic implications. J. Asian Earth Sci. 175, 121–127 (2017).ADS 

Google Scholar 
Low, S. L. et al. Oligocene Limnobiophyllum (Araceae) from the central Tibetan Plateau and its evolutionary and palaeoenvironmental implications. J. Syst. Palaeontol. 18, 415–431 (2020).
Google Scholar 
Bell, A. D. & Bryan, A. Plant Form: An Illustrated Guide to Flowering Plant Morphology (Timber Press, 2008).Zanne, A. E. et al. Three keys to the radiation of angiosperms into freezing environments. Nature 506, 89–92 (2014).ADS 
CAS 
PubMed 

Google Scholar 
Revell, L. J. phytools: an R package for phylogenetic comparative biology (and other things). Methods Ecol. Evol. 3, 217–223 (2012).
Google Scholar 
Paradis, E. & Schliep, K. ape 5.0: an environment for modern phylogenetics and evolutionary analyses in R. Bioinformatics. 35, 526–528 (2019).CAS 
PubMed 

Google Scholar 
Harmon, L. J., Weir, J. T., Brock, C. D., Glor, R. E. & Challenger, W. GEIGER: investigating evolutionary radiations. Bioinformatics. 24, 129–131 (2008).CAS 
PubMed 

Google Scholar 
Maddison, W. P. Confounding asymmetries in evolutionary diversification and character change. Evolution 60, 1743–1746 (2006).PubMed 

Google Scholar 
Forest, C. E., Molnar, P. & Emanuel, K. A. Palaeoaltimetry from energy conservation principles. Nature 374, 347–350 (1995).ADS 
CAS 

Google Scholar 
Valdes, P. J. et al. The BRIDGE HadCM3 family of climate models: HadCM3@ Bristol v1.0. Geosci. Model Dev. 10, 3715–3743 (2017).ADS 
CAS 

Google Scholar 
Gough, D. O. Solar interior structure and luminosity variations. Sol. Phys. 74, 21–34 (1981).ADS 
CAS 

Google Scholar 
Foster, G. L., Royer, D. L. & Lunt, D. J. Future climate forcing potentially without precedent in the last 420 million years. Nat. Commun. 8, 14845 (2017).ADS 
CAS 
PubMed 
PubMed Central 

Google Scholar 
Cox, P. M. Description of the “TRIFFID” Dynamic Global Vegetation Model. 1–16 (Met Office Hadley Centre, 2001).Cox, P., Huntingford, C. & Harding, R. A canopy conductance and photosynthesis model for use in a GCM land surface scheme. J. Hydrol. 212, 79–94 (1998).ADS 

Google Scholar 
McInerney, F. A., Strömberg, C. A. E. & White, J. W. C. The Neogene transition from C3 to C4 grasslands in North America stable carbon isotope ratios of fossil phytoliths. Paleobiology 37, 23–49 (2011).
Google Scholar 
Lu, H. Y. et al. Phytoliths as quantitative indicators for the reconstruction of past environmental conditions in China II: palaeoenvironmental reconstruction in the Loess Plateau. Quat. Sci. Rev. 25, 945–959 (2006).ADS 

Google Scholar  More

The copy numbers for 16S and ITS1 rRNA, and the sequencing depth for all samples are presented in Supplementary File 3 (qPCR data, Sequencing Rarefaction Curves). An average of 14,265.25 reads per housefly sample for the V4 16SrRNA and 16,149.4 reads per housefly sample for the ITS1 were retained after quality filtering. After quality filtering of the egg-laying substrate samples, an average of 10,371.75 reads were retained per sample for the V4 16SrRNA, and an average of 25,479.75 reads were retained per sample for the ITS1 region. The extracted DNA from newly emerged adult houseflies of the Spanish laboratory strain (12 samples in total, newly emerged adults, three replicates from four generations, strain SP100) returned a low copy number for the fungal ITS1 (qPCR data, Supplementary File 3) and a low number of acquired sequencing reads; they were therefore omitted from any further analysis of the fungal microbiota. In addition, the mitochondrial COI phylogeny showed that the Dutch wild-caught strain and the Dutch laboratory strain, which were sampled from the same locality at different times, are in close proximity and form a separate clade from the Spanish lab strain phylotypes (Supplementary File 2).The housefly microbiota alpha-diversity is determined by sampling environmentAbsolute richness (number of ASVs), Shannon index, and Phylogenetic diversity for all housefly strains and developmental stages are shown in Fig. 1. The highest bacterial alpha diversity was observed for the wild-caught housefly population GK0. Strain was an important factor for separating Shannon biodiversity levels both for newly emerged (F = 4.37, P  More