Philippa Kaur

Ant experimentsAnt coloniesSeven colonies of the garden ant L. niger collected from the Soka University and a nearby park were used in this study (Extended Data Table S1). They were placed in plastic cases (35 × 25 × 6 cm). Water was provided ad libitum. They were fed a sucrose solution, and were starved for 2–5 days before the start of the experiment. The colonies were queen-less colonies with 200–700 workers. Aqueous sucrose solution was used as a food resource (bait) in the experiments. The laboratory room where the experiments were performed and the ant colonies were kept was maintained at a temperature of 25–27 °C and a humidity of 60–70%. Artificial lights were also installed in this room.ApparatusWe used an apparatus, called “the main apparatus,” with two paths from the nest to the feeding site (length: 30 cm, width: 2 cm, height: 12 cm for the outward path and 15 cm for the return path) (Fig. 1). This apparatus could separate the outward path (bridge) from the inward path (bridge). Here, the outward path refers to that taken by ants from the nest to the feeding site, whereas the inward path refers to the path taken by ants from the feeding site to the nest.Figure 1The main apparatus used in the three experiments (the main experiment and the comparison experiments 1 and 2). Nests are connected to the experimental apparatus by a slope. In the main experiment, on the outward path, there is ant traffic from the nest to the feeding site on a pheromone trail, and on the inward path, there is ant traffic from the feeding site to the nest on a pheromone trail. In the comparison experiment 1, only a pheromone trail is present on both the outward and inward paths. In the comparison experiment 2, no pheromone trail or ant traffic is present on both the outward and inward paths.Full size imageTwo important features of this apparatus were as follows: firstly, it allowed ants to only enter the outward path from the nest. A rat-guard structure at the end of the inward path prevented the ants on the outward path from entering the inward path (Extended Data Fig. S1A). Secondly, we installed a vertical structure at the end of the outward path (height: 4 cm). After climbing the vertical structure, ants were not allowed to return to the outward path (Extended Data Fig. S1B). Moreover, we installed partitions on the feeding site, which also prevented ants from returning to the outward path after reaching the feeding site (Extended Data Fig. S1C). After entering the feeding site, ants had to pass through a narrow gap (width: 0.5 cm) created by the partition. No visual cues were offered as the apparatus was surrounded on all four sides by plastic walls.In this experiment, we made another apparatus for a single ant (target ant), which would be joining the ant trail on the main bridges (Extended Data Fig. S2). This apparatus, called “the confluence device,” was a detachable device that could be connected at right angles to the outward and inward bridges of the main apparatus. To connect this device to the outward bridge, we made the confluence path of this device under the inward bridge of the main apparatus, since the outward bridge was lower than the inward bridge. Thus, we made a slope on the outward confluence path connected to the outward bridge of the main apparatus. Further, because placing the ants directly on the sidewalk sometimes caused them to fall off the sidewalk owing to panic, we constructed a free space and a wall (height: 5 cm) in the middle of the confluence device on which the ants were placed calmly. Owing to this modification, we could let each target ant calm down and then access the main bridge whenever they wanted to. The apparatus used in this experiment was made of white plastic plates.Pheromone trail with ant trafficThis main experiment was limited to once a day for each colony. A sucrose solution was dripped into the feeding site. Target ants, which were walking on a plastic case as foragers, had been moved from their nests to another case immediately before a trail of (nontarget) ants was formed. Thus, dozens of ants were moved in advance to the case to be used as target ants. Subsequently, a trail of (nontarget) ants was formed from the nest to the main apparatus. Considering that it took some time for the ants that had finished foraging and returned to the nest to recruit their mates, the ants were left for approximately 40 min to an hour until a permanent ant trail was formed. It was difficult to form an ant trail immediately after the start of the experiment since no foraging pheromones could be produced in the first foraging trip on the outward path and since experienced foraging ants may make foraging pheromones on the outward path2,21,22. The target ants were allowed to enter bridges of the main apparatus after the establishment of a permanent ant trail. At that time, trails of individual target ants were started. Target ants were allowed to join at right angles to the path on the apparatus, one by one from the confluence device. Individual target ants were allowed to enter the main apparatus at four different points: (1) Left-Left (LL), located at the left side of the center of the outward path. The outward path was on the left side, whereas the inward path was on the right side for the experimenter when seen from the nest. (2) Left–Right (LR), located at the right side of the center of the outward path. (3) and (4) Right-Left (RL) and Right-Right (RR), located at the left and right sides of the center of the inward path, respectively (Fig. 2). We had set these four points to check if target ants tended to turn their body to a certain direction when entering the main bridges, regardless of the movement direction of the other ants. A video camera (Panasonic, AVCHD 30fps) was used to record the migration of ants to the feeding site or nest. Videos were taken from above, and target ants were used only once.Figure 2Four joining points (LL, LR, RL, and RR) and the confluence device (joining device). The confluence (joining) device was connected at right angles to the center of the outward and inward bridges of the main apparatus. Here, the LR version is shown as an example.Full size imageThe goal lines were set at 15 cm from the center of the main paths. We checked the side (nest side or feeding site side) from which a target ant passed the goal line.Pheromone trail with no ant trafficThis comparison experiment 1 was limited to once a day for each colony. Dozens of ants were moved in advance to another case to be used as target ants in a similar manner to the main experiment. The (nontarget) ants were left for about 40 min to an hour until a permanent ant trail was formed. Subsequently, we removed all the ants from the device. Then, target ants were allowed to enter on the side path one by one. In this case, we left the bait in place to control this experiment under the same condition as the main experiment. As the pheromone trail was created on the outward path as well as on the inward path, it was the only decision-making factor for the ants to join at the main path (outward/inward paths). We checked the side (nest side or feeding site side) from which a target ant passed the goal line in a similar manner to the main experiment.No pheromone trail or ant trafficThis comparison experiment 2 was limited to once a day for each colony. Dozens of ants were moved in advance to another case to be used as target ants in a similar manner to the main experiment. This experiment was conducted to investigate ant behavior under the following two conditions: (1) no ant trails and (2) no pheromones trails. The bait was in place in the same manner. We checked the side (nest side or feeding site side) from which a target ant passed the goal line in a similar manner to the main experiment. After each trial (the target ant passed the goal line), we wiped the apparatus with ethanol solution before the next target ant was allowed to enter the main paths.AnalysisThe goal lines were set at 15 cm from the center of the main paths. We checked which goal side the target ants reached the goal line on each trial. A reverse run referred to the goal to the nest on the outward path and the goal to the feeding site on the inward path. A normal run referred to the goal to the feeding site on the outward path and the goal to the nest on the inward path.In some cases of the main experiment, foraging (nontarget) ants that could not reach the feeding site on their outward path or could not return to the nest on their inward path would be against the ant flows. On the outward path, we considered that the ants conducted a “reverse flow” if the position of their heads was on the nest side compared with the position of their stomach. If not, we defined that the ants conducted a “normal flow” (Extended Data Fig. S3). On the inward path, we defined that the ants conducted a “reverse flow” if the position of their head was on the feeding site side compared with the position of their stomach. If not, we defined that the ants conducted a “normal flow” (Extended Data Fig. S3). We focused on the target ants that came in contact with ants with normal flow. Therefore, if an ant with reverse flow was located within 10 cm of the target ant, that trial was excluded from the analysis.Furthermore, we also evaluated if target ants coming in contact with foraging (nontarget) ants immediately after entering the trail would affect the goal choice. Therefore, we conducted an analysis focusing on the contact using the data from the main experiment. We examined whether or not the target ant made contact with other foraging ants until it passed a point 2 cm from the center of the path. As already mentioned, if the target ant came in contact with another ant moving against the normal flow of the ant trail, this contact was excluded from the counts. Moreover, we also excluded cases in which the body of target ants was on a point 2 cm from the center of the path by visual evaluation. Thus, we examined the goal choice of target ants by focusing on whether or not they came in contact with other ants immediately after joining the main bridges.We also conducted a preliminary experiment using a single path apparatus to investigate bi-directional trail behaviour. Please see the Extended Data File S1.Model descriptionThe models were coded using the C programming language. The model description follows the Overview, Design concepts, and Details protocol23,24.PurposeThe purpose of the model was to examine the mechanistic understanding of our findings. We adopted an action of target agents obtained from our ant experiments and compared it with another action of target agents on a trail that was contrary to the fact. To be more precise, target agents were allowed to obey an alignment rule in which they tended to move in the same direction with other agents. We named the former model as the reverse-rule model and the latter model as the alignment-rule model. By doing so, we could find the significance of our findings from ant experiments.Entities, state variables, and scalesWe developed two different models (reverse-rule model and alignment-rule model) that included two types of entities: agents and cells. The agent has the state variable Navigational state, which has two values: Navigational state = {wandering, foraging}. The cell has the state variable Pheromone; this value represents the amount of pheromones in each cell. We used a 2D lattice field and set a straight bridge with 61 cells × 5 cell sizes. We also set goal lines at x-coordinate =  − 30 and 30. If the agents reached coordinates satisfying their x-coordinate =  − 30 or 30, they were removed from the system. If the agents reached y-axis boundaries, their movement direction was restricted. Each trial continued until the target agent reached one of the two goal lines. However, trials were forcibly finished if the target agent never reached any goal line by t = 500-time steps. In total, we conducted 1000 trials.Process overview and schedulingAt the beginning of each trial, an artificial target ant (Navigational state = wandering) was introduced at the center of an artificial simulation field. Foraging agents (Navigational state = foraging) were randomly distributed on the simulation field in advance.Agents on the simulation field selected one direction from two directions (+ x and − x) on each time step and updated their positions. Briefly, an agent at coordinate (x, y) selected one direction from two directions (+ x and − x) and updated its position with one of the three coordinates—(x − 1, y), (x − 1, y + 1), or (x − 1, y − 1)—if it selected the − x direction, or—(x + 1, y), (x + 1, y + 1), or (x + 1, y − 1)—if it selected the + x direction by scanning pheromones on these three coordinates. For example, if an agent at coordinate (0, 2) decided to move in + x direction at one time, the position of this agent was replaced with one of (1, 3), (1, 2) and (1, 1) from (0, 2) by scanning pheromones on these three coordinates. The target agent selected the − x/ + x direction with equal probability on each time step until it met the foragers. In contrast, foraging agents tended to decide to move in the − x direction on each time step with a high probability and therefore they tended to select the − x direction for position updating. Foraging agents deposited pheromones before leaving the current cell (see submodel entitled “Position updating” and submodel entitled “Pheromone updating”). In contrast, the target agents did not deposit pheromones.Using above submodels, artificial ants sometimes met other agents. If the target agent (Navigational state = wandering) met the foragers (Navigational state = foraging), the target agent tended to select one direction from two directions (+ x and − x) on each time step thereafter with a high probability, which was dependent on which direction the met foragers came from. More strictly, in the reverse-rule model, the target agent tended to move in an opposite direction from the foragers if it met the foragers coming from the opposite direction. On the contrary, the target agent in the alignment-rule model tended to move in the same direction with foragers if it met the foragers moving in the same direction (see submodel entitled “The interaction between the target agent and foragers”). For example, in the reverse-rule model, if the target agent at coordinate (x, y), whose previous coordinate was (x − 1, y), met the forager coming from the opposite direction, whose previous coordinate was (x + 1, y), the target agent decided to move in + x direction on each time step thereafter with a high probability until similar events occurred.Design conceptThe mean goal time was the emergent property of the model. Sensing was important as the agents scanned the pheromone concentrations. Stochasticity was used to determine in which direction the agent moved and to select one cell using the pheromone concentrations.InitializationWe set a single agent (target agent) on the coordinate (0, 2) and its Navigational state was set to wandering (Extended Data Fig. S5A). We also set N foraging agents on the bridge whose Navigational state was set to foraging. Therefore, N + 1 agents were on the test field at the beginning of each trial. A target agent was the agent k = 0, whereas foraging agents were agents k = 1, 2, …, N. These foragers were randomly distributed on the bridge. Thus, x(k) (in) {n |− 30 ≤ n ≤ 30, n is an integer} and y(k) (in) {n | 0 ≤ n ≤ 4, n is an integer} for k  > 0.Foraging agents were set to move in the -x direction (Direction(k) for k  > 0 = − x). On the other hand, the target agent randomly chose one direction from two directions at the beginning of each trial (Direction(0) was set to + x or − x with equal probability). Herein, Direction(k) can be − x or + x, which implies bias in the movement direction. The parameter prob(k) indicates the probability of moving in Direction(k). The target agent selected the − x/+ x direction with equal probability on each time step until it met the foragers. Therefore, the parameter prob was set to 0.50 for the target ant (prob(0) = 0.50), whereas prob was set to 0.80 for foraging agents (prob(k) = 0.80 for k  > 0). The amount of pheromones on each cell was set to 1 at the beginning of each trial (pheromone(x, y) = 1) and the pheromone evaporation rate q was set to 0.99.The model descriptions are explained using submodels. A Submodel: the interaction between the target agent and foragers causes differences between two rules (the reverse-rule model and the alignment-rule model).SubmodelsSubmodel: the interaction between the target agent and foragersThe parameters Direction(0) and prob(0) were replaced with new ones whenever the following events occurred.In the reverse-rule model, for any agent k (k  > 0),Herein, (xt(k), yt(k)) indicates the x–y-coordinate for the agent k at time t. Furthermore, (xt(0), yt(0)) = (xt(k), yt(k)) means that the target agent and the agent k occupy the same cell at time t while (xt(0) − xt−1(0)) × (xt(k) − xt−1(k)) =  − 1 indicates that the target agent meets the agent k came from the opposite direction. The target agent replaces Direction(0) with an opposite direction from the forager k (see Extended data Fig. S5B).In the alignment-rule model, for any agent k (k  > 0),(xt(0) − xt−1(0)) × (xt(k) − xt−1(k)) = 1 indicates that the target agent meets the agent k came from the same direction. The target agent replaces Direction(0) with a same direction with the forager k (See Extended Data Fig. S5B).In the reverse-rule model, these events are driven from the experimental observations of real ants. Target ants appear to move against the trail and seem to move straight by contacting those other nestmates that come from the opposite direction. Also, target ants seem to select the reverse goal even if physical contact with ant nestmates does not occur immediately after entering the bridge. So, regarding parameter replacements, we did not consider the position at which the target agent met another agent. Note that foraging agents did not change these parameters until the end of each trial. Further, Direction(0) can be replaced with − x from + x and vice versa whenever the target agent meets foragers that come from the opposite direction.In the alignment-rule model, the target agent tends to move in the same direction with other agents. This is contrary to the experimental observations of real ants.Submodel: position updatingFor all k agents (k = 0–N), the movement direction and position updates are shown as follows (Extended Data Fig. S5C);Here, rnt(k) indicates a random number. Thus, rnt(k) (in) [0.00, 1.00].Prob(0) for the target agent is initially set to 0.50. Therefore, the target agent selects one direction from the two (− x and + x) on each time step randomly before the condition described in submodels—the interaction between the target agent and foragers is satisfied. On the other hand, foraging agents select − x direction with a high probability (= Prob(k)) on each time step. After selecting one direction from two (− x and + x), agents scan three cells in the direction of movement. Using pheromone concentrations on those three cells, they update their positions.If agents reach coordinates satisfying their y-coordinate = 4 or 0, those agents update their position by selecting not three but two coordinates since they are located on the edges of the bridge.Submodel: pheromone updatingForaging agents (k  > 0) deposited pheromones on the current cell when leaving that cell.Then, pheromones are evaporated using the evaporation rate q.For each time iteration, these submodels operated in the following order.STEP 1: The interaction between the target agent and foragers.STEP 2: Position updating.STEP 3: Pheromone updating.AnalysisTo check the accuracy of our model, we counted which goal side the target agent entered the goal line from using the reverse-rule model by setting N = 9. If the target agent passed the goal line at x-coordinate =  − 30 (30), we considered that it reached the normal (reverse) goal. Note that trials in which the target agent never reached any goal lines by t = 500 were excluded from this analysis. Furthermore, to investigate the adaptability of the reverse run mechanism, we examined the time until the target agent reached the goal lines using the reverse-rule model and the alignment-rule model. Herein, we set two different conditions with respect to the number of foraging agents (N = 4 and 9). More

Study subjects and measurement of skin physiological parametersWe analyzed skin microbiome and mycobiome from cheeks and foreheads of healthy younger (19–28 years old, Y-group) and older (60–63 years old, O-group) Korean women who were free from cutaneous disorders (Table 1 and Supplementary Table S1). All 61 subjects had been living in Seoul, Korea, for more than 3 years with normal skin conditions. We preferentially selected those who had sebum secretion greater than 30 arbitrary units and moisture greater than 50 arbitrary units in both groups. Among the measurements of moisture content, pH, sebum content, and transepidermal water loss (TEWL), only sebum and TEWL decreased significantly in the O-group compared to the Y-group in the cheeks (P = 2.25e−06, Wilcoxon rank-sum test; P = 0.019, Welch two-sample t test) and forehead (P = 1.33e−06, Wilcoxon rank-sum test; P = 0.003, Welch two-sample t test). Whereas no significant differences were found in the average values for moisture (cheeks: Y-group, 59.9; O-group, 56.6; forehead: Y-group, 61.1; O-group, 58.7) and pH (cheeks: Y-group, 6.0; O-group, 5.8; forehead: Y-group, 6.0; O-group, 5.6) between the two age groups.Table 1 Characteristics of subjects for aged related skin microbiome and mycobiome study.Full size tableComparisons in cheek and forehead microbiome and mycobiome between the two age groupsWe analyzed bacterial communities from 27 Y-group samples (cheeks, n = 13; forehead, n = 14) and 24 O-group samples (cheeks, n = 12; forehead, n = 12) and fungal communities from 28 Y-group samples (cheeks, n = 15; forehead, n = 13) and 32 O-group samples (cheeks, n = 16; forehead, n = 16), except for samples that were eliminated from the Illumina Mi-Seq sequencing due to low sequence reads (bacteria,  3. 0) (Fig. 4). Pathways belonging to the metabolism category were dominant in each age group. In the cheek of the Y-group, pathways involved in energy metabolism by bacteria, such as glycolysis/gluconeogenesis, citrate cycle, pentose phosphate pathway, fructose and mannose metabolism, galactose metabolism, d-alanine metabolism, and thiamine metabolism, were predominant, whereas in the cheek of the O-group, degradation-related pathways, such as fatty acid degradation, synthesis and degradation of ketone bodies, benzoate degradation, and chloroalkane and chloroalkene degradation, were predominant. In the forehead of the Y-group, glycolysis/gluconeogenesis, pentose phosphate pathway, fructose and mannose metabolism, galactose metabolism, d-glutamine and d-glutamate metabolism, d-alanine metabolism, and thiamine metabolism pathway were significantly more abundant, whereas in the forehead of the O-group, fatty acid degradation, synthesis and degradation of ketone bodies, valine/leucine and isoleucine degradation, and limonene/pinene degradation pathway were significantly more abundant.Figure 4Heat map for significantly different predicted functional pathways on (a) cheeks and (b) foreheads of Korean women by age based on LEfSe analysis (LDA score  > 3.0).Full size imageThe metabolism pathway for biotin, a water-soluble vitamin that is effective for skin health and essential for keratin production15, was more prevalent in the cheek and forehead of the Y-group. Interestingly, the metabolism pathway for lipoic acid, which is known to possess beneficial effects against skin aging and is used widely in cosmetic and dermatological products16,17, was significantly higher in the foreheads of the Y-group. We tracked the specific ASVs possessing these pathways, in both biotin metabolism and lipoic acid metabolism, Cutibacterium sp. (ASV2136 and ASV2130) and Staphylococcus sp. (ASV3008) were predicted to have the top three relative abundances in KOs. The relative abundances in biotin metabolism and lipoic acid metabolism of Cutibacterium sp. (ASV2136) were 24.9% and 26.1%, respectively. The relative abundances for each pathway for Staphylococcus sp. (ASV3008) were 10.2% and 18.7%, and for Cutibacterium sp. (ASV2130), they were 9.3% and 10.0%, respectively. We confirmed these two pathways in the genome of skin bacteria, C. acnes (Supplementary Fig. S2). These additional analyses support the reliability of the function in the skin environment of Cutibacterium. Interestingly, from the LEfSe result, Cutibacterium sp. (ASV2136) had a significantly higher abundance in the cheek and forehead microbiome of the Y-group. The pathway of biosynthesis of lipopolysaccharide, also known as bacterial endotoxins, showed higher abundance in the cheek and forehead microbiome of the O-group. The ASVs that contribute to inferring the LPS biosynthesis pathway were identified as Paraburkholderia sp. (ASV5030) and B. vesicularis (ASV4155). Also, pathways related to antibiotic biosynthesis (biosynthesis of vancomycin group antibiotics) and bacterial motility (bacterial chemotaxis and flagellar assembly; both belonging to the cellular processes category) were prominent in the cheek and forehead of the O-group. PICRUSt2 analysis implied that, regardless of skin site differences, the potential functions of the microbial community that compose the skin microbiome were similar according to age.Network analysis on cheek and forehead microbiome and mycobiomeWe performed SParse InversE Covariance estimation for Ecological Association Inference (SPIEC-EASI) analysis to evaluate the overall network of the skin microbes. The results of network density (D) on 81 cheek and 87 forehead ASVs, calculated using the ratio of the number of edges, showed higher network density in the skin microbiome of the Y-group (D = 0.015 and D = 0.001, in cheek and forehead, respectively) than the O-group (D = 0.007 and D = 0.007, respectively) (Fig. 5). To examine network correlation between bacteria and fungi, network density for Bacteria–Fungi (DBF) was calculated by the actual number of edges and a potential number of edges in a correlation ([bacterial nodes × fungal nodes]/2). We confirmed higher network density in the cheek of the Y-group (DBF = 0.008) than the O-group (DBF = 0) and edges of the major bacterial and fungal taxa, such as Staphylococcus sp. (ASV3008)—M. sympodialis (ASV500) and Roseomonas sp. (ASV4088)—M. restricta (ASV482), were observed in the cheek of the Y-group. In the forehead, edges of Methylobacterium sp. (ASV4314)—M. globosa (ASV454), Methylobacterium sp. (ASV4314)—Zygosaccharomyces rouxii (ASV208), and Venionella sp. (ASV3575)—M. sympodialis (ASV500) were observed in the Y-group, and edges of Cutibacterium sp. (ASV2107)—M. globosa (ASV461), Staphylococcus sp. (ASV3024)—M. arunalokei (ASV446), and Methylobacterium (ASV4314)—M. dermatis (ASV448) were observed in the O-group (DBF = 0.004). We found a network between bacteria and fungi with different kingdom levels in the skin microbiome, and especially, we confirmed that different genus or species level microbe was involved in the microbial network according to skin location and Y-, O-group.Figure 5Network analysis of the ASVs on (a) cheeks and (b) forehead of Korean women. Each node represents the ASV and the size of the node is based on relative abundance of each ASV. Color markings indicate the major taxa except for unidentified bacteria or fungi. Shapes represent the level of kingdom, Bacteria (bold) and Fungi (dotted line). The ASVs were selected for bacterial ASVs found in more than half of all samples on the cheeks and forehead, respectively, and for the fungal ASVs with a relative abundance of more than 0.1% in each of the cheeks and forehead samples. The D value is network density calculated using the ratio of the number of edges.Full size image More

Thank you for visiting nature.com. You are using a browser version with limited support for CSS. To obtain
the best experience, we recommend you use a more up to date browser (or turn off compatibility mode in
Internet Explorer). In the meantime, to ensure continued support, we are displaying the site without styles
and JavaScript. More

Epithemia isolation and cultureThe Epithemia cells were isolated from 0.5 L of seawater collected from depths of 25, 75, and 100 m in the North Pacific Subtropical Gyre (22°45′ N, 158°00′ W). Seawater was collected during the near-monthly Hawaii Ocean Time-series (HOT) expeditions to the long-term monitoring site Station ALOHA (water depth ca. 4800 m) in October 2014 (HOT cruise #266) and February–July 2019 (HOT cruises #310–313). Serial dilution (unialgal strains UHM3202, UHM3203, UHM3204) or micropipette isolation of single cells (clonal strains UHM3200, UHM3201, UHM3210, UHM3211) were used to establish the Epithemia cultures, which were grown in a seawater-based, low-nitrogen medium. Filtered (0.2 µm) and autoclaved, undiluted Station ALOHA seawater was amended with 2 μM EDTA, 50 nM ferric ammonium citrate, 7.5 μM phosphoric acid, trace metals (100 nM MnSO4, 10 nM ZnCl2, 10 nM Na2MoO4, 1 nM CoCl2, 1 nM NiCl2, 1 nM Na2SeO3), vitamins (50 μg/L inositol, 10 μg/L calcium pantothenate, 10 μg/L thiamin, 5 μg/L pyridoxine HCl, 5 μg/L nicotinic acid, 0.5 μg/L para-aminobenzoic acid, 0.1 μg/L folic acid, 0.05 μg/L biotin, 0.05 μg/L vitamin B12), and 106 μM Na2SiO3. Although not tested here, simpler formulations of diazotroph media such as PMP40 or RMP41 may also be suitable for growing Epithemia, when made with 100% seawater and adding Na2SiO3. The cultures were subsequently incubated at 24 °C on a 12:12 h light:dark cycle with 50–100 μmol quanta m−2 s−1 using cool white fluorescent bulbs. All E. pelagica and E. catenata symbioses were stable under these medium and incubation conditions. E. pelagica was successfully isolated from at least one of the three depths that were targeted during each sampling occasion.Morphological observationsEpithemia living and fixed cells were imaged by light and epifluorescence microscopy using a Nikon Eclipse 90i microscope at 40×–60× magnification. Diatom cell sizes were determined using >60 live, exponentially growing cells, imaged in either valve view (E. pelagica) or girdle view (E. catenata). Endosymbiont (spheroid body) cell sizes were averaged from DNA-stained cells for E. pelagica UHM3200 (n = 78) and E. catenata UHM3210 (n = 91), imaged by epifluorescence microscopy after preparing samples as follows: Epithemia cells were fixed in 4% glutaraldehyde for 30 min, pelleted at 1000 × g for 1 min, the supernatant was exchanged with 0.5% Triton X-100 (in autoclaved filtered seawater), samples were incubated for 10 min with gentle agitation, cells were then pelleted at 4000 × g for 1 min, supernatant was exchanged with autoclaved filtered seawater and fixed in 4% glutaraldehyde, and samples were stained with 1× final concentration of SYBR Gold nucleic acid stain (Invitrogen, cat. # S11494) for 2 h. For routine observations of endosymbionts (e.g., determining presence/absence and number per host cell), osmotic shock was used to disrupt the cell contents of diatom host cells and improve visualization of the endosymbionts. This was achieved by gently pelleting cells and exchanging the medium with either ultrapure water or 2–3 M NaCl solution, followed by immediate observation. While this is a simple technique for detecting and visualizing endosymbionts (Fig. 1c, f), it does not accurately represent the natural location of endosymbionts within the host cells, as seen when compared to fixed cell preparations for epifluorescence microscopy (Fig. 1n, o). To assess the presence of fluorescent photopigments in endosymbiont cells, live host cells were pelleted at 4000 × g for 5 min and crushed using a microcentrifuge tube pestle (SP Bel-Art, cat. # F19923-0000) to release the endosymbionts. The crushed pellet was resuspended in 75% glycerol containing live Synechococcus WH7803 cells (positive control for fluorescence), and samples were observed by epifluorescence microscopy using filter cubes appropriate for observing phycoerythrin (EX: 551/10, BS: 560, EM: 595/30) and chlorophyll (EX: 480/30, BS: 505, EM: 600LP).The loss of endosymbionts from Epithemia cultures (UHM3200 and UHM3210) was observed after propagating cells for four months in nitrogen-replete medium (K)18, where approximately 5–10% of the culture was transferred to fresh medium about every two weeks. Observations were only made at the end of the four-month period. Endosymbionts were not observed growing freely in these cultures, and the absence of endosymbionts within host cells was confirmed by the failure to observe spheroid bodies by light microscopy after osmotic shock of the diatoms, as well as a failure to amplify the endosymbiont SSU (16S rRNA) and nifH genes from cellular DNA extracts. PCR reactions were performed in parallel with DNA extracts from control cultures (grown in low-nitrogen medium), using the same template DNA amount (10 ng) and PCR conditions (see methods for Marker gene sequencing and phylogenetics).Ultrastructural observations by electron microscopy (EM) were conducted for E. pelagica UHM3200 and E. catenata UHM3210. EM preparations of diatoms typically involve the oxidative removal of organic matter to uncover the fine details of frustule ultrastructure. However, in the case of E. catenata, oxidatively cleaned cells lacked structural integrity, leading to collapsed frustules when dried and viewed by scanning EM (SEM). For this reason, both species were prepared for SEM with and without (Fig. 1a, d) the oxidative removal of organic matter, and cleaned E. catenata frustules were further analyzed by transmission EM (TEM). To remove organic matter, 100 mL of exponentially growing culture was pelleted by centrifugation at 1000 × g for 10 min and resuspended in 30% H2O2. Cells were boiled in H2O2 for 1–2 h, followed by rinsing cells six times in ultrapure water by sequential centrifugation at 1000 × g for 10 min and resuspension of cell pellets. Suspensions of the cleaned cells were dried on aluminum foil and mounted on aluminum stubs with double-sided copper tape. For some E. catenata SEM preparations, the cleaned frustules were dehydrated in an ethanol dilution series and exchanged into hexamethyldisilazane (HMDS) prior to drying on aluminum foil; this was to minimize the collapse of frustules resulting from drying. To prepare cells with organic matter intact, 25 mL of exponentially growing culture was mixed with an equal volume of fixative solution (5% glutaraldehyde, 0.2 M sodium cacodylate pH 7.2, 0.35 M sucrose, 10 mM CaCl2) and incubated overnight at 4 °C. Cells were gently filtered onto a 13 mm diameter 1.2 μm pore size polycarbonate membrane filter (Isopore, Millipore Sigma), washed with 0.1 M sodium cacodylate buffer (pH 7.4, 0.35 M sucrose), fixed with 1% osmium tetroxide in 0.1 M sodium cacodylate (pH 7.4), dehydrated in a graded ethanol series, and critical point dried. Filters were mounted on aluminum stubs with double-sided conductive carbon tape. All SEM stubs were sputter coated with Au/Pd, prior to observing on a Hitachi S-4800 field emission scanning electron microscope at the University of Hawai’i at Mānoa (UHM) Biological Electron Microscope Facility (BEMF). Cleaned E. catenata cells were prepared for TEM by drying a drop of sample on a formvar/carbon-coated grid and observing on a Hitachi HT7700 transmission electron microscope at UHM BEMF.Additional light microscopy of hydrogen-peroxide cleaned frustules was conducted for E. pelagica UHM3201 and E. catenata UHM3210. Samples were mounted in Naphrax (PhycoTech, Inc., cat. # P-Naphrax200) and observed at 100× using an Olympus BX41 Photomicroscope (Olympus America Inc., Center Valley, Pennsylvania) with differential interference contrast optics and an Olympus SC30 Digital Camera at California State University San Marcos.A key to the strains used in each micrograph is provided in Supplementary Table 2.Marker gene sequencing and phylogeneticsFor each Epithemia strain, 25–50 mL of culture was pelleted at 4000 × g for 10 min, and DNA was extracted from the pellet using the ZymoBIOMICS DNA Miniprep Kit (Zymo Research, cat. # D4300). Marker genes were amplified with the Expand High Fidelity PCR System (Roche, cat. # 4743733001), using conditions previously described for genes SSU encoding 18S rRNA (Euk328f/Euk329r)42, LSU encoding 28S rRNA (D1R/D2C)43, rbcL (rbcL66+/dp7−)44,45, psbC (psbC+/psbC−)44, and cob (Cob1f/Cob2r)21. For the endosymbionts, a partial sequence for the SSU (16S rRNA) gene was amplified using a primer set targeting unicellular cyanobacterial diazotrophs, CYA359F/Nitro821R46,47, and the nifH gene was amplified using new primers specific to the nifH of Cyanothece-like organisms, ESB-nifH-F (5′-TACGGAAAAGGCGGTATCGG-3′) and ESB-nifH-R (5′-CACCACCAAGRATACCGAAGTC-3′), with a 55 °C annealing temperature and 75 s extension time. All primers were synthesized by Integrated DNA Technologies (IDT). Amplified products were cloned and transformed into E. coli using the TOPO TA Cloning Kit for Sequencing (Invitrogen, cat. # K457501), and plated colonies were picked and grown in Circlegrow medium (MP Biomedicals, cat. # 113000132). Plasmids were extracted with the Zyppy Plasmid Miniprep kit (Zymo Research, cat. # D4019) and sequenced from the M13 vector primers using Sanger technology at GENEWIZ (South Plainfield, NJ). For the diatom SSU (18S rRNA) gene, sequencing reactions were also performed using the 502f and 1174r primers48.Phylogenetic trees (Fig. 2) were inferred using concatenated alignments for both diatom host genes (SSU encoding 18S rRNA, psbC, rbcL) and endosymbiont genes (SSU encoding 16S rRNA, nifH). For each gene, nucleotide sequences were aligned using MAFFT v7.45349 (L-INS-i method), and sites with gaps or missing data were removed. An appropriate nucleotide substitution model was selected for each gene alignment using jModelTest v2.1.1050. Bayesian majority consensus trees were inferred from the concatenated alignments using MrBayes v3.2.751 with two runs of 4–8 chains, until the average standard deviation of split frequencies dropped below 0.01. Maximum likelihood bootstrap values were generated for the Bayesian tree using RAxML v8.2.1252, implemented with 1000 iterations of rapid bootstrapping. To further analyze the phylogenetic position of the new Epithemia species in the broader context of Surirellales and Rhopalodiales diatoms, individual gene trees (SSU encoding 18S rRNA, LSU, rbcL, psbC, and cob; Supplementary Figs. 13–19) were constructed from sequences aligned using MAFFT (automatic detection method) and trimmed using trimAl v1.253 (gappyout method). rRNA gene phylogenies were also inferred using sequences aligned according to the global SILVA alignment for SSU and LSU genes using SINA54, which were either left untrimmed in the case of the LSU gene or trimmed to remove highly variable positions (SINA’s “012345” positional variability filter) and gappy positions (trimAL v1.2, gappyout method) in the case of the SSU gene. These trimming strategies were selected based on their ability to maximize the monophyly of the previously described Rhopalodiales clade and minimize the separation of known conspecific strains, such as the strains of E. pelagica described here. All gene phylogenies were inferred using the Bayesian methods described above. To investigate the level of support for constrained tree topologies placing E. catenata within or outside of the genus Epithemia and family Rhopalodiaceae, SH55 and AU56 statistical tests were performed in IQ-TREE 257 (implementing ModelFinder58) using all alignments from the individual gene trees (Supplementary Table 3).Given E. catenata’s unusual morphology, test trees were inferred with the inclusion of diatom sequences from orders Bacillariales (Nitzschia, Pseudo-nitzschia), Cymbellales (Didymosphenia), Naviculales (Amphiprora, Navicula, Pinnularia), and Thalassiophysales (Amphora, Halamphora, Thalassiophysa); however, E. catenata was consistently placed within Rhopalodiales, and these trees were not pursued further.An additional nifH phylogeny was constructed using all environmental sequences from NCBI’s non-redundant nucleotide (nt) database >300 bp and sharing >95% nucleotide sequence identity with EpSB and EcSB nifH sequences (Supplementary Fig. 23), including 51 environmental sequences from prior studies investigating marine diazotrophs34,59,60,61,62,63,64,65,66. Environmental nifH sequences were aligned to the previously generated nifH sequence alignment using MAFFT (automatic method detection and addfragments options), and the best-scoring maximum likelihood phylogeny was inferred using RAxML with 1000 iterations of rapid bootstrapping. NCBI accession numbers for all tree sequences are in the Source Data file.Analysis of Epithemia endosymbiont nifH sequences in environmental datasetsNucleotide sequences for EpSB and EcSB nifH were queried against NCBI’s non-redundant nucleotide (nt) database using webBLAST67 (megablast; https://blast.ncbi.nlm.nih.gov/) and SRA databases for nifH amplicon sequencing projects from the marine environment using the SRA Toolkit68 (dc-megablast, with database validation using vdb-validate; https://github.com/ncbi/sra-tools). Database hits with 98–100% nucleotide identity over an alignment of the entire subject sequence (BLAST alignment length = subject sequence length) were identified, and the associated sample’s latitude and longitude coordinates (where available) were mapped. Coordinates were also mapped for metagenome and metatranscriptome samples containing matches to unigene MATOU-v1_93255274 from the Marine Atlas of Tara Oceans Unigenes69, a unigene that shares 100% identity over the entire length of the EpSB UHM3202 nifH sequence and >99.4% identity with all other EpSB nifH sequences.The presence of EpSB and EcSB nifH sequences was examined in metagenomes prepared from sinking particles collected at 4000 m depth at Station ALOHA27,28. The sinking particles were collected during intervals of 12, 10, and 8 days during 2014, 2015, and 2016, respectively, using a McLane sediment trap equipped with a 21-sample bottle carousel. The presence of EpSB and EcSB nifH sequences in the metagenomes was assessed by blastn70, after first removing low quality bases from metagenomic reads using Trimmomatic v0.3971 (parameters: LEADING:20 TRAILING:20 MINLEN:100). For each sediment trap metagenome, the total number of reads matching EpSB or EcSB nifH nucleotide sequences with 100% identity were tallied and normalized to the total number of reads in the database. Only EpSB-matching reads were detected in this analysis.Quantitative PCRSpecific PCR primers were designed targeting a 102 bp region of E. pelagica’s LSU gene (Epel-LSU-F, 5′-GAAACCAGTGCAAGCCAAC-3′; Epel-LSU-R, 5′-AGGCCATTATCATCCCTTGTC-3′) and an 85 bp region EpSB’s nifH gene (EpSB-nifH-F, 5′-CACACTAAAGCACAAACTACC-3′; EpSB-nifH-R, 5′-CAAGTAGTACTTCGTCTAGCTC-3′) and were synthesized by IDT. Gene copy concentrations were quantified for Station ALOHA water samples (~2 L) collected by Niskin bottles at 5, 25, 45, 75, 100, 125, 150, and 175 m on January 16 and July 1 (except 5 m), 2014, during HOT cruises #259 and #264. Samples were filtered onto 25 mm diameter, 0.02 μm pore size aluminum oxide filters (Anotop; Whatman, cat. # WHA68092102) and stored at −80 °C until extracting DNA using the MasterPure Complete DNA and RNA Purification Kit (Epicentre, cat. # MC85200) according to Mueller et al.72. Briefly, a 3-mL syringe filled with 1 mL of tissue and cell lysis solution (MasterPure) containing 100 μg mL−1 proteinase K was attached to the outlet of the filter, and the filter inlet was sealed with a second 3-mL syringe. The lysis solution was pulled halfway through to saturate the filter membrane, and the entire assembly was incubated at 65 °C for 15 min while attached to a rotisserie in a hybridization oven rotating at ca. 16 rpm. The lysis buffer was then drawn fully into the inlet syringe, transferred to a microcentrifuge tube, and placed on ice. The remaining steps for protein precipitation and removal and nucleic acid precipitation were carried out following the manufacturer’s instructions. For each sample, DNA was resuspended in a final volume of 100 μL. Quantitative PCR (qPCR) was performed using the PowerTrack SYBR Green Master Mix system (Applied Biosystems, cat. # A46109) and run on an Eppendorf Mastercycler epgradient S realplex2 real-time PCR machine. Reactions (20 µL total volume) were prepared according to the manufacturer’s protocol, containing 500 nM of each primer. Sample reactions (four replicates) contained 2 μL of environmental DNA extract (24–76 ng DNA), while standards (three replicates) contained 2 μL of gBlocks Gene Fragments (IDT) that were prepared at 1, 2, 3, 4, 5, and 6 log gene copies/μL. The gBlocks Gene Fragments were 500 bp in length and encompassed the entire E. pelagica UHM3201 LSU sequence and positions 1–500 of the EpSB UHM3201 nifH sequence, respectively. The main cycling conditions consisted of an initial denaturation and enzyme activation step of 95 °C for 2 min, followed by 40 cycles of 95 °C for 5 s and 57 °C or 55 °C for 30 s for the LSU and nifH genes, respectively. Melting curves were analyzed to verify the specificity of the amplifications, and reactions containing Epithemia catenata DNA extract were included as negative controls. Reaction efficiencies were 104.23% and 95.15% for the LSU and nifH genes, respectively. The limit of detection for these assays was not empirically determined. gBlocks sequences, qPCR threshold cycle values, and conversion equations are provided in the Source Data file.Physiology experimentsThe daily patterns of N2 fixation were quantified for E. pelagica UHM3200 and E. catenata UHM3210 using two techniques: acetylene (C2H2) reduction to ethylene (C2H4) and argon induced dihydrogen (H2) production (AIHP). Both analyses were conducted using a gaseous flow-through system that quantified the relevant trace gas on the sample outlet line with a temporal resolution of 10 min73. To conduct the measurements, a 10-mL subsample of each Epithemia culture was placed in a 20-mL borosilicate vial and closed using gas-tight rubber stoppers and crimp seals. Separate bottles were used for H2 production and C2H2 reduction. During the experimental period, the temperature was maintained at 25 ± 0.2 °C using a benchtop incubator (Incu-Shaker; Benchmark Scientific) and light exposure was 200 μmol photons m−2 s−1 at wavelengths of 380–780 nm with a 12:12 h square light:dark cycle (Prime HD+; Aqua Illumination). To conduct the AIHP method, the sample vial containing the culture was flushed with a high purity gas mixture consisting of argon (makeup gas; 80%), oxygen (20%), and carbon dioxide (0.04%). In the absence of N2, all of the electrons that would have been used to reduce N2 to NH3 are diverted to H2 production, thereby providing a measure of Total Nitrogenase Activity (TNA). The C2H2 reduction assay also represents a measure of TNA. Our analytical set-up introduced C2H2 at a 1% addition (vol/vol) to the high purity air with a total flow rate (13 mL min−1) identical to the AIHP method. The gas emissions were analyzed using separate reductive trace gas analyzers that were optimized for the quantification of H2 and C2H4. To verify the observed daily patterns in N2 fixation, 15N2 assimilation measurements were conducted on triplicate samples of Epithemia cultures at targeted time points. Five milliliters of 15N-enriched seawater was added to the subsamples, which were subsequently crimp sealed and incubated for a 2 h period with the same light and temperature conditions as the daily gas measurements. At the end of the incubation, the contents of each vial were filtered onto a pre-combusted glass fiber filter. The concentration and isotopic composition (δ15N) of particulate nitrogen for incubated and non-incubated (i.e., natural abundance) samples was measured using an elemental analyzer/isotope ratio mass spectrometer (Carlo-Erba EA NC2500 coupled with a ThermoFinnigan Delta Plus XP). For each of the described analyses, cell-specific rates were calculated based on the average of triplicate cell concentration measurements, obtained from cell samples preserved at 4 °C with Lugol’s iodine solution and quantified within a week using a Sedgwick-Rafter counting chamber (Electron Microscopy Sciences, cat. # 68050-52). All rate measurement data is provided in the Source Data file.Reporting summaryFurther information on research design is available in the Nature Research Reporting Summary linked to this article. More

All of the statistical analyses were conducted using the R programming language version 4.0.5 within the RStudio IDE version 1.4.124,25. Data visualization and processing were performed with the ‘tidyverse’ collection, ‘foreach’ and ‘doParallel’ R packages26,27,28. Geographical Information System (GIS) operations on raster and vector files were conducted using the ‘sf’, ‘exactextractr’ and ‘raster’ R packages29,30,31.Data sources and pre-processingBiodiversity dataWe used the North American Breeding Bird Survey dataset as our source of biodiversity data due to its long temporal coverage and spatial extent14,32. The BBS is composed of bird species abundance records collected since 1966 from over 4,000 survey routes across the countries of Mexico, USA and Canada. For this study we focused solely on routes in the USA, due to their longer time dimension. Data collection follows public access roads that are 24.5 miles long (circa 39.2 km) using a point count protocol whereby routes are surveyed every half-mile (800 m) for a total of 50 stops. At each stop, observers stand for 3 min and record the species and the abundance of every bird seen or heard within 400 m of their location. The routes are surveyed by volunteers with experience in bird observation, and surveys are conducted from late April to July to capture the peak of the breeding season.We selected the years 2001 and 2016 as the two timepoints of our analysis. This 15-year timeframe corresponded to the longest possible timespan for which land cover data products were available at high spatial resolution18. Before analysis, we subset the BBS dataset by removing routes that had incomplete survey lengths (less than 50 point count stops, indicated by the RouteTypeDetailID field value being less than 2 in the extracted BBS dataset) or that were surveyed under adverse weather conditions such as high wind and rain (as indicated by the Run Protocol ID field being equal to 1), which could affect bird occurrence and detectability. Following this filtering process, the total number of BBS routes analysed was 960 (Extended Data Fig. 1).For higher precision when inferring the relationship between avian diversity and environmental variables, we subdivided each route into five segments of equal length, consisting of 10 count locations each. This approach was motivated by the need to more closely associate bird communities with the land cover composition in the area in which they are found. To minimize the spatial autocorrelation between adjacent segments and avoid overlaps in landscapes analysed, we filtered the data to keep only the first, third and fifth segment of each route. These segments therefore formed our sampling unit used in all analyses.We recognized that environmental conditions and stochastic trends in populations could introduce variability in biodiversity calculated from bird community data. We therefore extracted, for each segment and each species, the average population count across a 3-year period centred on our two timepoints (2000, 2001, 2002 and 2015, 2016, 2017)33. We then calculated the mean abundance of each species across these 3 years.The effect of observer experience34,35,36 was accounted for by sourcing the observer ID responsible for each route at each timepoint and including it as a random effect in the legacy model (see ‘Model development’ section). We also controlled for the time of day as it is plausible to expect visibility and avian species activity patterns to vary between early morning and later parts of the day. Time of day for each segment was calculated by averaging across the start and end time data entries associated with each route, and then including this as a covariate in both the legacy and equilibrium models (see ‘Model development’ section). However, we did not model detectability issues associated with traffic noise and disturbance for two reasons. First, all BBS surveys are conducted along public access roads with a vehicle, so the disturbance is expected to be similar across sites. Second, previous studies have found no clear evidence for noise being the main cause for reduced bird abundance near roads37.Following these procedures, our processed BBS dataset included entries of mean abundances of each species for a total of 2,880 segments, corresponding to segment 1, 3 and 5 of 960 routes (Extended Data Figs. 1 and 2). For each segment, at each timepoint we calculated different measures of alpha diversity following the Hill numbers framework38. We then selected to use the effective number of species at q = 1, calculated as the exponential of the Shannon–Wiener Index38. The effective number of species at q = 1 sits at the theoretical half-way point between the classic species richness measure that accounts only for the absolute number of species (q = 0) and the Berger-Parker dominance index (q = infinity), which instead only reflects the most common species. Thus, the effective number of species is a robust alternative to species richness, which does not take account of species rarity or detectability and can thus lead to biased biodiversity estimates16,17.Land cover and environmental dataLand cover data for the US for our focal years of 2001 and 2016 were sourced from the open-access NLCD CONUS products developed by the US Geological Survey (USGS)18,39. The NLCD products are high-resolution (30 m pixel dimensions) classified raster files covering the land area of the whole USA. This dataset provides us with the opportunity to look at finely gridded spatio-temporal changes in a landscape over a relatively long timeframe of 15 years, while utilizing data collected and analysed with the same methods (for example, land use classification algorithms).To reduce the number of potentially collinear explanatory variables included in our models, we aggregated the land cover variables provided by the NLCD dataset. We summarized these to five land cover categories: ‘urban’ (an aggregate of the Developed-Open Space (subclass 21), Developed-Low Intensity (22), Developed-Medium Intensity (23) and Developed-High Intensity classes); ‘forest’ (an aggregate of the Deciduous Forest (41), Evergreen Forest (42) and Mixed Forest (43) classes); ‘grassland’ (an aggregate of the Shrub (52), Grassland/Herbaceous (71) and Pasture/Hay (81) classes); ‘cropland’ (cultivated Crops (82) subclass) and ‘wetland’ (an aggregate of the Woody Wetland (90) and Herbaceous Wetland (95) classes). The Perennial Ice/Snow (12), Open Water (11) and Barren Land (31) classes were excluded from the analysis as they were very uncommon in our dataset. The distribution and total area of the land cover categories across the US are shown in Supplementary Figs. 1 and 2. Temperature data were sourced from the 30 arc-seconds gridded PRISM climate database19 and were extracted as the mean across May and June for each group of years from which bird abundances were taken.We first sampled the landscape surrounding each segment using a range of buffer shapes and sizes, and then selected the buffer type on the basis of the capacity of each buffer type to explain the response variable. The types of buffers that we explored were: a circular buffer around the centroid of the polygon defined by the vertices of each segment (4,000 m radius) and a series of three buffers around the segment line (500 m, 2,000 m and 4,000 m radius). The best fit was given by the smallest buffer size of 500 m, shown in Extended Data Fig. 2, which also coincides with the BBS protocol effective counting distance of 400 m and more closely reflects the size of bird territories14. Land cover variables were computed as a percentage of the total buffer area. Change in percentage points for each land cover type between the 2 years was computed by subtracting the values at the two timepoints. A change product is also provided by the USGS databases40, but it does not meet our needs because it considers land cover changes based on a ranking. Nonetheless, a comparison of urban land cover change between the timepoints showed a similar result (Supplementary Fig. 4). Land cover data were processed geospatially using the NAD 83 Conus Albers Coordinate Reference Systems projection, EPSG 5070.Model developmentTheoretical backgroundWe developed a statistical model that conceptualized extinction debts and colonization credits by combining the following two concepts: (1) the settled biodiversity of avian communities in a given landscape composition (that is, a system at equilibrium) and (2) the lagged response in the species diversity in a given landscape due to recent land cover changes (that is, a system moving to a new equilibrium). We reasoned that, given enough time, and with no further changes in land cover, the effective number of species at a given location would eventually equilibrate. The equilibrium distribution of the effective number of species emerges with the waning of the legacy effect of previous landscape compositions in encouraging or impeding the recruitment and survival of particular species. We did not model these ecological mechanisms directly, but instead expressed the equilibrium of the effective number of species, and the rate of approach to this equilibrium, as empirical functions of environmental covariates. It is important to keep in mind that during a finite time interval following environmental change, it is possible that our observations of effective number of species represent a system in a transitory state towards its new equilibrium. Yet, environmental changes may occur at rates that never allow the system to equilibrate. Although the equilibration processes are latent (that is, not amenable to direct observation), the combination of equilibrium and temporal legacy components into an integrated model, applied to a dataset with extensive environmental replication (due to spatial expansiveness), has allowed us to retrieve distributions for all relevant model parameters (see below).Model overviewThe observed effective number of species Rs,t at site s in year t for t = t1,t2 is modelled as a normally distributed variate with mean μs,t and standard deviation σ$$R_{s,t} approx mathrm{Normal}left( {mu _{s,t},sigma } right)$$
(1)
We assume that, under landscape change, the system is in a state of flux and that the data are from observations witnessing the transition between two (unattained) equilibria. The expected state of the system at any given point in time, μs,t, was formulated as a mixture of past and future equilibrium distributions (that is, a weighted average of the two distributions, where the weights are given by the complementary proportions ω and 1 − ω)$$mu _{s,t} = fleft( {x_{s,t_2};beta } right)omega left( {{Delta}x_{s,t_1,t_2};gamma } right) + fleft( {x_{s,t_1};beta } right)left( {1 – omega left( {{Delta}x_{s,t_1,t_2};gamma } right)} right)$$
(2)
Here, the function f describes the equilibrium distribution of the effective number of species as a function of the configuration of the local environment, captured in covariates xs,t. The weighting function ω depends on covariates ys,t derived from the difference in the local land cover between 2016 and 2001 (that is, it is a function of the land cover change that has taken place). The mixture weights ω and (1 − ω) determine the relative importance of the two equilibrium distributions (past or current). If ω = 1, the interpretation is that the new equilibrium distribution has been completely attained, and thus the current (2016) effective number of species is entirely explained by the current (2016) land cover. Conversely, if ω = 0, the current effective number of species is entirely explained by the past (2001) land cover. The vectors of parameters β and γ, presented in equation (2), are inferred from model fitting.We also augmented equation (2) with a function g of static covariates and random effects z that we expect to have an impact on the effective number of species. Thus, the model comprised equilibrium components, a temporal legacy component and static covariates:$$mu _{s,t} = fleft( {x_{s,t_2};beta } right)omega left( {{Delta}x_{s,t_1,t_2};gamma } right) + fleft( {x_{s,t_1};beta } right)left( {1 – omega left( {{Delta}x_{s,t_1,t_2};gamma } right)} right) + gleft( {z_s;alpha } right)$$
(3)
in which (fleft( {x_{s,t};beta } right)) are the equilibrium components for the two timepoints, (omega left( {Delta x_{s,t_1,t_2};gamma } right)) is the temporal legacy component, and (gleft( {z_s;alpha } right)) is the function that captures the static covariates and random effects, with α being the estimated static covariates parameter effects.Equilibrium componentsWe defined the equilibrium distribution of the effective number of species at a given timepoint as a function (fleft( {x_{s,t};beta } right)) of land cover. This function describes the expected effective number of species at location s, given sufficient time for the community to adapt to the given land cover composition. We now describe this function in more detail.The equilibrium component was formulated as a log-linear model comprising a total of I = 5 environmental covariates (the percentage cover of five landscape classes: urban, forest, grassland, wetland and cropland), using 2nd-order polynomial terms, captured by the coefficient j, to account for optima in effective number of species along each of the five environmental dimensions:$$fleft( {x_{s,t}} right) = {mathrm{exp}}left( {beta _0 + mathop {sum }limits_{i = 1}^{I = 5} mathop {sum }limits_{j = 1}^{J = 2} beta _{i,j}x_{i,s,t}^j} right)$$
(4)
In equation (4), the β parameters capture the effect of covariates on the equilibrium and are assumed to be the same for each environmental composition. A simplifying assumption necessary for the application of this model is that the effective number of species had equilibrated at the first timepoint. As data become available for more years in the future, the influence of this assumption on the model results will diminish and more accuracy will be achievable with multiple timepoints.To allow for conditionality in the effects of one land cover variable on the response of the effective number of species to another land cover variable, we extended this function with pairwise interaction terms k between all the linear terms for land cover variables and pairwise linear-quadratic terms, as follows:$$fleft( {x_{s,t}} right) = {mathrm{exp}}left( {beta _0 + mathop {sum }limits_{i = 1}^{I = 5} mathop {sum }limits_{j = 1}^{J = 2} beta _{0,i,j,}x_{i,s,t}^j + mathop {sum }limits_{i = 1}^{I = 4} mathop {sum }limits_{k = i + 1}^{K = 5} beta _{1,i,k}x_{i,s,t}x_{k,s,t}} right)$$
(5)
Temporal legacy componentThe main covariates, ({Delta}x_{i,s,t_1,t_2}), for the part of the model that captures the temporal legacy, (omega left( {{Delta}x_{s,t_1,t_2};gamma } right)), are derived from the change in land cover (({Delta}x_{i,s} = x_{i,s,t_2} – x_{i,s,t_1})) between the two timepoints$$x_{i,s,} = left{ {begin{array}{*{20}{l}} {x_{1,i,s} = left| {{Delta}x_{i,s}} right|,} hfill & {x_{2,i,s} = 0,} hfill & {{mathrm{if}},{Delta}x_{i,s} < 0} hfill \ {x_{1,i,s} = 0,} hfill & {x_{2,i,s} = {Delta}x_{i,s},} hfill & {{mathrm{otherwise}}} hfill end{array}} right.$$ (6) where ({Delta}x_{s,t,z}) is a vector of the ith environmental change variable (that is, urban, forest, grassland, wetland, cropland) at site s and for directionality z. The effect of these covariates on the mixture weights is given by:$$omega left( {{Delta}x_{s,t_1,t_2};gamma } right) = {mathrm{exp}}left( {mathop {sum }limits_{i = 1}^{I = 5} - gamma _{i,z}{Delta}x_{z,s,i}} right)$$ (7) This formulation weights the contribution that the environmental variables at the two timepoints have on the current effective number of species, as a function of the magnitude and directionality of change in each type of land cover covariate. The γ parameters, and subsequently the temporal legacy component, are allowed via the inclusion of the environmental change data ({Delta}x_{s,t,z}), to account for the distance between the land cover at the two timepoints, therefore quantifying how far the initial community would need to travel to reach equilibrium in 2016 as a function of the type, magnitude and directionality of change. It should be noted that our model, in equation (3), is only implicitly related to the speed with which the effective number of species reacts to environmental changes. Instead, it quantifies how much further it would still have to travel to reach the expected equilibrium associated with the current configuration of the landscape.Static covariatesAs described in model equation (3), we included a function of static covariates to which we can expect the effective number of species to respond without lags relating to the past landscape. We added a linear and quadratic fixed effect for temperature in 2016 to control for any difference in the effective number of species related to climatic characteristics and to allow for a parabolic relationship to be expressed (optima either at mean or extremes values). We also controlled for the heterogeneity of a landscape by including the effective number of land cover types, computed in the same way as the effective number of species, as a fixed effect40. A fixed effect for time of day, reflecting the time at which each segment was surveyed, was included to correct for differences in species detectability between early morning and later parts of the day41. An observer-level random effect was also added to control for variation between observers35,36 and partly account for between-route variation, given that we would expect observers who collect data from multiple routes to do so within a relatively small area. Spatial autocorrelation of the effective number of species was tested for all segments at once and by different radiuses for neighbour inclusion (500 m, 1,000 m, 5,000 m, 10,000 m, 100,000 m), using the Moran’s I statistic42. Spatial autocorrelation was not corrected for because Moran’s I was not significant at any spatial scale (P  > 0.05). Pseudo-replication between neighbouring segments was avoided by considering segments 1, 3 and 5, whose land cover buffers did not overlap (Extended Data Fig. 2).Model fittingThe model was fitted within a Bayesian framework using a Hamiltonian Markov chain Monte Carlo algorithm implemented in the STAN programming language43 version 4.3.0 and the ‘cmdstanr’ R package version 2.26.144.We ran 4 chains, sampling for 1,000 iterations with a burn-in period of 500 iterations each. These numbers of iterations were sufficient to achieve chain convergence. The STAN sampling was run on four parallel threads on a multi-core Intel i7 – 8750H processor with a maximum clock speed of 4.1 GHz.For the purposes of Bayesian inference, all slope parameters associated with the equilibrium component equation (5) and the static additive terms were assigned an unbiased prior (beta _{i,j} approx Nleft( {0,1} right),{mathrm{and}},z_s approx Nleft( {0,1} right),) where N is normal, with the aim of shrinking the parameter estimated towards 0 (that is, no covariate effect). A gamma distributed prior, with shape and rate 0.001, was assigned to the standard deviation of the random effect. For the following known and expected relationships, we also truncated the range of parameter values by bounding the upper or lower limits of the prior/posterior distributions. Intercept and standard deviation of the observer random effect were bounded below by 0. Linear effects for the environmental covariates and temperature were bounded below at 0, while their quadratic counterparts were bounded above at 0. Interaction terms were not limited. The temporal legacy component parameters were given a uniform (U) prior (gamma _isim Uleft( {0,1} right)), bounded between 0 and 1 to act as a weighting proportion between the present and the past. The upper bound on the gamma parameters to 1 does not bias us towards an increased contribution of the past land cover, but instead provides a more conservative approach.Model diagnostics were conducted by assessing chain convergence visually through trace plots, as well as statistically by employing the Gelman-Rubin test, which compares the estimated between-chain and within-chain variances45. Chain autocorrelation and the associated effective sample size were also monitored. In the case of low effective sample size, the chains were extended until the effective sample size exceeded a threshold value of 400. The marginal posterior distribution for each parameter was visualized via a density plot to check for multimodality.Model selection was conducted to inform choice of the size and shape of the land cover buffer around each sampled segment. We did so by comparing values of the Watanabe-Akaike Information Criterion leave-one-out (WAIC)-loo information criterion46 of four different models, each computed using land cover data calculated with two different buffer options of various sizes: a circular buffer around the centroid of the polygon defined by the vertices of each segment (4,000 m radius) and a series of buffers around the segment line (500 m, 2,000 m and 4,000 m radius). This approach was implemented through the ‘loo’ R package version 2.1, which provides an improvement on the original WAIC by including diagnostic measures around the point-wise log-likelihood value estimated around each sample draw47.Visualization of model predictionsA map of the USA (Fig. 1) was produced to represent the predicted extinction debts and colonization credits (that is, positive or negative distance in the effective number of species from the expected equilibria). The map was produced on a hexagonal grid at a spatial resolution of 10 km vertex-to-opposite-vertex, with each hexagon covering a total of 86 km2. Values of extinction debt and colonization credit were calculated by subtracting the predicted effective number of species produced by the model (equation 3) from the predicted effective number of species at equilibrium in 2016 (that is, when the legacy component equals 1). To correctly propagate and represent uncertainty in the extinction debts and colonization credits presented, this process was repeated 1,000 times for predictions originating from different draws from the posterior distribution. Uncertainty in the form of the geometric coefficient of variation, calculated as (root {2} of {{e^{left( {mathrm{log}left( {sigma + 1} right)^2} right)} – 1}}) where σ is the standard deviation, is mapped in Extended Data Fig. 4a. Extended Data Fig. 4 also includes a copy of Fig. 1 (Extended Data Fig. 4b) for reference, alongside upper (Extended Data Fig. 4c) and lower (Extended Data Fig. 4d) credible intervals.Over/underestimation values of biodiversity that could arise by neglecting debts and credits were computed as the difference between the effective numbers of species predicted by the equilibrium model and the legacy model, multiplied by 100 and then divided by the predicted effective number of species under the legacy model. This calculation results in a percentage measurement of the extent to which (in relative terms) the current effective number of species under- or overestimates the diversity that a given landscape can sustain at equilibrium.To further validate our predicted extinction debts and colonization credits, we compared the direction of the expected changes with the recorded difference in effective numbers of species between 2016 and 2019 (the latest year for which data are available). To do so, we sourced bird abundances from the North American BBS dataset14,32 for the year 2019 and conducted the same data processing as described above for the other two timepoints. We then conducted a Pearson correlation test to assess how well the observed change followed the model-predicted one. We are nevertheless aware that a 3-year timeframe is unlikely to be large enough for debts and credits to fully manifest.Plots were also generated to describe the behaviour of the mixture weight, ω (equation 7), which captures the contribution (weighting) of the landscape composition in determining the effective numbers of species at the two timepoints (Fig. 2 in the main text). Values of ω across the whole spectrum of plausible land cover change values (that is, from −100 to +100) were simulated by averaging over 10,000 draws from the posterior distribution of each γ parameter. Credible intervals were measured by taking the 95% range of the 10,000 draws.Explaining spatial variation in debts and creditsThe extinction debts and colonization credits predicted for the contiguous USA states were further modelled to identify which past land cover changes were the main drivers of the delayed biodiversity changes in USA bird communities. We considered the values of debts or credits associated with the 92,000 individual 86 km2 hexagons (Fig. 1) as a response variable. We then specified a Gaussian linear model including the magnitude of each land cover change as explanatory covariates. Positive and negative changes in each covariate were treated as separate linear components to differentiate their effects. The model was fitted to 1,000 sets of debts and credits, each originating from predictions based on independent draws from the posterior distribution. For each generalized linear model (GLM) fit, we then subsequently sampled each parameter distribution another 1,000 times and extracted the summarized parameter estimates from a total of 100,000 values. Model coefficients and their resulting uncertainty from the above process are presented in Fig. 4 and in more detail as part of Supplementary Table 3.Reporting SummaryFurther information on research design is available in the Nature Research Reporting Summary linked to this article. More

Pollen-inferred landscape change and pre-industrial demographyRecently, data derived from tree rings or ice cores have been employed to approximate changes in human economic activity related to past epidemics, as well as to warfare and climatic variability46,47. However, none of these proxies is directly related to human demography or provides a basis to estimate variation in the Black Death’s mortality on a regional scale across Europe (to date only a single archaeological study using pottery as a proxy for demographic change on the national level, focusing on just a single country—England—has appeared48).In recent years, pollen data have been proven to be closely related to demographic variability. Most importantly, detailed comparisons of historical documentary data on population trends and landscape changes as revealed by pollen data have been carried out on a local scale and a close link between changes in European pollen data and changes in European local demography over the past millennium has been demonstrated on multiple occasions, that is, during the period and region of our concern here49,50. A strong link between long-term demographic trends as visible in regional settlement numbers and macro-changes in land cover (deforestation/afforestation) have also been confirmed for ancient Greece51. Additionally, a recent publication successfully employed pollen data to test the extent of the mortality associated with the sixteenth-century Spanish and Portuguese empires’ colonization of tropical regions in the Americas and Asia52. However, as of now there is no method to quantify past demographic trends in absolute numbers based on palaeoecological data. Consequently, we also focus in this paper on relative changes in historical societies’ populations and test the now common idea that the Black Death caused enormous mortality across Europe (with many scholars now arguing for a mortality exceeding 30% and upwards of 50% of the population within a few years) (see also Fig. 1). Using our BDP approach, we conclude this hypothesis is not maintainable. Our evidence for demography-related landscape changes (or lack thereof) negates it.Our main indicator is cereal pollen. In pre-industrial economies, rural labour availability (hence rural population levels) and the spatial scale of cereal cultivation were directly related. An increase in the extent and intensity of cereal cultivation—as reflected in pollen data—would have required not only a predilection and demand for cereals, but also greater availability of labour and thus population growth or significant immigration. The maintenance of existing agricultural activity, in turn, would have required relatively stable population levels53,54,55. The uniform ~50% mortality postulated for the Black Death across Europe should have resulted in a large and significant decline of cereal cultivation and parallel forest regrowth across Europe, as previously demonstrated for mid-fourteenth-century Sweden26 and singular sites in some regions of western Europe56. This result agrees with the fact that Black Death mortality could be high among people at productive age, as illustrated for England57,58. Moreover, even in the case of England, a comparatively commercialized and adaptive rural economy in mid-fourteenth-century Europe, the loss of 50% of the population led to a significant decline in the total area under cultivation (as documented by heterogeneous written sources)59. In Italy, another well-developed economy at that time, the expansion of large estates following the Black Death also did not compensate for the general loss of cereal productivity60. This effect, high mortality driving arable contraction, must have been yet more pronounced in more subsistence-oriented and less adaptive economies, with limited surplus production, such as in regions of the Iberian Peninsula, Germany, Sweden and particularly east-central Europe. Importantly, palaeoecological evidence for arable contraction may be indicative, to some extent, of not only rural population decline but also urban population decline in the region, as there is evidence in some areas, following the pandemic, of rural-to-urban migration, of country-dwellers repopulating urban centres10. Possibly less common was intraregional rural migration, as marginal lands were abandoned for better quality soils, which were more likely to remain under cultivation26,61.Therefore, cereal pollen remains our most potent pollen indicator related to demographic changes in pre-industrial European societies. Other pollen indicators, reflecting rewilding and reforestation (secondary ecological succession) of cereal fields abandoned as a result of significant mortality, or the transformation of cereal fields into pastures, which required less rural labour and thus also could have been a response to high plague mortality, play a secondary role in our analysis and provide further support for our conclusions.BDP data collectionExisting online palynological databases (the European Pollen Database (EPD)62 www.europeanpollendatabase.net, and the Czech Quaternary Palynological Database (PALYCZ)63, https://botany.natur.cuni.cz/palycz/), as well as personal contacts of the study authors and a systematic publication search were employed to identify palynological sites in Europe reaching the required chronological and resolution quality for the study of the last millennium. In order to enable statistical analysis, we included only sites clustered in well-defined historical-geographical regions, excluding isolated sites even if the quality of a site’s data was very good. Data of sufficient quality and amount from regions for which the Black Death is well-studied, notably central and northern England and the Low Countries, is not presently available; to the best of our knowledge, for each of these regions there currently is not more than a single isolated site56, which does not allow for the application of statistical approaches.In total, 261 pollen records with the average temporal resolution of 58 years and 14C-age control (or varve chronology), have been collected. The age–depth models of the sequences have been provided by authors in original publications, by the EPD or developed through the Clam package (version 2.3.4) of R software for the purpose of this study. The analytical protocol for pollen extraction and identification is reported in the original publications. The Pollen Sum includes all the terrestrial taxa with some exceptions based on the selection done in the original publications. The full list of sequences, exclusions from the Pollen Sum, age-depth models and full references are reported in Supplementary Data 1.The taxa list has been normalized by applying the EPD nomenclature. In this respect, the general name Cichorioideae includes Asteraceae subf. Cichorioideae of the EPD and PALYCZ nomenclatures, which primarily refers to the fenestrate pollen of the Cichorieae tribe64. Ericaceae groups Arbutus unedo, Calluna vulgaris, Vaccinium and different Erica pollen types, whereas deciduous Quercus comprehends both Q. robur and Q. cerris pollen types65. Rosaceae refers to both tree and herb species of the family. Finally, Rumex includes R. acetosa type, R. acetosella, R. crispus type, Rumex/Oxyria and Urtica groups U. dioica type and U. pilulifera.BDP summary pollen indicatorsIn order to connect changes visible in the pollen data to human demographic trajectories, we assembled four summary pollen indicators that describe specific landscapes related to human activity. They reflect different degrees of demographic pressure on the landscape (cereal cultivation, pastoral activities, which are less-labour intensive than cereal cultivation, abandonment and rewilding) as well as different durations of land abandonment that might have occurred post-Black Death. Our indicators account for the fact that Europe is a continent rich in natural heritage, with a wide range of landscapes and habitats and a remarkable wealth of flora and fauna, shaped by climate, geomorphology and human activity. In order to ensure uniform interpretation of the indicators, we relied on criteria that can be applied to all European landscapes regardless of their local specificity. Cereals and herding are directly related to human activities and are barely influenced by spatial differences. More complex is the succession of natural plants with their ecological behaviour and inter-species competition. For this reason, we relied on existing quantitative indicators of plant ecology.The Ellenberg L – light availability indicator66 provides a measure of sunlight availability in woodlands and consequently of tree-canopy thickness, reflecting the scale of the natural regeneration of woodland vegetation after cultivation or pasture activities61. Nonetheless, ecological studies have suggested that geographic and climatic variability between different European regions can influence the Ellenberg indicator system67,68,69,70,71. The original indicators were primarily designed for Central Europe58, but several studies developed Ellenberg indicators for other regions, reflecting the specific ecology of the selected taxa (British Isles;72 Czech Republic;73 Greece;74 Italy;75 Sweden76). Plants with L values between 5 and 8 are listed in the fast succession indicator, the ones with L values ranging from 1 to 4 are included in the slow succession indicator. The result is the following list:1) Cereals: only cultivated cereals have been included: Avena/Triticum type, Cerealia type, Hordeum type, Secale. 2) Herding includes pastoral indicators linked to the redistribution of human pressure: Artemisia, Cichorioideae, Plantago lanceolata type, Plantago major/media type, Polygonum aviculare type, Rumex, Trifolium type, Urtica, Vicia type. 3) Fast Succession comprises indicators of relatively recent reforestation of cultivated land after abandonment: Alnus, Betula, Corylus, Ericaceae, Fraxinus ornus, Juniperus, Picea, Pinus, Populus, deciduous Quercus, Rosaceae. 4) Slow Succession includes indicators of secondary succession established after several decades of abandonment: Abies, Carpinus betulus, Fagus, Fraxinus, Ostrya/Carpinus orientalis, Quercus ilex type.In order to validate the indicators overcoming the regional limits of Ellenberg values, a different subdivision has been provided following the Niinemets and Valladares shade tolerance scale for woody species of the Northern Hemisphere77. The subdivision of taxa in the Fast and Slow succession indicators remains the same with only three changes: Fraxinus ornus and Picea move from Fast to Slow succession and Fraxinus from Slow to Fast succession. Extended Data Figs. 1 and 2 show that the two groupings yield the same results, which confirms the reliability of our indicators. There is only one clear exception (Russia), with one more region where smaller-scale diversion occurs for only one indicator, Slow Succession (Norway). The different indicator behaviour results from the different attribution of Picea in our two sets of succession indicators: at high latitude, Picea characterizes the final stage of the ecological succession and hence its different attribution results in different summary indicator values in Russia for the two stages of ecological succession, fast and slow.Please note our summary indicators are not designed to reflect the entirety of the landscape and reconstruct all of its different components. Rather, they are a means of approximating changes in the landscape related to the types of human activities, and their intensity, as much as they relate to demographic changes in human populations using and inhabiting these landscapes.BDP analytical statistical and spatial methodsTo control for local specificity, pollen percentages of every taxon from each pollen site were standardized. From the taxa percentage in a given year the arithmetic mean calculated for the observations from the period 1250–1450 was subtracted and the result divided by the standard deviation for the 1250–1450 period. Standardized taxa results were assembled for each site into four BDP summary indicators. Since each indicator has different numbers of taxa, the sum of standardized taxa values calculated for a given year and site was divided by the number of taxa in the indicator. For the purposes of replication, this standardized pollen dataset, comprising the four indicators for each sample from each site, is available as Supplementary Data 2.This dataset has been analysed in two ways, statistically and spatially.For the statistical approach, standardized regional indices of landscape transformation were created for each region by calculating the average value for all sites within the region, for each of the subperiods analysed in the study (1250–1350 and 1351–1450; 1301–1350 and 1351–1400; 1325–1350 and 1351–1375). Differences between means for each subperiod were measured by the use of the bootstrapping based on 10,000 resamples. The 90% and 95% confidence intervals were estimated with the bias-corrected and accelerated method (BCa)78. These results are visualized in Fig. 5 for the comparison of the subperiods of 1250–1350 versus 1351–1450, and in Supplementary Figs. 4 and 6 for the comparison of the subperiods of 1300–1350 versus 1351–1400 and 1325–1350 versus 1351–1375, respectively.For the spatial approach, we employed the Bayesian model AverageR developed within the Pandora and IsoMemo initiatives (https://pandoraapp.earth/) to map the distribution of pollen indices across Europe. AverageR is a generalized additive model that has been described previously79. It relies on a thin-plate regression spline80 to predict new, unseen data using the following model:$$Y_{i} = {g}( {{{{mathrm{longitude}}}},{{{mathrm{latitude}}}}} ) + {varepsilon}_{i}$$Where Yi is the independent variable for site i; g(longitude, latitude) is the spline smoother; and εi ∼ N(0, σε) is the error term.The spline smoother can be written as X × β where X is a fixed design matrix and β is the parameter vector. Surface smoothing is controlled by employing a Bayesian smoothing parameter estimated from the data and trades-off bias against variance to make the optimal prediction75. This parameter β is assumed to follow a normal distribution: β ~ N(0, 1 /δ × λ × P), where P is a so-called penalty matrix of the thin plate regression spline, which penalizes second derivatives81. The δ parameter is by default set to 1 but this can be adjusted to suit smoothing needs for each application. In our study δ was set at 0.9 to match the preferred spatial scale of analysis for our dataset (approx. 250 to 500 km).AverageR was employed to generate smoothed surfaces for three sets of temporal bins (1250–1350 versus 1351–1450, as well as 1300–1350 versus 1351–1400 and 1325–1350 versus 1351–1375) and for the four BDP indicators (Supplementary Figs. 3, 5 and 7). For the same indicator the difference between the two temporal bins was plotted (Fig. 5; Supplementary Figs. 4 and 6).Reporting SummaryFurther information on research design is available in the Nature Research Reporting Summary linked to this article. More

Erwin, D. H. Curr. Biol. 25, R930–R940 (2015).CAS 
Article 

Google Scholar 
Rabosky, D. L. & Hurlbert, A. H. Am. Nat. 185, 572–583 (2015).Article 

Google Scholar 
Hughes, M., Gerber, S. & Wills, M. A. Proc. Natl Acad. Sci. USA 110, 13875–13879 (2013).CAS 
Article 

Google Scholar 
Novack-Gottshall, P. M. et al. Nat. Ecol. Evol. https://doi.org/10.1038/s41559-021-01656-0 (2022).Article 

Google Scholar 
Wright, D. F. Sci. Rep. 7, 13845 (2015).
Google Scholar 
Deline, B. in Elements of Paleontology hcxz (Cambridge Univ. Press, 2021).Ausich, W. I., Kammer, T. W., Rhenberg, E. C. & Wright, D. F. Palaeontology 58, 937–952 (2015).Article 

Google Scholar 
Zamora, S. & Rahman, I. A. Palaeontology 57, 1105–1119 (2015).Article 

Google Scholar 
Sheffield, S. L. & Sumrall, C. D. Palaeontology 62, 163–173 (2018).Article 

Google Scholar 
Deline, B. et al. Curr. Biol. 30, 1672–1679 (2020).CAS 
Article 

Google Scholar 
Novack-Gottshall, P. M. Paleobiology 33, 273–294 (2007).Article 

Google Scholar 
Anderson, P. S. Paleobiology. 35, 321–342 (2009).Article 

Google Scholar 
Cole, S. R. & Hopkins, M. J. Sci. Adv. 7, eabf4027 (2021).Article 

Google Scholar 
Hopkins, M. J. & St. John, K. Syst. Biol. 70, 1163–1180 (2021).Article 

Google Scholar 
Blomberg, S. P., Garland, T. Jr & Ives, A. R. Evolution 57, 717–745 (2003).Article 

Google Scholar 
Wright, A., Wagner, P. J., & Wright, D. F. in Elements of Paleontology hcx2 (Cambridge Univ. Press, 2021).Cole, S. R., Wright, D. F. & Ausich, W. I. Palaeogeogr. Palaeoclimatol. Palaeoecol. 521, 82–98 (2019).Article 

Google Scholar  More

Isolation and cultivationIn total, 39 mixotrophic flagellates were investigated (Table 1). Thirty-three isolates were enriched and isolated from euphotic zone samples at Station ALOHA (22° 45′ N, 158° 00′ W) in February and May 2019. To select for mixotrophic grazers of Prochlorococcus, whole seawater was amended with K medium [28] (1/20 final concentration), and live Procholorococcus (MIT9301) was added as prey (~5 × 106 cells mL−1 final concentration). Enriched seawater samples were incubated under ~70 µmol photons m−2 s−1 irradiance on a 12 h:12 h light:dark cycle and monitored by microscopy daily up to five days. Each day, samples were serially diluted to extinction (9 dilution steps, 12 replicates per dilution) in 96-well plates in nutrient-reduced K medium (1/20 concentration) with a constant background of Prochlorococcus cells. Wells at the highest dilution showing growth of putative grazers were subjected to 3–6 further rounds of dilution to extinction. Four additional mixotrophs were isolated in full K medium using water from earlier cruises, and two (dictyochophyte strains UHM3021 described in [24] and UHM3050) were enriched in K minus nitrogen medium (K-N) without Prochlorococcus enrichment. All isolates were rendered unialgal, but not axenic, and maintained at 24 °C in K-N medium (~0.2 µM N) amended with Prochlorococcus prey, under the same light conditions as above. Dense Prochlorococcus cells grown in Pro99 medium [29], were harvested and concentrated through gentle centrifugation at 2000 RCF for 5 minutes and resuspended in fresh K-N medium to minimize nutrient carryover. To ensure their long-term accessibility, the isolates used in this study are being transferred to the National Center for Marine Algae and Microbiota at the Bigelow Laboratory for Ocean Science, East Boothbay, ME, USA (ncma.bigelow.org).Table 1 Mixotrophic flagellates investigated in this study.Full size table18S rRNA gene sequencing and phylogenetic analysisCells were harvested by centrifuging 25–50 mL dense cultures at 3000 RCF for 10 min at 4 °C. Genomic DNA was extracted from the pellets using the ZymoBIOMICS DNA Kit (Zymo Research, Irvine, CA, USA). A near-full-length section of the eukaryotic small-subunit ribosomal RNA (18S rRNA) gene was amplified by PCR with the Roche Expand High Fidelity PCR System (Sigma-Aldrich, St. Louis, MO, USA) using either forward primer 5′-ACCTGGTTGATCCTGCCAG-3′ and reverse primer 5′-TGATCCTTCYGCAGGTTCAC-3′ [30], or Euk63F 5′-ACGCTTGTCTCAAAGATTA-3 and Euk1818R 5′-ACGGAAACCTTGTTACGA-3′ [31]. Amplicons were purified using spin columns (DNA Clean & Concentrator-25; Zymo Research, Irvine, CA, USA) and sequenced (Sanger) using the same PCR amplification primers and an additional reverse primer 1174R, 5′-CCCGTGTTGAGTCAAA-3′ [32], when necessary to connect two ends. For phylogenetic analyses, similar sequences were retrieved from the PR2 database [33] based on BLAST similarity, and two environmental homologs (GenBank Acc. FJ537342 and FJ537336) were retrieved from NCBI GenBank for the undescribed haptophyte taxon, which was not affiliated with any reference sequence from the PR2 database. Sequence alignments including 39 isolates, 29 reference and 2 outgroup taxa were created with MAFFT v7.450 using the G-INS-i algorithm [34] in Geneious R11.1.5 (http://www.geneious.com) [35]. Terminal sites that lacked data for any of the sequences were trimmed and any sites with greater than 25% gaps were removed from the alignment, which generated a total sequence length of 1617 bases. Phylogenetic analysis was performed using MrBayes v3.2.6 in Geneious R11.1.5 [36] with two runs of four chains for 1,000,000 generations, subsampling every 200 generations with burn-in length 100,000, under the GTR substitution model. The Bayesian majority consensus tree was further edited within iTOL v5 [37]. All 18 S rRNA gene sequences were deposited in GenBank with accession numbers MZ611704–MZ611740; MN615710–MN615711.Microscopic observationThe average diameter of flagellates in the exponential growth phase (n = 20 cells per strain) was measured by transmitted light microscopy using image analysis software (NIS-Elements AR, Nikon, Minato City, Tokyo, Japan) calibrated with a stage micrometer. Equivalent spherical diameter (ESD) and biovolumes were calculated assuming spherical cells. Chloroplasts were visualized by autofluorescence under epifluorescence microscopy. An average ESD of 0.64 µm was used for Prochlorococcus prey [24].Visual evidence of phagocytosis was obtained by adding fluorescent beads (0.5 μm YG Fluoresbrite Microspheres; Polysciences) to each culture. Samples post incubation (~2 h) were fixed with an equal volume of 4% ice-cold glutaraldehyde, and subsamples (20 µL) were mounted on a glass slide under a coverslip. Paired images captured using epifluorescence and transmitted light microscopy (Olympus BX51 with Leica DFC 7000 T color digital camera) were overlain to identify cells with ingested beads.Grazing experimentsLong-term grazing experiments were conducted for all 39 grazers, and 31 were used to quantify grazing rates based on rates of disappearance of Prochlorococcus cells, which persist but do not readily grow in K-N medium [24]. Rates were not calculated for eight isolates (seven Florenciella and one DictyX) because they were sampled at a lower frequency. Fifteen isolates representing all genera (or approximately genus-level clades) were examined in more detail by replicating grazing experiments two times (marked in bold in Supplementary Table S1), while the remaining sixteen isolates were tested once to survey within- and across- genus variation. Prior to the experiments, all grazer cultures were maintained/acclimated in the experimental medium (K-N with prey). Experiments were initialized by inoculating late-exponential-phase grazers into fresh K-N medium at a final concentration of ~103 flagellates mL−1, and adding live, unstained prey at a final concentration of 2–3 × 106 Prochlorococcus mL−1. Grazers were incubated for 3–8 days in total, depending on how fast prey were ingested. To minimize carryover of dissolved nutrients and prey growth in the grazing experiments, Prochlorococcus were grown to stationary phase in Pro99 medium, then pelleted (2000 RCF for 3–5 min) and resuspended in fresh K-N medium prior to addition to control and experimental cultures. Control cultures of grazer without added prey and prey without grazers were included during each grazing experiment to confirm that grazer growth and prey removal were attributable to grazing. Cell concentrations of prey and grazers were measured every 12–24 h by flow cytometry of glutaraldehyde-fixed samples at final concentration of 0.5% (Attune NxT; Thermo Fisher Scientific, Waltham, MA, USA, and CytoFLEX; Beckman Coulter Life Sciences, Indianapolis, IN, USA). Populations were distinguished based on light scatter and pigment autofluorescence and occasionally confirmed with DNA fluorescence (stained post-sampling with DNA stain SYBR Green I). Ambient bacteria concentrations, monitored using DNA stains, were ≤10% of the added Prochlorococcus at the start and during most periods of all grazing experiments.Grazing rates and biovolume conversion efficiencyWe calculated ingestion rates, I (prey grazer−1 h−1) for each grazer as:$$I = frac{{P_t – P_{t + 1}}}{{G_{{{{{rm{avg}}}}}}(T_{t + 1} – T_t)}}$$
(1)
where (P_t) and (P_{t + 1}) are the prey abundance at sampling interval t and t + 1 (cells mL−1), Gavg is the arithmetic mean grazer abundance (cells mL−1) over the time interval, and ((T_{t + 1} – T_t)) is the time (h) between two sampling intervals. Clearance rates (C, nL grazer−1 h−1) were calculated by dividing the ingestion rate by average prey concentration over the same interval, and specific clearance rate (body volume grazer−1 h−1) was calculated by dividing the clearance rate by cellular biovolume (µm3) of each grazer. Equation (1) uses a linear approximation of prey and grazer trajectories over the sampling interval, which was appropriate for our data where change could be relatively linear, concave-up, or concave-down. Other commonly used ingestion rate calculations assume exponential prey decline [38] and/or exponential grazer growth [39].Clearance rates over time in each experiment were assessed visually to obtain a representative series of rates that minimized potential influence of modest prey growth/decline observed in control cultures (Supplementary Fig. S1), as well as potential slowing of ingestion as the grazer neared carrying capacity or depleted prey to a low concentration. For each experiment a contiguous set of relatively constant rates were used to calculate a mean clearance rate. This assessment sometimes excluded the first 12–24 h, but not when removal rates were particularly fast. Intervals when the grazer neared carrying capacity were also often excluded, if grazing rates slowed down. To assess whether clearance rates increased as prey were depleted we plotted clearance rate and ingestion rate as a function of prey concentration (Supplementary Figs. S2 and S3). In general, there was no relationship between clearance rate and prey concentration, and ingestion increased linearly with prey concentration. These patterns imply that the prey concentrations in this experiment did not saturate the ingestion rates of these grazers. Under non-saturating prey concentrations the average clearance rate over an experiment should be a good estimate of the maximum clearance rate (Cmax). Consistent with this interpretation, functional responses fit to these experiments yielded Cmax estimates that were similar to the reported average clearance rates. Because these experiments were not designed with a sufficient range of prey density to esimate functional responses, we do not report the Cmax estimates, but we note here that our reported average clearance rates may be useful as approximate Cmax numbers in future work.Six grazers representing three classes (three dictyochophytes, two haptophytes, and one chrysophyte) were further investigated to determine functional grazing responses using a wide range of initial prey densities (105–107 cells mL−1). Functional responses were modeled using the Holling type II curve, (I = frac{{I_{{{{{rm{max}}}}}}P}}{{P + frac{{I_{{{{{rm{max}}}}}}}}{{C_{{{{{rm{max}}}}}}}}}}), where I is the ingestion rate over a sampling interval (Eq. 1), Imax is the maximum ingestion rate, Cmax is the maximum clearance rate and P is the arithmetic mean prey density between two sampling points. This curve was fit to ingestion rate data using maximum likelihood with R package bbmle [40].For 31 isolates we calculated the amount of grazer biovolume created per prey biovolume consumed, using data from the same grazing experiments used to calculate grazing rates. It was calculated based on the following formula:$$E = frac{{(F_{{{{{rm{f}}}}}} – F_{{{{{rm{i}}}}}})}}{{(P_{{{{{rm{i}}}}}} – P_{{{{{rm{f}}}}}})}}frac{{B_{{{{{rm{F}}}}}}}}{{B_{{{{{rm{p}}}}}}}}$$
(2)
where Ff and Fi are the final and initial flagellate concentrations, Pi and Pf are initial and final prey concentrations in each culture, and BF and BP are the cellular biovolume of prey and grazer. We refer to the quantity E as the biovolume conversion efficiency, and we use it as an indicator of physiological differences among diverse mixotrophs. Note that biovolume conversion efficiency can be greater than 1, if prey have greater nutrient:biovolume than the grazer.Quantitative PCRReal-time, quantitative PCR (qPCR) was performed to quantify the 18S rRNA gene abundances of representative mixotroph groups discriminated at approximately the genus level, including Florenciella, Rhizochromulina and another undescribed clade within the class Dictyochophyceae; Chrysochromulina and another undescribed clade within the division Haptophyta; clade H in the class Chrysophyceae; and Triparma eleuthera and Triparma mediterranea in the class Bolidophyceae. Primers (Supplementary Table S2) were designed to target a short region (95–176 bases) of the 18S rRNA gene and meet basic criteria (≤2 °C difference in melting temperature between members of a pair, %G + C content between 45 and 65%, ≤1 degenerate position per primer, no predicted primer dimers). Sequences considered targets for a given primer set had ≤1 mismatch across both primers, which included all or most known members within the corresponding targeted clade. Members in the nearest non-targeted clade had ≥3 mismatches distributed across both primers. Efficiency and specificity of the synthesized primers (IDT Inc., Coralville, IA, USA) was tested by ensuring there was specific amplification (qPCR followed by melting curve analysis and gel electrophoresis) when using DNA from cultures within the targeted group and no amplification when using DNA from cultures close to, but outside of the targeted group (Supplementary Table S3). Empirical observations of amplification success using control cultures were used to infer whether species known only by environmental sequences were likely to amplify with a given primer set (Supplementary Fig. S4).In situ gene abundances were quantified in water samples collected from Station ALOHA at 5, 25, 45, 75, 100, 125, 150, and 175 m, during HOT cruise numbers 259 (Jan), 262 (Apr), 264 (Jul), and 266 (Oct) of 2014. Seawater (ca. 2 L) was filtered through 0.02 μm pore-size, aluminum oxide filters (Whatman Anotop, Sigma-Aldrich, Saint Louis, MO, USA) and stored at −80 °C. Genomic DNA of both grazer cultures and environmental samples was extracted (MasterPure Complete DNA and RNA Purification Kit; Epicentre) as described elsewhere [41]. Four replicated PCR reactions (10 μL) were carried out for each sample except for Triparma (duplicates) and consisted of 5 μl of 2× PowerTrack SYBR Green Master Mix (Thermo Fisher Scientific, USA), 10 ng environmental DNA, 500 nM of each primer, and nuclease-free water. Reactions were run on an Eppendorf Mastercycler epgradient S realplex2 real-time PCR instrument. Each run contained fresh serial dilutions (1–6 log gene copies) of target-specific, 750-bp synthetic standards (gBlocks, IDT) prepared in triplicate. The cycling program included an initial denaturation step of 95 °C for 2 min, followed by 40 cycles of 95 °C for 5 s and 55 °C for 30 s. Specificity of amplification was checked with a melting curve run immediately after the PCR program and occasionally, by gel electrophoresis. Amplification efficiencies ranged from 95% to 106% for all the primers.To convert gene copies to cell numbers, 18S rRNA gene copy number per cell−1 was determined for representative isolates in the seven targeted genera/clades. Known quantities of cultured cells (106–107 cells) from each isolate with 2–8 replicates were pelleted at 4000 RCF for 15 min at 4 °C. DNA was extracted from the pelleted cells (MasterPure Complete DNA and RNA Purification Kit, Lucigen), quantified by fluorometry (Qubit, Invitrogen) and the extract volume adjusted to achieve a DNA concentration of 10 ng µL−1. The expected number of eukaryotic cells µL−1 of extract was calculated as the difference between the total cells in the sample prior to centrifugation and in the supernatant afterward (as determined by flow cytometry) divided by the final extract volume. Copy number of the 18S rRNA genes µL−1 of extract was determined by qPCR with the appropriate group-specific primers. The resulting value of gene copies µL−1 was divided by the equivalent number of eukaryotic cells µL−1 in the extract (assuming 100% extraction efficiency) to derive minimum estimates of gene copies cell−1. An average value for representatives within each genus/clade (1–5 isolates) was used to calculate in situ cell concentrations for the genus. These derived in situ abundances were compared to flow cytometric counts of total photosynthetic picoeukaryotes at Station ALOHA obtained from the Hawai’i Ocean Time-series Data Organization and Graphical System (https://hahana.soest.hawaii.edu/hot/hot-dogs/).Global distribution revealed through Tara Oceans 18S rRNA metabarcodesTo estimate the relative abundance of the OTUs closely related to our diverse isolates on a broader geographic scale, we searched the 18S rRNA-V9 sequence data from the 0.8–5 µm fraction of surface water sampled at 40 stations by the Tara Oceans project (http://taraoceans.sb-roscoff.fr/EukDiv/). Reads for ‘Tara lineages’ with highest similarity (E-value < 10−15) to each of our targeted clades (Supplementary Table S1) were expressed as a fraction of total reads excluding dinoflagellates but included all other Tara Oceans phytoplankton ‘taxogroups’: Bacillariophyta, Bolidophyceae, Chlorarachnea, Chlorophyceae, Chrysophyceae/Synurophyceae, Cryptophyta, Dictyochophyceae, Euglenida, Glaucocystophyta, Haptophyta, Mamiellophyceae, Other Archaeplastida, Other Chlorophyta, Pelagophyceae, Phaeophyceae, Pinguiophyceae, Prasino-Clade-7, Pyramimonadales, Raphidophyceae, Rhodophyta and Trebouxiophyceae. Dinoflagellates were excluded because of the difficulty in assigning phototrophic vs. heterotrophic status to all taxa, and because nearly all dinoflagellate reads were from a single, poorly annotated OTU that was also highly abundant in larger size fractions. More