Philippa Kaur

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

Analytical targetsTwo species of undulation motion rays with different pectoral fin shapes bred in KAIYUKAN were analyzed: sharpnose stingray Dasyatis acutirostra and pitted stingray Dasyatis matsubarai (Fig. 1a,b). Blender 2.7925 was used to construct stingray models from pictures26,27 as accurately as possible; Blender is a free and open-source 3D creation suite used to make realistic characters for movies, etc. Detailed information on how to construct models using Blender is provided in our previous paper28. To focus on the effects of pectoral fin movements, we did not consider the body’s shape as in the previous studies12,29. The height and disk width (WD) of all models were set to 0.01 m and 0.44 m, respectively, considering the previous studies30,31. The disk length of each model was determined from WD, referring to the aspect ratio of the rays’ photographs26,27; the disk length (LD) of D. acutirostra and D. matsubarai are 0.348 m and 0.344 m, respectively.Figure 1Analytical targets and description of motion. (a) Analytical model of D. acutirostra. (b) Analytical model of D. matsubarai. (c) Description of motion, (d) the relationship between any two points on the surface before and after the deformation.Full size imageMotionThe motion was given to satisfy the following equations:$$z = left{ {begin{array}{*{20}l} {A;sin left( {omega left( {t – kTleft( {frac{{angl{text{e}}left( {x_{i} ,y_{i} } right) – 10^{{text{o}}} }}{{All;angl{text{e}}}} – theta } right)} right)h_{1} h_{2} } right.} hfill & {left( {10^{{text{o}}} le angl{text{e}}left( {x_{i} ,y_{i} } right) le 170^{{text{o}}} } right)} hfill \ {A;sin left( {omega left( {t – kTleft( {frac{{350^{{text{o}}} – left( {angl{text{e}}left( {x_{i} ,y_{i} } right) – 10^{{text{o}}} } right)}}{{All;angl{text{e}}}}} right)} right)h_{1} h_{2} } right.} hfill & {left( {190^{{text{o}}} le angl{text{e}}left( {x_{i} ,y_{i} } right) le 350^{{text{o}}} } right)} hfill \ end{array} } right.$$
(1)
$$begin{array}{c}{h}_{1}=a{r}_{i}^{3}+b{r}_{i}^{2}+c{r}_{i}end{array}$$
(2)
$$h_{2} = left{ {begin{array}{*{20}l} {dleft( {angleleft( {x_{i} ,y_{i} } right) – 10^{ circ } } right)^{2} + eleft( {angleleft( {x_{i} ,y_{i} } right) – 10^{ circ } } right)} hfill & {left( {10^{ circ } le angleleft( {x_{i} ,y_{i} } right) le 170^{ circ } } right)} hfill \ {dleft( {350^{ circ } – left( {angleleft( {x_{i} ,y_{i} } right) – 10^{ circ } } right)} right)^{2} + eleft( {350^{ circ } – left( {angleleft( {x_{i} ,y_{i} } right) – 10^{ circ } } right)} right)} hfill & {left( {190^{ circ } le angleleft( {x_{i} ,y_{i} } right) le 350^{ circ } } right)} hfill \ end{array} } right.$$
(3)
$$begin{array}{c}{left({r}_{i}-{r}_{i-1}right)}^{2}+{left({z}_{i}-{z}_{i-1}right)}^{2}={left({r}_{i}^{mathrm{^{prime}}}-{r}_{i-1}^{mathrm{^{prime}}}right)}^{2}+{left({z}_{i}^{mathrm{^{prime}}}-{z}_{i-1}^{mathrm{^{prime}}}right)}^{2}end{array}$$
(4)
$$begin{array}{c}angleleft({x}_{i},{y}_{i}right)= angleleft({x}_{i}^{^{prime}},{y}_{i}^{^{prime}}right).end{array}$$
(5)
Equation (1) represents the amount of movement of the model surface in the z-axis direction, where (A) is the amplitude of the pectoral fin tip, (omega) is the angular velocity, (t) is time, (k) is the wavenumber, (T) is the period, angle(({x}_{i},{y}_{i})) is the angle made by the line connecting the center of rotation and any point (({x}_{i},{y}_{i})) on the model surface with the x-axis, Allangle is the range where the motion is given (160°), and (theta) is the phase difference between the movements of the right and left pectoral fins (Fig. 1c). ({h}_{1}) is the weighting from the center to the radial direction: it is necessary to set the amplitude at the ray’s center to zero and increase the amplitude toward the pectoral fin tip (a = 119.786, b = -7.957, c = 0.498). ({h}_{2}) is the weighting in the circumferential direction: it is necessary to increase the amplitude from the anterior to the tip of the pectoral fin and decrease the amplitude from the tip of the pectoral fin to the posterior (d = − 1.563 × 10–4, e = 0.025). Equation (4) is the condition in which the distance between two neighboring points in the same radial direction is equal before and after the movement (Fig. 1d). (r) is the distance between the center of rotation and any point (({x}_{i},{y}_{i})), defined as (sqrt{{x}_{i}^{2}+{y}_{i}^{2}}). Equation (5) defines angle(({x}_{i},{y}_{i})) as being constant before and after the move (Fig. 1c). The variables after the move are marked with ‘. Variables used in the analysis are A = 0.089 m, T = 0.499, k = 1.270, and ω = 12.599 rad/s. Videos of the created motion from the front and the side are shown in the “Supplement” (Supplement Movies 3, 4).Analytical conditionsAnalysis cases were conducted with eight conditions: two types of pectoral fin shape (Fig. 1a,b) and four types of phase difference (0 (T), 0.25 (T), 0.5 (T), and 0.75 (T)). These conditions were set for investigating the effects of phase differences between left and right pectoral fin movements on swimming and how these effects vary with pectoral fin shape.Numerical methodsA CFD simulation of the ray models in the water flow was performed using OPENFOAM, an open-source finite volume method CFD toolbox32, to calculate the forces acting on the rays in each axial direction and the moment around each axis. The governing equations were the continuity equation and the three-dimensional incompressible Reynolds-averaged Navier–Stokes equation, expressed by:$$begin{array}{*{20}c} {nabla cdot u = 0} \ end{array}$$
(6)
$$begin{array}{*{20}c} {frac{partial u}{{partial {text{t}}}} + nabla cdot left( {uu} right) = – nabla p + nabla cdot left( {vnabla u} right) + nabla cdot left[ {nu left{ {left( {nabla u} right)^{T} – frac{1}{3}nabla cdot uI} right}} right], } \ end{array}$$
(7)
where (u) is the velocity vector, t is the time, p is the static pressure divided by the reference density, (nu) is the kinematic viscosity, and I is the unit tensor. The Reynolds number was defined regarding the previous studies3 as:$$begin{array}{c}{R}_{e}=frac{U{L}_{D}}{nu },end{array}$$
(8)
where (U) (/ms) is the given flow speed3 (1.5 × LD/ms), ({L}_{D}) (m) is the length of the ray models, and (nu) is the kinematic viscosity of water at 20 °C (1.0 × 10–6 m2/s). The Reynolds number in this study is 1.8 × 105; considering this, we used the k–ω shear stress turbulence model33,34. The k–ω shear stress turbulence model is a type of Reynolds-averaged Navier–Stokes equation (RANS) turbulence model that is widely used to calculate for the fish swimming flow35,36,37. The overset grid method was used in this study; it is a generic implementation of overset meshes. For both static and dynamic cases, cell-to-cell mapping between multiple, disconnected mesh regions is employed to generate a composite domain38,39. This method permits complex mesh motions and interactions without the penalties associated with deforming meshes. The process is described in detail by Noack40. The calculation volume was 5.4 WD in length, 5.4 WD in height, and 5.4 WD in width (Fig. 2a,b). A hexahedral volume mesh was created using the snappyHexMesh of OPENFOAM. The fluid region was divided into two parts: the overset region and the background region (Fig. 2a,b). The overset region moves and transforms to match the motion of the ray and was made with fine meshes around the analysis target and coarse meshes in the outlying areas; a one-layer boundary layer mesh was created around the analysis target. The overset region shape is an ellipsoid (Fig. 2a,b). The minimum mesh volume is 7.3 × 10–10 (m3), and the maximum mesh volume is 2.6 × 10–2 (m3). The total number of meshes was 9.0 × 105. At the outlet boundary, the average static relative pressure was set to 0 Pa. The surfaces of the fish model were formed into non-slip surfaces.Figure 2Meshes for CFD simulation and differences in force between different meshes. (a) Meshes at the coronal plane of the whole fluid region. (b) Frontal cross-section of the fluid region at the green line in (a). The red region is the overset region. (c,d) Comparison of the instantaneous drag coefficient and the moment coefficient around the y-axis of D. matsubarai between the coarse, fine, and dense mesh.Full size imageThe drag coefficient ({C}_{D}left(tright)), the lateral force coefficient ({C}_{l}left(tright)), the lift coefficient ({C}_{L}left(tright)), the moment coefficient around the x-axis ({C}_{mx}left(tright)), the moment coefficient around the y-axis ({C}_{my}left(tright)) and the moment coefficient around the z-axis ({C}_{mz}left(tright)) were calculated as:$$begin{array}{c}{C}_{D}left(tright)=frac{Dleft(tright)}{frac{1}{2}rho {U}^{2}{L}_{D}{W}_{D}}end{array}$$
(9)
$$begin{array}{c}{C}_{l}left(tright)=frac{lleft(tright)}{frac{1}{2}rho {U}^{2}{L}_{D}{W}_{D}}end{array}$$
(10)
$$begin{array}{c}{C}_{L}left(tright)=frac{Lleft(tright)}{frac{1}{2}rho {U}^{2}{L}_{D}{W}_{D}}end{array}$$
(11)
$$begin{array}{c}{C}_{mx}left(tright)=frac{{M}_{psi }left(tright)}{frac{1}{2}rho {U}^{2}{L}_{D}^{2}{W}_{D}}end{array}$$
(12)
$$begin{array}{c}{C}_{my}left(tright)=frac{{M}_{phi }left(tright)}{frac{1}{2}rho {U}^{2}{L}_{D}^{2}{W}_{D}}end{array}$$
(13)
$$begin{array}{c}{c}_{mz}left(tright)=frac{{M}_{theta }left(tright)}{frac{1}{2}rho {U}^{2}{L}_{D}^{2}{W}_{D}},end{array}$$
(14)
where (Dleft(tright)) is the calculated drag, (lleft(tright)) is the calculated lateral force, (Lleft(tright)) is the calculated lift, ({M}_{psi }left(tright)) is the calculated moment around the x-axis, ({M}_{phi }left(tright)) is the calculated moment around the y-axis, ({M}_{theta }left(tright)) is the calculated moment around the z-axis, and (rho) (kg/m3) is the density of water at 20 °C (998 kg/m3). As shown in a previous study41. the propulsive efficiency (eta) is defined as the ratio of output power ({P}_{o}) to input power ({P}_{e}) which can be written as:$$begin{array}{c}{P}_{o}left(tright)=frac{1}{T}{int }_{0}^{T}Dleft(tright)Udtend{array}$$
(15)
$$begin{array}{c}{P}_{e}left(tright)=frac{1}{T}{int }_{0}^{T}left[Dleft(tright)dot{x}left(tright)+lleft(tright)dot{y}left(tright)+Lleft(tright)dot{z}left(tright)right]dtend{array}$$
(16)
$$begin{array}{c}eta =frac{{P}_{o}}{{P}_{e}}.end{array}$$
(17)
An in-house program calculated the forces acting on rays in each axial direction and the moment around each axis. The numerical method’s validity and reliability were verified by comparing previous experimental and numerical analytical studies of heaving and pitching on airfoil naca001341. A high degree of similarity to previous studies was confirmed; the mean difference in the propulsive efficiency from the previous study of analysis was 5%, and the difference from the previous study of the experiment was 9%. Detailed information such as mesh, length, and velocity, of this analysis method’s verification is provided in the “Supplement”.A grid sensitivity study was conducted using three meshes: coarse, fine, and dense. The coarse mesh has 8.1 × 105 elements, the fine mesh has 9.0 × 105 elements, and the dense mesh has 9.9 × 105 elements. The analysis was conducted using a condition with no phase difference of D. matsubarai. As shown in Fig. 2c,d, the drag coefficient and the moment coefficient around the y-axis are almost the same when the mesh is fine and when the mesh is dense. The mean drag and propulsive efficiency error of fine and coarse meshes are 2.7% and 3.5%, respectively. The fine mesh was used in all simulation cases considering accuracy. We used the same meshes for all cases. 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

Clare, E., Barber, B., Sweeney, B., Hebert, P. & Fenton, M. Eating local: Influences of habitat on the diet of little brown bats (Myotis lucifugus). Mol. Ecol. 20, 1772–1780 (2011).CAS 
PubMed 
Article 

Google Scholar 
Larter, N. C. & Gates, C. C. Diet and habitat selection of wood bison in relation to seasonal changes in forage quantity and quality. Can. J. Zool. 69, 2677–2685 (1991).Article 

Google Scholar 
Veloso, C. & Bozinovic, F. Dietary and digestive constraints on basal energy metabolism in a small herbivorous rodent. Ecology 74, 2003–2010 (1993).Article 

Google Scholar 
Hawke, T., Bates, H., Hand, S., Archer, M. & Broome, L. Dietary analysis of an uncharacteristic population of the Mountain Pygmy-possum (Burramys parvus) in the Kosciuszko National Park, New South Wales, Australia. PeerJ 7, e6307 (2019).PubMed 
PubMed Central 
Article 

Google Scholar 
Pearce-Higgins, J. W. Using diet to assess the sensitivity of northern and upland birds to climate change. Clim. Res. 45, 119–130 (2010).Article 

Google Scholar 
Eitzinger, B. et al. Assessing changes in arthropod predator-prey interactions through DNA-based gut content analysis—variable environment, stable diet. Mol. Ecol. 28, 266–280 (2019).CAS 
PubMed 
Article 

Google Scholar 
Edgar, G. J. Predator-prey interactions in seagrass beds. II. Distribution and diet of the blue manna crab Portunus pelagicus Linnaeus at Cliff Head, Western Australia. J. Exp. Mar. Biol. Ecol. 139, 23–32 (1990).Article 

Google Scholar 
Beck, J. L., Peek, J. M. & Strand, E. K. Estimates of elk summer range nutritional carrying capacity constrained by probabilities of habitat selection. J. Wildl. Manag. 70, 283–294 (2006).Article 

Google Scholar 
DeYoung, R. W., Hellgren, E. C., Fulbright, T. E., Robbins, W. F. Jr. & Humphreys, I. D. Modeling nutritional carrying capacity for translocated desert bighorn sheep in western Texas. Restor. Ecol. 8, 57–65 (2000).Article 

Google Scholar 
Hua, L. et al. Captive breeding of pangolins: current status, problems and future prospects. Zookeys 507, 99–114 (2015).Article 

Google Scholar 
Nielsen, J. M., Clare, E. L., Hayden, B., Brett, M. T. & Kratina, P. Diet tracing in ecology: Method comparison and selection. Methods Ecol. Evol. 9, 278–291 (2017).Article 

Google Scholar 
Galimberti, A. et al. DNA barcoding as a new tool for food traceability. Food Res. Int. 50, 55–63 (2013).CAS 
Article 

Google Scholar 
Soininen, E. M. et al. Shedding new light on the diet of Norwegian lemmings: DNA metabarcoding of stomach content. Polar Biol. 36, 1069–1076 (2013).Article 

Google Scholar 
Rees, G. N., Shackleton, M. E., Watson, G. O., Dwyer, G. K. & Stoffels, R. J. Metabarcoding demonstrates dietary niche partitioning in two coexisting blackfish species. Mar. Freshw. Res. 71(4), 512–517 (2019).Article 

Google Scholar 
Taberlet, P., Coissac, E., Pompanon, F., Brochmann, C. & Willerslev, E. Towards next-generation biodiversity assessment using DNA metabarcoding. Mol. Ecol. 21, 2045–2050 (2012).CAS 
PubMed 
Article 

Google Scholar 
Aylagas, E., Borja, Á., Irigoien, X. & Rodriguez-Ezpeleta, N. Benchmarking DNA metabarcoding for biodiversity-based monitoring and assessment. Front. Mar. Sci. 3, 96 (2016).
Google Scholar 
De Barba, M. et al. DNA metabarcoding multiplexing and validation of data accuracy for diet assessment: Application to omnivorous diet. Mol. Ecol. Resour. 14, 306–323 (2014).PubMed 
Article 

Google Scholar 
Kartzinel, T. R. et al. DNA metabarcoding illuminates dietary niche partitioning by African large herbivores. Proc. Natl. Acad. Sci. 112, 8019–8024 (2015).CAS 
PubMed 
PubMed Central 
ADS 
Article 

Google Scholar 
Lopes, C. et al. DNA metabarcoding diet analysis for species with parapatric vs sympatric distribution: A case study on subterranean rodents. Heredity 114, 525–536 (2015).CAS 
PubMed 
PubMed Central 
Article 

Google Scholar 
Guillerault, N., Bouletreau, S., Iribar, A., Valentini, A. & Santoul, F. Application of DNA metabarcoding on faeces to identify European catfish Silurus glanis diet. J. Fish Biol. 90, 2214–2219 (2017).CAS 
PubMed 
Article 

Google Scholar 
Jakubavivciute, E., Bergström, U., Eklöf, J. S., Haenel, Q. & Bourlat, S. J. DNA metabarcoding reveals diverse diet of the three-spined stickleback in a coastal ecosystem. PLoS ONE 12, e0186929 (2017).Article 

Google Scholar 
Grant, T. & Fanning, D. Platypus 4th edn. (CSIRO Publishing, 2007).Book 

Google Scholar 
Hawke, T. et al. Long-term movements and activity patterns of platypus on regulated rivers. Sci. Rep. 11, 1–11 (2021).MathSciNet 
Article 

Google Scholar 
Gregory, J., Iggo, A., McIntyre, A. & Proske, U. Receptors in the bill of the platypus. J. Physiol. 400, 349 (1988).CAS 
PubMed 
PubMed Central 
Article 

Google Scholar 
McLachlan-Troup, T., Dickman, C. & Grant, T. Diet and dietary selectivity of the platypus in relation to season, sex and macroinvertebrate assemblages. J. Zool. 280, 237–246 (2010).Article 

Google Scholar 
Harrop, C. & Hume, I. Digestive tract and digestive function in monotremes and nonmacropod marsupials. Compar. Physiol. Primitive Mamm. 4, 63–77 (1980).
Google Scholar 
Klamt, M., Davis, J. A., Thompson, R. M., Marchant, R. & Grant, T. R. Trophic relationships of the platypus: Insights from stable isotope and cheek pouch dietary analyses. Mar. Freshw. Res. 67, 1196–1204 (2016).CAS 
Article 

Google Scholar 
Faragher, R., Grant, T. & Carrick, F. Food of the platypus (Ornithorhynchus anatinus) with notes on the food of brown trout (Salmo trutta) in the Shoalhaven River, NSW. Austral. J. Ecol. 4, 171–179 (1979).Article 

Google Scholar 
Grant, T. R. Food of the platypus, Ornithorhynchus anatinus (Ornithorhynchidae: Monotremata) from various water bodies in New South Wales. Aust. Mammal. 5, 135–136 (1982).
Google Scholar 
Marchant, R. & Grant, T. The productivity of the macroinvertebrate prey of the platypus in the upper Shoalhaven River, New South Wales. Mar. Freshw. Res. 66, 1128–1137 (2015).Article 

Google Scholar 
Krueger, B., Hunter, S. & Serena, M. Husbandry, diet and behaviour of platypus Ornithorhynchus anatinus at Healesville Sanctuary. Int. Zoo Yearbook 31, 64–71 (1992).Article 

Google Scholar 
Thomas, J. L., Handasyde, K. A., Temple-Smith, P. & Parrott, M. L. Seasonal changes in food selection and nutrition of captive platypuses (Ornithorhynchus anatinus). Aust. J. Zool. 65, 319–327 (2018).Article 

Google Scholar 
Hawke, T., Bino, G. & Kingsford, R. T. Damming insights: Variable impacts and implications of river regulation on platypus populations. Aquat. Conserv. Mar. Freshwat. Ecosyst. 31, 504–519 (2021).Article 

Google Scholar 
Bino, G., Kingsford, R. T., Grant, T., Taylor, M. D. & Vogelnest, L. Use of implanted acoustic tags to assess platypus movement behaviour across spatial and temporal scales. Sci. Rep. 8, 5117 (2018).PubMed 
PubMed Central 
ADS 
Article 

Google Scholar 
Chinnadurai, S. K., Strahl-Heldreth, D., Fiorello, C. V. & Harms, C. A. Best-Practice guidelines for field-based surgery and anesthesia of free-ranging wildlife. I. Anesthesia and Analgesia. J. Wildl. Dis. 52(2 Suppl), S14–27. https://doi.org/10.7589/52.2S.S14 (2016).PubMed 
Article 

Google Scholar 
Fiorello, C. V., Harms, C. A., Chinnadurai, S. K. & Strahl-Heldreth, D. Best-Practice guidelines for field-based surgery and anesthesia on free-ranging wildlife. Ii. Surgery. J. Wildl. Dis. 52(2 Suppl), S28–39. https://doi.org/10.7589/52.2S.S28 (2016).PubMed 
Article 

Google Scholar 
Vogelnest, L. & Woods, R. Medicine of Australian mammals: CSIRO Publishing (2008).Geller, J., Meyer, C., Parker, M. & Hawk, H. Redesign of PCR primers for mitochondrial cytochrome c oxidase subunit I for marine invertebrates and application in all-taxa biotic surveys. Mol. Ecol. Resour. 13, 851–861 (2013).CAS 
PubMed 
Article 

Google Scholar 
Leray, M. et al. A new versatile primer set targeting a short fragment of the mitochondrial COI region for metabarcoding metazoan diversity: Application for characterizing coral reef fish gut contents. Front. Zool. 10, 1–14 (2013).Article 

Google Scholar 
Greenfield, P. Greenfield hybrid analysis pipeline (GHAP). v1 (CSIRO, 2017).
Google Scholar 
Shackleton, M. et al. How does molecular taxonomy for deriving river health indices correlate with traditional morphological taxonomy?. Ecol. Indic. 125, 107537 (2021).Article 

Google Scholar 
Hebert, P. D., Cywinska, A., Ball, S. L. & Dewaard, J. R. Biological identifications through DNA barcodes. Proc. R. Soc. Lond. Ser. B Biol. Sci. 270, 313–321 (2003).CAS 
Article 

Google Scholar 
Ostell, J. & Sayers, E. W. Dennis A. Benson, Mark Cavanaugh, Karen Clark, Ilene Karsch-Mizrachi, David J. Lipman.Wilcoxon, F. Individual comparisons by ranking methods. Biometrics Bulletin, 1. In Breakthroughs in Statistics, 196–202 (Springer, 1992).Oksanen, J. et al. Package “vegan”. Community ecology package, version. Vol 2, No. 9, 1–295. (2013).R Development Core Team. R: A Language and Environment for Statistical Computing (R Foundation for Statistical Computing, 2021).
Google Scholar  More

A workflow for the calculation of CSI is presented in Fig. 1c. For all the analyses, we used the R v. 4.0.3 software environment20 implemented in RStudio v. 1.4.1103. The scripts used for each methodological step are available at the Figshare repository21. After data download from primary sources (PaleoClim and WorldClim), specifically for the CSI-future map set we performed an initial step aimed to obtain individual bioclimatic variables for each future time period for the four SSPs (Fig. 1b). To achieve this, the median values of nine GCMs were calculated in functions compiled in raster R package22 for each individual bioclimatic variable (see a few exceptions of number of GCMs used in Table 2).Table 2 General circulation models (GCM) used to construct the future map sets.Full size tableThe standard deviation (SD) was estimated as a measure of the amount of variation or dispersion along time series, from which the resulting output maps showed the places where climate conditions remained constant or variable across the temporal periods considered (Fig. 1a,b). The SD, as a way to identify stable/unstable climatic areas, was previously used in other climatic or evolutionary studies4,14. To compute the SD output rasters, we applied the mosaic function setting “fun = sd” from raster R package, calculating the SD for each pixel in the 12 time period rasters for CSI-past and five times for CSI-future, independently for each variable. The mosaic function was also used for the range calculation, with “fun = min” and “fun = max” to obtain the minimum and maximum values of input rasters, respectively, with a further step for subtracting maximum to minimum values.Specifically, for CSI-past, as it includes several time periods with sea-level dropping below the present level (T1, T3, T5, T6, T7, T8, T9; Fig. 1a), we applied a mask of the current land surface, i.e. taking the T12 (Anthropocene) as a template. With this additional step, we were able to remove those pixels (grid cells) currently under the sea but that were once emerged. Most of these pixels, however, were only emerged during the LGM (ca. 21 ka), thus having values for bioclimatic variables for just a single time period (instead of the 12 routinely used for the variability estimation). The inclusion of these areas would result in highly climatically stable regions (low SD values; Supplementary Fig. 1), but this would be an obviously biased result. In contrast, we did not remove those areas affected by the sea-level rising periods, as only three periods contained “NoData” values (T2, T4, T10; Fig. 1a). However, to take this fact into consideration, we created a raster file in which these areas submerged during warm periods are indicated (see Supplementary Fig. 1). Finally, for both CSI-past and CSI-future, the resulting SD values were normalized to values between 0 and 1, with 0 representing completely stable areas and 1 the most unstable ones.The next step was focused on the selection of a relatively uncorrelated set of variables for each map set. We used the removeCollinearity function from virtualspecies R package23 that estimates the correlation value among pairs of variables from a given number of random sample points (10,000 in present case) according to a given method (Pearson for the present case) and a threshold of statistic selected (r  > 0.8 as a cut-off value). The function removeCollinearity returns a list of uncorrelated variables according to the settings specified, randomly selecting just one variable from groups of correlated ones (see Table 1 for a complete list of variables used for each map set). As we compiled estimates of variability independently for each variable and map set (e.g. SD bio1 past, SD bio2 past, etc.), each user can define his own CSI, selecting the more interesting variables according to the case of study.The final CSI maps were obtained by summing the SD values of the variables selected and the subsequent outputs normalized (0 to 1) (Figs. 2–4). Histogram plots were represented with ggplot2 R package24 and maps were exported with ArcGIS v.10.2.2 (Esri, Redlands, California, USA 2014). The histograms were computed for these final CSI maps, which represent the frequency and distribution of CSI values. We presented the final CSI maps with two different colour ramp schemes with ArcGIS. The first consisted of defining equal interval breaks from 0 to 1. The second was based on defining 32 categories with different value breaks for past and future map sets according to the value frequency shown by the histogram plot, i.e. the category with the highest CSI values (no. 32) was 0.71–1 in the past map set and 0.356–1 in the future map set.Fig. 2Maps of Climate Stability Index (CSI) values for the past map set from Pliocene (3.3 Ma) to present (1979–2013), at 2.5 arc-min grid resolution. Colours range from blue for low standard deviation (SD) values, which represents areas with low climatic fluctuations (i.e, low values of CSI) during the period Pliocene–present, to red for high SD values, which shows areas where high climatic fluctuations would have taken place (i.e., high values of CSI). On the upper map, the colour ramp shows equal interval breaks. The histogram with frequency and distribution of CSI values is also shown. On the lower map, the colour ramp has been manually adjusted to a defined set of break values (see details in the text).Full size imageFig. 3Maps of Climate Stability Index (CSI) values for the future conditions (Shared Socioeconomic Pathways: SSP1-2.6, SSP2-4.5, SSP3-7.0, and SSP5-8.5) from present (1970–2000) to future (2100), at 2.5 arc-min grid resolution. Colours range from blue for low standard deviation (SD) values, which represents areas with low climatic fluctuations (i.e, low values of CSI) from present to future, to red for high SD values, which shows areas where high climatic fluctuations would have taken place (i.e., high values of CSI). The colour ramp shows equal interval breaks. The histogram with frequency and distribution of CSI values is also shown for each future scenario.Full size imageFig. 4Maps of Climate Stability Index (CSI) values for the future conditions (Shared Socioeconomic Pathways: SSP1-2.6, SSP2-4.5, SSP3-7.0, and SSP5-8.5) from present (1970–2000) to future (2100), at 2.5 arc-min grid resolution. Colours range from blue for low standard deviation (SD) values, which represents areas with low climatic fluctuations (i.e, low values of CSI) from present to future, to red for high SD values, which shows areas where high climatic fluctuations would have taken place (i.e., high values of CSI). The colour ramp has been manually adjusted to a defined set of break values (see details in the text).Full size image More

Study siteAn open field small-scale release of a GM strain of Anopheles mosquitoes was carried-out in July 2019 in the village of Bana in Western Burkina Faso (see Supplementary Fig. 5). The study was granted regulatory authorisation from the National Biosafety Agency (NBA) (order No. 2018-453/MESRSI/SG/ANB of 10 August 2018 authorising the controlled release of genetically modified sterile male mosquitoes) and institutional ethical permission from the Institutional Ethics Committee for Research in Health Sciences: CEIRES (No. A-003/2019-CEIRES granted on January 9th 2019) and a programme of engagement established community acceptance. Details of the extensive stakeholder and communication processes and activities that were conducted in preparation of this release will be published elsewhere. The village of Bana is located in Western Burkina Faso (12°36′00″N, 3°28′59″W), 23 km west of the city of Bobo-Dioulasso.Bana has two main inhabited agglomerations of similar size: Bana Centre (administrative area) and Bana Market (economic area), separated by a 1.5 km unpopulated land band, crossed by a small semi-permanent river and a forest (see Supplementary Fig. 5). In its entirety, the village comprises about 130 compounds for about 759 inhabitants (local census, IRSS 2014). This region is characterised by two seasons: a wet season from June to September and a dry season from November to April. The mean annual rainfall in the village is about 800 mm and the mean temperature is about 27 °C (22–32 °C)52.Study designThe study design followed the format of an MRR experiment with an intensive period of recaptures followed by several months of monitoring to confirm the disappearance of the transgene. Both the period (July) and design (MRR-like experiment) were informed by previous baseline entomological studies and MRR experiments conducted in the same village41,52. Given the low population size expected in July and to avoid over-sampling, a lower recapture effort (reduction of daily swarm sampling number) was implemented than in previous MRR studies performed in the same area.41 The month of July corresponds to the start of the rainy season, when regular rains and cooler weather promote mosquito survival, and the target population of A. coluzzii is at a much lower level than later in the rainy season41,52. In July, plant coverage is still sparse and males tend to seek refuge inside houses and can be captured in good numbers via indoor sampling52.GM sterile strain maintenanceThe mosquito strain used in the experiment was the genetically modified mosquito Anopheles coluzzii sterile male strain referred to as Ac(DSM)2 (for Anopheles coluzzii Dominant Sterile male strain 2). This strain is the product of local introgression (series of backcrosses) of the original Ag(DSM)2 (dominant sterile male on Anopheles gambiae G3 mosquitoes strain 2) with a local A. coluzzii wild-type (WT) colony (female DSM-carrier crossed with male WT)34. The importation of Ag(DSM)2 in Burkina Faso, introgression with local wild type background and maintenance were conducted under regulatory authorisation from the National Biosafety Agency (N°000002/MRSI/SG/ANB of October 21th 2016). The wild-type A. coluzzii strain used for introgression and maintenance of Ac(DSM)2 was colonised in July 2014 from gravid female adults collected in village 7 of the Kou valley (VK7) in western Burkina Faso. Both colonies were maintained in a dedicated ACL2 (Arthropod Containment Level 2) insectary located within the IRSS main campus at Bobo-Dioulasso, Burkina Faso.For general stock-keeping purposes, Ac(DSM)2 was reared in a dedicated and highly secured climate-controlled room at a temperature fixed at 27.4 °C (±0.2, 95% Confidence intervals) and a relative humidity of 76.3% (±3.2, 95% CIs). Rearing rooms have natural light via windows and were supplemented with an artificial lighting regime of LD 12/12 h photoperiod, including dusk (1 h) and dawn (1 h). Larvae were reared in plastic trays (20 × 30 cm) with 1 l of deionized water and fed with an optimised larvae diet regime53. When mosquito larvae reached their level 3 instar (L3) larvae stage they were sorted manually between transgenic and non-transgenic mosquito larvae using a fluorescent stereomicroscope (Olympus SZX7, 2-8 Honduras street, London, United Kingdom) and put in separated trays to continue their development till pupation. At the pupal stage a second round of sorting occurred to separate male and female (sexing) from both strains. The sexing was done under a basic stereomicroscope (Olympus SZX7 basic, 2–8 Honduras street, London, United Kingdom) using a thin soft brush. Pupae from each strain and sex were placed in small plastic cups inside separate fresh adult cages to emerge. Adults were kept in 30 × 30 × 30 cm insect cages (produced locally) and continuously supplied with 10% (w/v) glucose solution (made with deionized water). Each generation, adult female transgenic mosquitoes were mated with male mosquitoes from the wild-type colony and blood-fed with fresh rabbit’s blood, using a membrane feeder (Hemotek® feeder, Hemotek Ltd, Blackburn United Kingdom). Gravid females were allowed to oviposit in plastic Petri dishes containing a wet sponge covered with filter paper. Eggs were collected and hatched in plastic trays. First instar larvae (L1) were then redistributed into several trays to keep similar larvae abundance (about 250 L1 larvae per tray).In accordance with Mendelian inheritance, stock-maintenance crosses between Ac(DSM)2 females and wild type colony males are expected to generate ~50% hemizygous transgenic male and female progeny referred to as Ac(DSM)2 and 50% non-transgenic sibling with a wild-type phenotype referred to as WT-Ac(DSM)2. That the actual phenotypic proportions matched the expected ratio was checked at each generation a part of standard procedures of colony maintenance.Production, sexing, marking and transport of release mosquitoesReleased males were derived from the 41st backcross generation from strain importation. Assuming Mendelian inheritance, the proportion of residual non-local genetic background after so many generations would be negligible (= 0.541).In rearing the release mosquito cohort, some changes were made in the stock-keeping procedure to maximise the fitness of male mosquitoes to be released. These changes aimed to minimise male mosquito handling during the entire process (rearing, sorting, marking and transport). Crucially, no transgenic versus non-transgenic sorting was done at larval stage resulting in a mix of transgenic and non-transgenic sibling males in the release generation. Additionally, to minimise the number of transgenic female mosquitoes released during the study, male versus female sexing was done at both pupae (initial) and adult (complementary) stages leading to a very high sorting accuracy (over 99.5%). Pupae sexing followed the procedure described for stock maintenance. Next, adult sexing focused on removing the few females resulting from errors in pupal sexing. It consisted of removing those rare females from male mosquito cages through inspection by eye of cages and in using a heat source to attract females. Once spotted, these were removed from male release cages using a mouth aspirator.After pupae sexing, male pupae were placed in 25 × 25 × 25 cm emergence cages (made locally and specially designed to fit dimensions of the secured coolboxes used for secure transportation) at a density of ~1400 pupae per cage. Following adult emergence, and over the following days, the cages were inspected by eye daily to check for and remove any females that had not been detected during the pupae sexing process. This procedure led to a total of 15,384 male mosquitoes aged 3–7 days have emerged in 15 cages and ready for marking and release purposes. Screening of ~50 males randomly picked from each emergence cage was conducted in the ACL2 insectary and revealed a slight bias in favour of WT-Ac(DSM)2 sibling males while Ac(DSM)2 male represented 43.3% (39.7–46.9, 95% CIs) of all emerged males. Based on this genotypic ratio, it was estimated that the male release cohort was equivalent to about 6659 transgenic male mosquitoes Ac(DSM)2 and 8725 non-transgenic sibling mosquitoes called WT-Ac(DSM)2 sibling. All males were kept untouched and in the same cages throughout the whole process until being released.The marking process was performed inside the ACL2 insectary facility, and was carried out the day before field release to allow enough time for mosquito recovery, rest and feeding. The environmental conditions were similar to those used during mosquito production. The mosquitoes were marked directly in their cages by using a cloud dye dusting technique. Aside from being fast, this highly efficient marking procedure (100% of mosquitoes successfully marked) was developed to allow the dust-marking of males in their original emergence cages, thereby avoiding male handling and damage during the marking process. This marking technique consisted of injecting pressurised red fluorescent colour powder (Bioquip® Gladwick Rancho Dominguez, CA 90220, USA; Ref: 1162R) into the cages by using a 5 ml syringe and needle to create a cloud of powder. The cages were wrapped with aluminium foil on all sides to prevent the dust from escaping through the meshed walls. Forceful injection of small amounts of powder from different sides of the cages through the aluminium cover and side netting created a dense cloud of fluorescent powder inside the cages to mark all the mosquitoes. Following marking, sugar-water was available ad-libitum to all marked mosquitoes until field release.About 2 h before the release time, the marked mosquitoes within the mosquito cages were transferred from the IRSS insectary to the release site in Bana village. Before leaving the IRSS insectary, the mosquito cages were covered by a second layer of mosquito net for security purposes. The cages were then wrapped with damp towels and placed in lockable cool boxes dedicated to their transport into the field. After having been secured, the cool boxes containing marked mosquitoes were transported to the release site. The entire process complied carefully with all regulatory requirements related to the permissions received for maintenance, handling and the release of these genetically modified organisms in Burkina Faso.Release phaseAll marked mosquitoes were released on the same day at around 5 pm (about one hour before swarming) in the centre of Bana village by opening the travel cages and allowing free exodus. Mosquitoes that did not leave were counted and subtracted from the released total (n = 534, 3.5%). Taking into account mortality and based on the ratio of Ac(DSM)2 and their siblings previously established, a total of 14,850 male mosquitoes were effectively released, with estimated numbers of 6428 hemizygous transgenic male A. coluzzii mosquitoes Ac(DSM)2 and 8422 non-transgenic WT-Ac(DSM)2 siblings.Recapture phaseMosquito recapture activities started the same day of release (about 2 h after mosquito release) and took place daily for a period of 20 days after release. Two different recapture methods were used: swarm collections using sweep nets (SWN) and pesticides spray catches (PSC) inside houses.Swarm sampling started on the evening of the release day using a well-established sweep net collection method47,54. Previous surveys in the same village41 had allowed mapping of swarm location or natural markers where swarming repeatedly occurs. To ensure sampling across the whole study area, a stratified randomised sampling procedure was used to select and sample 15 mosquito swarms daily at dusk using the sweep net collection method. The area of Bana village and Bana Marché were divided in six and four zones, respectively. Zone 1 and 2 in Bana Village are areas of high swarm abundance and the design ensured that these were not over-represented in swarm collections. Each evening, the teams of capturers set-out to collect up to five swarms per zones depending on swarm availability (swarms are fewer and smaller in early July than later in the month). All mosquitoes captured in the swarms were transported in their sweep nets to the field laboratory and frozen until the next morning for processing. At this stage, a random sample of 15 swarms each day was picked for dust screening and genetic analyses.Pyrethroid spray catches started the morning following the release and continued for 19 days. A set of 20 compounds were sampled each day. The sampling design followed that established in baseline studies leading to the release and in previous MRR studies41. Ten of the compounds were selected completely randomly and the other ten are a fixed set of compounds distributed regularly across the whole village. For each compound selected, a single room (1 sleeping room) within one of the house of compound was chosen for sampling. Although some compounds were selected more than once during the recapture period days, a different room (from a different house inside the same compound when applicable) was selected and no room was sampled twice during the survey period.Pyrethroid spray catches started the morning following the release and continued for 19 days. A set of 20 compounds were selected randomly each day. For each compound selected, a single room (sleeping room) was chosen for sampling. Although some compounds were selected more than once during the seven days, a different room (from a different house inside the same compound when applicable) was selected and no room was sampled twice during the survey period.Captured mosquitoes were identified morphologically in the field using adult anopheline morphological identification keys developed by Holstein55 and a field stereomicroscope (Perfex Sciences® Zoom Pro, Reference: S0852Z5 Toulouse, France). All An. gambiae s.l. mosquitoes were counted, checked for fluorescent dust marking using a Biofinder portable ultraviolet illuminator (Vansky, Shenzhen, China) and preserved in 80% ethanol. The identification of each marked mosquito was confirmed independently by two well-trained members of the staff before conservation in individual 1.5 ml storage microtubes for further analysis. The non-dusted wild Anopheles mosquitoes were pooled (10 individuals per tube) and stored in similar conditions. The location of each collection was recorded and mapped using a GPS (Garmin GPS) device, series GPSMAP®62.2.3. For all recaptured mosquitoes, we calculated the straight line distance from the release point to the recapture location using a Euclidean dispersal distance56. In the present case, the space was assimilated to a two-dimensional orthogonal axis system where xl and yl represent the coordinates of the release point and xr and yr represent the coordinates of the recapture point56. Calculation of the estimated flight distance of the mosquitoes then used the following formula:$${EFD}=sqrt{{left({x}_{r}-{x}_{l}right)}^{2}+{left({y}_{r}-{y}_{l}right)}^{2}}$$
(1)
Ac(DSM)2 male identificationMolecular analysis of recaptured marked mosquitoes was performed by PCR, to identify the Ac(DSM)2 strain and distinguish them from their non-transgenic WT-Ac(DSM)2 siblings. This PCR analysis consisted of detecting the integration of the eGFP::I-PpoI of the DSM transgene which characterised the transgenic mosquito strain Ac(DSM)2. In addition, a molecular species-diagnostic was performed concomitantly using the PCR technique based on the detection of SINE 200× locus57 and this PCR served as a control for DNA integrity. Each mosquito was split into two parts (abdomen and thorax) using forceps. The abdomen was used for the PCR and processed for DNA extraction using ‘squish’ buffer (PCR reaction buffer). The thorax was stored in 80% ethanol at −20 °C. For each mosquito analysed, the same DNA extract was used for both eGFP::I-PpoI transgene detection (identification of Ac(DSM)2 transgenic mosquito) and SINE 200X locus detection (for specie identification and DNA quality control). The Ac(DSM)2 construct was detected using the primers: pBacR-fwd [ATCGGTCTGTATATCGAGGTTTATT] and pBacR-Rev [CTCTAATATTTTGCCAAATGAAGTGCC] targeting the piggyBacR region required for insertion of the transgene. PCR reactions used the Gotaq® PCR kit (GoTaq® G2 Flexi DNA Polymerase, reference: M829B, Promega Corporation, 2800 Woods Hollow Road·Madison, WI 53711-5399, USA).Monitoring of Ac(DSM)2 non-persistenceMonthly mosquito collections were carried out using PSC and swarm sampling to confirm the disappearance of the Ac(DSM)2 transgene from the release site. Monitoring collections were conducted monthly for seven months. This period of monitoring was justified by the regulatory requirement of describing the Ac(DSM)2 disappearance through failure to detect the Ac(DSM)2 transgene for a minimum period of three consecutive months and with high statistical power. During each month of survey, a randomised selection of 20 houses (one room per house) and 20 swarms was sampled. All collected mosquitoes were identified morphologically using identification keys and a field stereomicroscope. Mosquitoes from A. gambiae complex were counted and preserved in 80% (v/v) ethanol for subsequent molecular identification. Each month, a representative sample of collected mosquitoes (up to 300 when available, from both PSC and swarm sampling) was analysed using the Ac(DSM)2-specific and species-specific PCR diagnostics described above to detect whether any A. gambiae s.l. mosquitoes were carrying the DSM transgene.Bayesian inference of mosquito survival, movement and population sizeWe fitted the recapture data to a diffusion model to further investigate dispersal and survival of the marked Ac(DSM)2 and their sibling males, and also to estimate the number of mosquitoes in the background population. This model assumes that the released mosquitoes tend to move in a random manner, meaning they repeatedly take short randomly directed flights that are independent of one another and of the environment. As described below, however, our estimation procedure does also allow for small additional movements where mosquitoes are attracted into nearby swarms at swarming time (dusk), or nearby houses for resting behaviour.We write the diffusion equation as$${partial }_{t}u=D{partial }_{x}^{2}u,$$
(2)
where (u(x,t)) is the probability density of the location of a single marked mosquito at location (x) and time (t), conditional on the individual being alive, and (D) is the diffusion coefficient. Assuming a point release at time (t=0), the above equation has solution$$uleft(r,tright)=frac{{e}^{-frac{{r}^{2}}{4Dt}}}{4,pi D,t}$$
(3)
where (r) is the distance from the release point. We next assume that the released mosquitoes have a constant survival probability of (s) per day, so that the expected number of extant released mosquitoes on day (d) is (R{s}^{d}) where (R) is the number that were released. The expected number of released mosquitoes in a small area ({dA}) is then given by$$qleft(r,dright)=R{s}^{d}frac{{e}^{-frac{{r}^{2}}{4Dd}}}{4,pi D,d}{dA}.$$
(4)
We take three further steps to convert this equation for (q(r,t)) into a likelihood function for the spatio-temporal distribution of recaptures of either Ac(DSM)2 or their sibling males. First, we pool the recaptures on a given day, and made by a given method (either swarm sampling or PSC), by partitioning the study area into annuli centred on the release location. These annuli are the recapture regions, and the expected number of extant marked mosquitoes in a given annulus is the integral of (q(r,d)) over that annulus. This step, therefore, averages out the expected number of marked mosquitoes from the inner to the outer radius of each annulus, and the annulus widths set the scale at which small movements towards swarms or houses, where mosquitoes may be recaptured, are assumed to occur in addition to random movements that underpin the diffusion model. We set the width of each annulus to 50 m, based on our judgement that this distance balances the capacity to separate recaptures at different distances (this capacity reduces with width), with the confidence that movements towards swarms or houses will largely remain within annuli (this confidence increases with width).Second, we assume the observation probability of mosquitoes in a given sample (representing an annulus, capture method, and day), is the number of unmarked mosquitoes in the sample divided by the (unknown) unmarked population size in that annulus. The unmarked population is assumed to have a uniform density, that we will infer alongside the mobility and survival parameters. Finally, we assume the number of marked mosquitoes in a given sample is Poisson-distributed around the expected number.For the data from each recapture method, we used the likelihood function to sample a posterior distribution for the diffusion coefficients and survival rates of the two types of released male mosquitoes, and the density of the unmarked population. We assumed uniform priors with respect to all five parameters and used a Markov chain Monte Carlo algorithm based on Metropolis-Hastings sampling to sample the posterior distribution directly from the log-likelihood. For each analysis (swarm or PSC), we sampled for 100,000 iterations, of which we discarded the initial 20,000 as a transient and thinned the remainder by 100, giving 800 samples in total.Statistical analysisData were analysed using the software JMP 14 (SAS Institute, Inc.). All data were checked for deviations from normality and heterogeneity, and analyses were conducted using parametric and non-parametric methods as appropriate. General linear modelling with Poisson distribution was used to describe male recaptures as a function of genotype and time. Kruskall-Wallis and Mann-Whitney test was used to describe respectively male participation in swarm and Euclidian dispersal distances. General linear modelling with Poisson distribution was used to describe male recaptures as a function of genotype and time. Estimates of population size, survival, and mobility were calculated using a Bayesian approach as described above.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

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