Designed nanomolar small-molecule inhibitors of Ena/VASP EVH1 interaction impair invasion and extravasation of breast cancer cells

Battling metastasis through inhibition of cell motility is considered a promising approach to support cancer therapies. In this context, Ena/VASP-depending signaling pathways, in particular interactions with their EVH1 domains, are promising targets for pharmaceutical intervention. However, protein–protein interactions involving proline-rich segments are notoriously difficult to address by small molecules. Hence, structure-based design efforts in combination with the chemical synthesis of additional molecular entities are required. Building on a previously developed nonpeptidic micromolar inhibitor, we determined 22 crystal structures of ENAH EVH1 in complex with inhibitors and rationally extended our library of conformationally defined proline-derived modules (ProMs) to succeed in developing a nanomolar inhibitor (Kd=120 nM,MW=734 Da). In contrast to the previous inhibitor, the optimized compounds reduced extravasation of invasive breast cancer cells in a zebrafish model. This study represents an example of successful, structure-guided development of low molecular weight inhibitors specifically and selectively addressing a proline-rich sequence-recognizing domain that is characterized by a shallow epitope lacking defined binding pockets. The evolved high-affinity inhibitor may now serve as a tool in validating the basic therapeutic concept, i.e., the suppression of cancer metastasis by inhibiting a crucial protein–protein interaction involved in actin filament processing and cell migration.

Metastasis is a complex multistep process (1, 2) employing, among others, mechanisms governing actin cytoskeleton dynamics involving integrin signaling and actin regulatory proteins (3–5). So far, all approved antimetastatic drugs antagonize integrins (6) or inhibit downstream kinases (7, 8) (SI Appendix, Fig. S1). In the metastatic setting however, these drugs appear to have only limited success (9–13) and 5-y survival is not increasing satisfactorily (14, 15), making new approaches in antimetastatic drug development essential to meet this urgent medical need.The enabled/vasodilator stimulated phosphoprotein protein family (Ena/VASP) acts as a crucial hub in cell migration by linking actin filaments to invadopodia and focal adhesions (16–22). Due to their role in the transformation of benign lesions into invasive and metastatic cancer, Ena/VASP proteins are discussed as part of the invasive signature and as a marker of breast carcinogenesis (23–25). At the advanced tumor stage, the protein family is overexpressed (26–28), which has been shown to increase migration speed in vivo (29) and to potentiate invasiveness (30). Yet, no sufficiently potent probes to interfere with Ena/VASP in vivo have been reported.The three vertebrate Ena/VASP family members, enabled homolog (ENAH), VASP, and Ena-VASP-like (EVL), share a tripartite structural organization in which two Ena/VASP homology domains (EVH1 and EVH2) are separated by a more divergent proline-rich central part. Interactions of the EVH2 domain are involved in the elongation and protection of barbed-end actin filaments from capping proteins and tetramerization (31, 32). EVH1 folds into a structured globular domain that interacts with proteins at focal adhesions (33), the leading edge (34, 35), and invadopodia (36, 37) by recognizing the motif [F/W/L/Y]PxϕP (35, 38) (ϕ hydrophobic, x any; SI Appendix, Fig. S3) in poly-L-proline type II helix (PPII) conformation.In the course of our research into small molecules as potential inhibitors of protein–protein interactions (39) we recently in silico designed and stereo-selectively synthesized scaffolds, coined ProMs, which mimic pairs of prolines in PPII conformation (40). The modular combination of different ProMs thereby allowed us to generate nonpeptidic secondary-structure mimetics that fulfill the steric requirements of the addressed proline-rich motif-recognizing domain (41–47). For the EVH1 domain, our proof-of-concept study yielded a canonically binding, nontoxic, cell-membrane-permeable, 706-Da inhibitor 1 (Fig. 1A) composed of two different ProM scaffolds and 2-chloro-(L)-phenylalanine (2-Cl-Phe) (40). While the synthetic inhibitor 1 represents the compound with the highest reported affinity toward Ena/VASP EVH1 domains, a further improvement was required for in vivo experiments. Here we report successful structure-based optimization of inhibitor 1 based on 22 high-resolution crystal structures of ENAH EVH1 in complex with different inhibitors (SI Appendix, Tables S1–S6), including the well-resolved C-terminal binding epitope TEDEL of ActA from Listeria monocytogenes (48). Newly identified interaction sites adjacent to the C terminus of 1 were addressed by in silico designed and stereo-selectively synthesized modifications of the ProM-1 scaffold (Fig. 1A). While drastically increasing the affinity against a rather flat protein surface we conserved structural simplicity, low molecular weight, nontoxicity, and cell-membrane permeability. Potent compounds against Ena/VASP were shown to also act in vivo, i.e., by inhibiting cancer cell extravasation in zebrafish at only 1 μM, thereby paving the way for future preclinical studies.Open in a separate windowFig. 1.(A) Structure of the first-generation Ena/VASP EVH1 inhibitor 1. All compositions share the N-acetylated 2-chloro-phenylalanine unit (blue) attached to a central ProM-2 scaffold (red). Esterification of the C terminus renders the inhibitors cell-membrane permeable (40). (B) General (modular) architecture of nonpeptidic, conformationally preorganized inhibitors used in this study. Structural variation (pink) was achieved by replacing the C-terminal ProM-1 unit (green) by ProM-9, ProM-13, ProM-12, ProM-15, or ProM-17 (Table 1).     

Superelastic oxide micropillars enabled by surface tension–modulated 90° domain switching with excellent fatigue resistance

Superelastic materials capable of recovering large nonlinear strains are ideal for a variety of applications in morphing structures, reconfigurable systems, and robots. However, making oxide materials superelastic has been a long-standing challenge due to their intrinsic brittleness. Here, we fabricate ferroelectric BaTiO₃ (BTO) micropillars that not only are superelastic but also possess excellent fatigue resistance, lasting over 1 million cycles without accumulating residual strains or noticeable variation in stress–strain curves. Phase field simulations reveal that the large recoverable strains of BTO micropillars arise from surface tension–modulated 90° domain switching and thus are size dependent, while the small energy barrier and ultralow energy dissipation are responsible for their unprecedented cyclic stability among superelastic materials. This work demonstrates a general strategy to realize superelastic and fatigue-resistant domain switching in ferroelectric oxides for many potential applications.

Superelastic materials are capable of recovering large amount of nonlinear “plastic” strains, way beyond their linear elastic regimes (1–4). They are ideal for a variety of applications from morphing structures, reconfigurable systems, to robots (5–8). The effects have traditionally been associated with macroscopically compliant/ductile rubbers (2) or microscopically phase-transforming shape memory alloys (SMAs) (7–11). The only macroscopically brittle oxide recently discovered to be superelastic is ZrO₂-based micropillars or particles (12–20), which is realized via austenite-martensite phase transformation similar to SMAs. Although ultimate strengths approaching the theoretical limit have been demonstrated in nanoscale samples (21, 22), long fatigue life is elusive, which is arguably more important for most applications. As a matter of fact, poor fatigue life has been a long-standing challenge for oxide ceramics in general (23, 24). Even for ductile SMAs that enjoy excellent fatigue life, irrecoverable residual strains gradually accumulate over cycling, leading to substantial variations in stress–strain curves at different cycles (9, 10, 25). We overcome these difficulties by reporting superelastic barium titanate (BaTiO₃ [BTO]) micropillars enabled by surface tension–modulated 90° domain switching, which exhibit excellent fatigue resistance, while bulk BTO crystals or ceramics are rather brittle. The demonstration of over one million cycles of loading and unloading without accumulating residual strains or noticeable variation in stress–strain curves is unprecedented among superelastic materials.BTO is a ferroelectric oxide exhibiting modest piezoelectric strains around 0.1 to 0.2% (26) and fracture toughness of ∼1 MPa ⋅ m^1/2, and thus it is quite brittle (27). Considerable research efforts have been devoted to enhancing its electric field–induced strain via 90° ferroelectric domain switching (28–30). However, the process is often irreversible, and external mechanisms such as restoring force (28, 29) and internal mechanisms such as defect pinning (30) have to be invoked to make the electrostrain recoverable. Nevertheless, it hints at the possibility of BTO being made superelastic by taking advantage of the stress-induced 90° domain switching (6). Earlier works suggest that surface tension induces an in-plane compressive stress that favors the axial polarization in one-dimensional ferroelectrics at small size (31, 32), which may provide the necessary restoring mechanism for the stress-switched domains. Thus, if a compressive axial force is applied, reversible domain switching may occur during unloading, leading to superelasticity. To verify this hypothesis, we fabricated single-crystalline BTO micropillars from [001]-oriented bulk crystals (SI Appendix, Fig. S1A) via focused ion beam (FIB), as detailed in Materials and Methods and SI Appendix, Fig. S1B. The diameters (Φ) of the micropillars range from 0.5 μm to 5 μm, with their height to diameter ratio fixed at 3. No visible defects can be seen from the scanning electron microscopy (SEM) images of these micropillars shown in Fig. 1 A–D, and their surfaces appear to be quite smooth, suggesting that no apparent damages are induced by FIB.Open in a separate windowFig. 1.Superelastic BTO micropillars below a critical size. (A–D) SEM images of the micropillars with Φ = 5, 3, 2, and 0.5 μm. (E–G) The first and second cycles of stress–strain curves for BTO micropillars with Φ = 5, 2, and 0.5 μm. (H) S_r/S_max and ΔW/W_max during the first cycle for BTO micropillars of different diameters. Here, S_r and S_max denote the residual strain and the maximum strain (SI Appendix, Fig. S6A), while ΔW and W_max are energy dissipated and stored in the first cycle, respectively (SI Appendix, Fig. S6F).     

Staphylococcal protein A inhibits complement activation by interfering with IgG hexamer formation

Immunoglobulin (Ig) G molecules are essential players in the human immune response against bacterial infections. An important effector of IgG-dependent immunity is the induction of complement activation, a reaction that triggers a variety of responses that help kill bacteria. Antibody-dependent complement activation is promoted by the organization of target-bound IgGs into hexamers that are held together via noncovalent Fc-Fc interactions. Here we show that staphylococcal protein A (SpA), an important virulence factor and vaccine candidate of Staphylococcus aureus, effectively blocks IgG hexamerization and subsequent complement activation. Using native mass spectrometry and high-speed atomic force microscopy, we demonstrate that SpA blocks IgG hexamerization through competitive binding to the Fc-Fc interaction interface on IgG monomers. In concordance, we show that SpA interferes with the formation of (IgG)₆:C1q complexes and prevents downstream complement activation on the surface of S. aureus. Finally, we demonstrate that IgG3 antibodies against S. aureus can potently induce complement activation and opsonophagocytic killing even in the presence of SpA. Together, our findings identify SpA as an immune evasion protein that specifically blocks IgG hexamerization.

Antibodies play a key role in the human immune response against bacterial infections. While antibodies can bind and neutralize bacterial virulence factors, they can also signal to components of the innate immune system and induce bacterial killing. To do so, antibodies bind bacterial cells via their variable (Fab) region and subsequently trigger Fc-mediated effector functions (1). The complement system, a large network of plasma proteins, forms an important effector of antibody-dependent immune protection against invading bacteria. An activated complement cascade results in efficient decoration of bacteria with C3-derived molecules that are essential to trigger highly effective phagocytic uptake via complement receptors on phagocytes. Furthermore, complement generates chemoattractants and induces direct killing of gram-negative bacteria. Because effective complement activation is an important effector mechanism of therapeutic antibodies in cancer (2), the ability of complement to kill bacteria could also be exploited for antibacterial therapies against (antibiotic-resistant) pathogens (3–5).The antibody-driven, “classical” complement pathway is initiated when circulating C1 complexes are recruited to antibody-labeled target surfaces (6). The most abundant antibody isotype in serum is immunoglobulin (Ig) G, which is subdivided into subclasses IgG1, IgG2, IgG3, and IgG4 in order of decreasing abundance. IgG antibodies can bind surface antigens via their Fab regions and subsequently recruit C1 via their Fc region (SI Appendix, Fig. S1A). The C1 complex consists of three large units: C1q, C1r, and C1s. C1q comprises the antibody recognition unit of the C1 complex and is composed of six globular heads connected by collagen-like stalks. On binding of C1q, its associated proteases C1r and C1s are activated to cleave other complement proteins that together form enzymes on the surface that catalyze the covalent deposition of C3b molecules onto the bacterial surface (SI Appendix, Fig. S1A). C3b molecules are recognized by complement receptors on phagocytes (neutrophils, macrophages), which engulf and digest bacteria intracellularly. The deposition of C3b also results in amplification of the complement cascade and activation of downstream complement effector functions.In recent years, it has become clear that efficient binding of C1 to target-bound IgG molecules requires IgGs to form ordered hexameric ring structures (7, 8). Cryo-electron tomography and atomic force microscopy studies revealed that the six globular heads of C1q can simultaneously bind to each of the six IgG molecules that form a hexameric binding platform (7) (SI Appendix, Fig. S1B). The formation of these hexamers is induced on antibody binding to surface-bound antigens and driven by noncovalent interactions between the Fc regions of neighboring IgG molecules (9) (Fig. 1A).Open in a separate windowFig. 1.The Ig-binding domains of SpA bind to residues of the IgG-Fc region that are involved in IgG hexamerization. (A) IgG hexamer crystal packing of IgG1-b12 (Protein Data Bank [PDB] ID code 1HZH). A single IgG is depicted in gray, and the IgG-Fc domain is enclosed in the dashed box. (B) Schematic representation of SpA organization. SpA consists of a signal sequence, five Ig-binding domains (E, D, A, B, and C), an Xr region (octapeptide repeats variable in number), and a cell wall attachment and sorting region that includes a constant Xc region, the LPETG motif, a hydrophobic anchor, and positively charged residues. (C) Sequence alignment of the five highly homologous Ig-binding domains of SpA. The amino acid residues conserved in all five domains are highlighted in green. The residues involved in the interaction with the Fc region of IgG are shown in pink. (D) Space-filling presentation depicting the Fc domain of IgG1-b12 and its interaction with SpA-B (PDB ID code 1FC2; the complementary Fc docking domain for SpA-B is hidden). The residues involved in Fc-Fc interactions required to form the IgG hexameric ring are depicted in green, and the crystal structure of SpA-B is shown in red.Interestingly, some bacteria produce IgG-binding molecules that recognize the Fc domain of IgGs (10). The best known of these is staphylococcal protein A (SpA), a 42-kDa protein that has a high affinity for the Fc region of IgG and thus is commonly used as a tool in affinity chromatography to purify monoclonal antibodies. SpA is produced by Staphylococcus aureus, an important human pathogen that is the main cause of serious hospital-acquired infections, such as bacteremia, sepsis, and endocarditis (11). Due to the dramatic increase in antibiotic resistance and the lack of proper vaccines, physicians are frequently left with no useful or suboptimal alternatives when treating these infections.SpA is considered an important virulence factor (12, 13) and vaccine candidate (14, 15). The protein is abundantly present on the bacterial cell wall (16–18) but is also released in the extracellular environment (19, 20). SpA is composed of a signal sequence, five sequential Ig-binding domains (denoted E, D, A, B, and C), an Xr region and a cell wall attachment and sorting region (21) (Fig. 1B). Each of the five repeating Ig-binding domains adopts a three-helical structure that can bind to the Fc region of IgG via helices I and II (22) and the Fab region of the VH3 type family of antibodies via helices II and III (23–25). The binding of SpA to Ig-Fc regions is considered to protect S. aureus from phagocytic killing (14) while cross-linking of Ig-Fab regions triggers the proliferation and apoptotic collapse of B cells (26). The five Ig-binding domains are highly homologous, sharing 74 to 91% of their amino acid sequence relative to the A domain (27) (Fig. 1C). It has been demonstrated that the binding interface of the B domain of SpA (SpA-B) and IgG1-Fc involves 11 amino acid residues of SpA-B and 9 residues of IgG1-Fc (22). Interestingly, SpA binds to all IgG subclasses except IgG3, due to a substitution in one of the nine Fc-contact residues in IgG3 (His⁴³⁵ in IgG1 becomes Arg⁴³⁵ in IgG3) (28). The residue Arg⁴³⁵ has been suggested to cause steric hindrance to SpA when binding to the IgG3-Fc (29). The reported crystal structure of SpA-B and its IgG-Fc interaction site are depicted in Fig. 1D. Of note, SpA binds to the same interface where the IgG Fc-Fc interactions take place to form the hexameric IgG platform required for complement activation.Here we investigated the impact of IgG-Fc binding properties of SpA on the assembly of IgG molecules into hexamers both in solution and on antigenic surfaces. We show that SpA blocks the formation of the hexameric C1q binding platform and as a result inhibits IgG-dependent complement activation and opsonophagocytic killing (OPK) of S. aureus. Our data provide an important contribution to the understanding of molecular mechanisms of complement evasion, which is crucial for the intelligent design of new therapeutic strategies to tackle infectious diseases.     

Persistent polyamorphism in the chiton tooth: From a new biomineral to inks for additive manufacturing

Engineering structures that bridge between elements with disparate mechanical properties are a significant challenge. Organisms reap synergy by creating complex shapes that are intricately graded. For instance, the wear-resistant cusp of the chiton radula tooth works in concert with progressively softer microarchitectural units as the mollusk grazes on and erodes rock. Herein, we focus on the stylus that connects the ultrahard and stiff tooth head to the flexible radula membrane. Using techniques that are especially suited to probe the rich chemistry of iron at high spatial resolution, in particular synchrotron Mössbauer and X-ray absorption spectroscopy, we find that the upper stylus of Cryptochiton stelleri is in fact a mineralized tissue. Remarkably, the inorganic phase is nano disperse santabarbaraite, an amorphous ferric hydroxyphosphate that has not been observed as a biomineral. The presence of two persistent polyamorphic phases, amorphous ferric phosphate and santabarbaraite, in close proximity, is a unique aspect that demonstrates the level of control over phase transformations in C. stelleri dentition. The stylus is a highly graded material in that its mineral content and mechanical properties vary by a factor of 3 to 8 over distances of a few hundred micrometers, seamlessly bridging between the soft radula and the hard tooth head. The use of amorphous phases that are low in iron and high in water content may be key to increasing the specific strength of the stylus. Finally, we show that we can distill these insights into design criteria for inks for additive manufacturing of highly tunable chitosan-based composites.

Biominerals are broadly used by organisms to reinforce structural materials, enabling for instance locomotion, feeding, and defense but also finding application in sensing (1). A defining principle of mineralized tissues is their composite nature, reaping synergy from the combination of a soft macromolecular matrix and a hard, inorganic mineral phase (2, 3). Organisms functionally grade such composites by precisely controlling the phase, size, shape, orientation, dispersion, and spatial distribution of mineral nanoparticles. Harnessing the biological capability to create composites that combine complex shape with mechanical properties that are graded intricately yet over several orders of magnitude in range is of interest for a broad range of functional materials—for instance, for soft robotics (4).Chitons, a class of marine mollusks, are best known for the extreme hardness, rich chemistry, and intricate phase assemblage of their radula teeth (Fig. 1 and SI Appendix, Fig. S1, for recent reviews see refs. 5 and 6; for a cross-cutting review that includes chiton biominerals, see ref. 7). However, the mechanical system of their dentition is not only significantly more complex but also based on a continuous organic phase that is differentially reinforced. The chiton therefore serves as an excellent model system to study biological mechanisms and design principles.Open in a separate windowFig. 1.Radula teeth of C. stelleri. (A) Ventral aspect of C. stelleri [Image credit: Linda Schroeder (photographer)]. (B) Mouth and protruding anterior end of radula (ra) (8). Image credit: VicHigh Marine/David Young. (C) Mosaic image of the entire radula showing all stages of development, including deposition of the organic scaffold (stage I), infiltration of the cusp with ferrihydrite (stage II), conversion to magnetite (stage III), mineralization of the core (stage IV), and mature teeth (stage V). (D) SEM image of the anterior end of the radula with mature teeth. Major lateral teeth consist of the tricuspidate head (hd, 400 × 350 × 150 µm³) and the stylus (st, L-shape 1,400 × 1,400 × 350 µm³). The stylus anchors teeth on the thin (∼100 µm; SI Appendix, Fig. S10C) and flexible radula membrane (rm). The stylus canal (stc) runs along the length of the stylus but terminates below the head. (E) Rendering of a virtual section of a tooth head (hd) and upper stylus (st) generated from a 3D reconstruction of the normalized linear attenuation coefficient (LAC) as determined by synchrotron microcomputed tomography. Note the LAC is highest for the outer magnetite layer (ml) of the head, intermediate for AFP-based composite of the core (co), and rather low in the stylus. Typical for Cryptochiton type teeth, the core is exposed in a window (wi) in the magnetite layer on the trailing (anterior) face of the tooth.The chiton radula is slender ribbon with transverse rows of teeth (Fig. 1 B–D) (9). Two outsize major lateral teeth make contact with the substrate during grazing. Major lateral teeth consist of a tooth head and stylus (Fig. 1 D and E and SI Appendix, Figs. S1 and S10) that is anchored on the radula membrane (10, 11). The stylus canal runs along the length of the hollow stylus and terminates in a dead end below the junction zone between stylus and cusp. During the feeding stroke, the radula slides over a curved supporting surface and bends in two orthogonal directions. This results in a characteristic scraping and sweeping motion of radula teeth (12). The loss of entire rows of teeth due to wear is compensated for by synthesis of new teeth at the posterior end (13, 14). Newly formed teeth mature in several stages as they are transported toward the anterior end by the radula membrane (15). As a result, their entire development can be observed in one animal (Fig. 1C).The function of the radula requires highly disparate material properties. The tooth complex is based on a continuous organic scaffold comprised of semicrystalline, partially deacetylated α-chitin and protein and may be crosslinked by tanning reactions (16). The tooth head is comprised of a highly mineralized cusp with exceptional hardness, wear resistance, and self-sharpening properties (17). The cusp is supported by a softer core. In Cryptochiton-type teeth, magnetite covers the entire posterior surface (leading edge) of the cusp and all but a rectangular window on the anterior surface (trailing edge) (5). The biomineral of the tooth core, amorphous ferric phosphate (AFP), is thus exposed in the window (Fig. 1E and SI Appendix, Fig. S1) (18). Use of at least two biominerals is typical for chiton teeth, and a considerable number has been identified in the core of different species (see ref. 5 and references therein). Additionally, ferrihydrite occurs as a transient precursor phase in stage II of radula development (Fig. 1C) (18).The entire head is mounted on the stylus, an L-shaped chitinous tissue that is integral to the complex movement of the tooth head during the feeding stroke (Fig. 1D and SI Appendix, Fig. S10 A and B) (19). The stylus connects the tooth head to the radula and orients the tooth with respect to the substratum (20). The radula membrane (Fig. 1D and SI Appendix, Fig. S10C) has lower stiffness to accommodate complex shape changes but at the same time must be strong, tough, and resistant to fatigue to survive the cycling bending and unbending during feeding. As a consequence of these vastly different requirements, we expect that mechanical properties are strongly graded. This is well documented for the cusp and core of radula teeth but less explored for the stylus and radula (17, 21, 22).The stylus and the radula membrane are generally referred to as unmineralized tissues, even though the presence of transient mineral has been suggested and the junction zone does mineralize in later stages (5, 22). However, the chemical form in which iron appears in the junction zone and stylus remains unclear. We therefore set out to map the redox state and chemical environment of iron in the tooth of Cryptochiton stelleri with the long-term goal to trace these parameters over the development of the radula and thus gain insight into the mechanisms at play. This would then provide a foothold on the way to designing bio-inspired syntheses.Given the complex shape (Fig. 1 D and E) and small size of chiton teeth and the extraordinarily rich chemistry of iron, mapping of multiple, often poorly crystalline phases requires techniques that combine high spatial resolution with sensitivity for subtle differences in oxidation state and coordination geometry of iron. Synchrotron Mössbauer spectroscopy (SMS) recently emerged as a powerful tool that combines high spatial resolution with the deep chemical insights offered by classical Mössbauer spectroscopy (23–25). Herein, we report on our discovery of a biomineral in the mature upper stylus of C. stelleri using SMS and correlative imaging and spectroscopy techniques.     

The number of catalytic cycles in an enzyme’s lifetime and why it matters to metabolic engineering

Metabolic engineering uses enzymes as parts to build biosystems for specified tasks. Although a part’s working life and failure modes are key engineering performance indicators, this is not yet so in metabolic engineering because it is not known how long enzymes remain functional in vivo or whether cumulative deterioration (wear-out), sudden random failure, or other causes drive replacement. Consequently, enzymes cannot be engineered to extend life and cut the high energy costs of replacement. Guided by catalyst engineering, we adopted catalytic cycles until replacement (CCR) as a metric for enzyme functional life span in vivo. CCR is the number of catalytic cycles that an enzyme mediates in vivo before failure or replacement, i.e., metabolic flux rate/protein turnover rate. We used estimated fluxes and measured protein turnover rates to calculate CCRs for ∼100–200 enzymes each from Lactococcus lactis, yeast, and Arabidopsis. CCRs in these organisms had similar ranges (<10³ to >10⁷) but different median values (3–4 × 10⁴ in L. lactis and yeast versus 4 × 10⁵ in Arabidopsis). In all organisms, enzymes whose substrates, products, or mechanisms can attack reactive amino acid residues had significantly lower median CCR values than other enzymes. Taken with literature on mechanism-based inactivation, the latter finding supports the proposal that 1) random active-site damage by reaction chemistry is an important cause of enzyme failure, and 2) reactive noncatalytic residues in the active-site region are likely contributors to damage susceptibility. Enzyme engineering to raise CCRs and lower replacement costs may thus be both beneficial and feasible.

As the synthetic biology revolution brings engineering principles and practices into the life sciences, biomolecules are being rethought as component parts that are used to build new biosystems and improve existing ones (1–3). Enzymes—the working parts of metabolic systems—are targets for this rethinking and are increasingly being repurposed by rational design and directed evolution (4).Substrate specificity, catalytic efficiency, and expression level are common performance specifications for enzyme parts in metabolic engineering, but life span is not, despite its centrality in other engineering fields. Knowing an engineering component’s life span (how long it lasts in service) is critical to preventing system failures and optimizing maintenance schedules (5). Failure metrics such as “mean time to failure” (6) are consequently used widely in engineering, which distinguishes three types of failures: early, wear-out, and random or stochastic. All three have counterparts in enzymes operating in vivo (Fig. 1A) (7–18), but wear-out and random failures (Fig. 1A, red font) are most relevant to length of working life.Open in a separate windowFig. 1.The engineering concept of component failure and its application to enzymes in vivo. (A) The types of failure in manufactured components and their counterparts in enzymes operating in vivo. (B) Schematic representation of the time dependence of the hazard rate and the cumulative probability (increasing color density) that an individual component will have failed.In manufactured systems, wear-out failures are caused by cumulative deterioration processes or by use-dependent wear (Fig. 1A). Like all proteins, enzymes are subject to cumulative deterioration from oxidation, racemization, or other chemical events (“protein fatigue”) that can affect any part of the molecule and degrade its function (9–11). However, use-dependent wear-out has no equivalent in enzymes, i.e., enzyme performance is not progressively degraded by operation of the catalytic cycle in the way a bearing is worn down a little each time it turns (Fig. 1A). Rather, a random catalytic misfire or a chemical attack by a substrate or product on a vulnerable residue in the active-site region can instantly inactivate an enzyme, whatever its age (14–18). Such failures thus have a constant hazard rate and are random or stochastic, like the abrupt failure of a transistor due to a current surge (Fig. 1A).Although the hazard of random failure does not depend on a part’s age, the cumulative probability that any individual part will experience a random failure increases with time (Fig. 1B). Given long enough, certain types of enzyme molecule may thus be doomed to have a terminal, catalysis-related accident. Such self-inflicted inactivation processes are important considerations for industrial enzymes (i.e., enzymes used ex vivo as reagents) and the number of catalytic cycles that each enzyme molecule carries out in its lifetime—often called “total turnover number”—is a key industrial performance criterion (19–21).The number of catalytic cycles mediated before self-inactivation could also be key to in vivo enzyme performance. Recent proteomic evidence points to damage from the reaction catalyzed as a major mode of enzyme failure and to the possibility that some reactions do more damage than others. Thus, in the bacterium Lactococcus lactis, a fivefold increase in growth rate was accompanied by a sevenfold increase in protein turnover rate (22). This near proportionality implies that L. lactis enzymes catalyze a similar number of reactions in their lifetimes, whatever the growth rate. This fits with reaction-related damage as a cause of failure: The faster the growth, the more flux through reactions, the more damage to enzymes, and the sooner enzymes fail. Similarly, protein turnover in yeast was faster when enzymes were in active use (23). Furthermore, in L. lactis, yeast, and Arabidopsis, the fastest turning-over metabolic enzymes include many with reactive substrates, products, or intermediates (SI Appendix, Table S1) (22–24), i.e., with a high risk of spontaneous chemical damage to the active site.The rates at which enzyme proteins are degraded and resynthesized are critical to the cellular energy economy because such turnover can consume about half the maintenance energy budget in microbes and plants (22, 25–27). High enzyme protein turnover rates therefore potentially reduce the productivity of biosystems ranging from microbial fermentations to crops (26, 28, 29). Consistent with such reduction, fast protein turnover is associated with low biomass yield in yeast (27) and with low growth rate in Arabidopsis (30). Also, slowing the turnover of abundant, fast-turnover enzymes is predicted to substantially increase growth rate and biomass yield in plants (26, 31) and other organisms (32).Rational design or directed evolution can now be used to tune protein turnover rates (33–35). However, before setting out to reduce enzyme turnover it is essential to define target enzymes and to understand why they turn over fast in the first place. Accordingly, here we calculate and compare the life spans of enzymes from three kingdoms using the criterion of “catalytic cycles until replacement” (CCR) (33), defined as the moles of substrate converted per mole of enzyme before the enzyme is replaced, i.e., the following:CCR=Metabolic flux rateEnzyme replacement rate.[1]CCR is the in vivo equivalent of the ex vivo “total turnover number” mentioned above but is a preferable term as it avoids confusion with the term “turnover number,” a synonym in enzymology for k_cat (20). CCR is envisioned as a potential constant, with reaction wear-and-tear being matched with degradation rates to maintain CCR as a factor hardwired to the structural and (bio)chemical stability of a given enzyme (33). We then compare each enzyme’s CCR to its reaction chemistry and across kingdoms to find shared attributes underlying CCR values. Our findings imply that CCRs are commonly influenced by random collateral damage from the reaction catalyzed and that enzymes could be engineered to reduce this damage and its attendant enzyme replacement costs. More generally, the findings point to catalysis-related accidents as a sizeable but underrecognized cause of enzyme failure and replacement.     

Selective constraints on global plankton dispersal

Marine microbial communities are highly interconnected assemblages of organisms shaped by ecological drift, natural selection, and dispersal. The relative strength of these forces determines how ecosystems respond to environmental gradients, how much diversity is resident in a community or population at any given time, and how populations reorganize and evolve in response to environmental perturbations. In this study, we introduce a globally resolved population–genetic ocean model in order to examine the interplay of dispersal, selection, and adaptive evolution and their effects on community assembly and global biogeography. We find that environmental selection places strong constraints on global dispersal, even in the face of extremely high assumed rates of adaptation. Changing the relative strengths of dispersal, selection, and adaptation has pronounced effects on community assembly in the model and suggests that barriers to dispersal play a key role in the structuring of marine communities, enhancing global biodiversity and the importance of local historical contingencies.

Ocean microbial biogeography is determined by the balance of two opposing forces: dispersal by the ocean currents and selection by the local environment (1). In the limit where global dispersal is fast relative to population turnover, environmental conditions alone should be sufficient to predict the presence or absence of a particular species from any given location on Earth (2, 3). This is the view encapsulated in the hypothesis of Baas-Becking (4) that “everything is everywhere, but the environment selects.” On the other hand, if global dispersal is slow relative to population turnover, limited connectivity between ocean regions will tend to reinforce chance differences between isolated communities (5, 6), with geographically isolated but otherwise similar environments displaying significant differences in taxonomic composition.Over evolutionary timescales, the balance of dispersal and selection will affect community assembly [through diversification and mass effects (7)], ecosystem function (through biogeochemical cycling), and ultimately, the resilience of marine ecosystems to environmental change (8). Therefore, understanding the mechanisms that lead to niche diversification and biogeographic structure in microbial communities is a fundamental pursuit of marine microbial research. A central question is to what degree are biogeographic patterns attributable to local selection based on contemporary environmental factors or to independent stochastic processes occurring in geographically isolated regions (SI Appendix, Fig. S1) (1).Recent analysis of metagenomic data (Fig. 1) (9) has shown that large-scale trends in community composition are correlated both with environmental variables and with geographic distance, with distinct clusters emerging along environmental gradients and among the most rapidly connected sites, suggesting that both history and environment play important roles. When sample sites are clustered based on metagenomic pairwise β-diversity (SI Appendix), there is discernible ecological similarity among sites within the same ocean basins (Fig. 1A), although we also see geographically proximate sites clustered far apart and sites from geographically remote locations clustered together (Fig. 1B and SI Appendix, Figs. S4 and S5). These broad patterns appear to reflect both geographic proximity and environmental selection (9). Nonetheless, it can be difficult to assign causal mechanisms, and the drivers of observed biogeography thus remain uncertain.Open in a separate windowFig. 1.Taxonomic community similarity clusters in the 0.22- to 3-μm size fraction across Tara Oceans sites (replotted using data from ref. 9). (A) Community similarity is shown with colors by projecting the Taxonomic Jaccard dissimilarity matrix into the “rgb” (red-green-blue) color space using the t-SNE (t-distributed stochastic neighbour embedding) dimension-reduction algorithm (10). (B) Links between community similarity clusters (dimensionless x and y coordinates) and spatial location (colors corresponding to ocean basins). SI Appendix, Fig. S1 has an interpretation of B.The roles of selection and dispersal have both been examined using global-scale models but typically, with one in isolation from the other. On one hand, population dynamic models have focused on the role of selection from among a universal background of candidate species (11), in line with the Baas-Becking (4) hypothesis. On the other hand, a number of studies have addressed the question of global gene flow in oceanic microbial communities, using particle tracking models to assess connectivity through the surface waters (6, 12), but these have typically assumed ecological neutrality (5) and have thus ignored the role of selection. While some studies find that the ocean surface is very rapidly connected on timescales of decades or less (12), others suggest that current rates of passive dispersal are insufficient to overcome biogeographic differences created by chance mutations occurring in geographically isolated regions of the ocean (6).In order to distinguish between the biogeographic effects of selection and dispersal, we need a framework that accounts for both processes together. In this paper, we develop a population genetic model representing taxonomic and phenotypic diversity within a single clonally reproducing plankton population, embedded within an empirically constrained representation of the ocean circulation (13). In contrast to previous studies, our model accounts for population size, stochastic demography, natural selection, adaptation, and transport through the ocean interior (we find that dispersal pathways restricted to the ocean surface are artificially sensitive to fluid convergence and divergence). With a more realistic transport term accounting for dispersal at all depths, we find that varying the degree of selection and adaptation leads to very different model outcomes in terms of community biogeography and global connectivity. We show that selection based on thermal niches acts as a major constraint on dispersal, with the clear effects on biogeographic organization at the global scale.     

Selective cysteine-to-selenocysteine changes in a [NiFe]-hydrogenase confirm a special position for catalysis and oxygen tolerance

In [NiFe]-hydrogenases, the active-site Ni is coordinated by four cysteine-S ligands (Cys; C), two of which are bridging to the Fe(CO)(CN)₂ fragment. Substitution of a single Cys residue by selenocysteine (Sec; U) occurs occasionally in nature. Using a recent method for site-specific Sec incorporation into proteins, each of the four Ni-coordinating cysteine residues in the oxygen-tolerant Escherichia coli [NiFe]-hydrogenase-1 (Hyd-1) has been replaced by U to identify its importance for enzyme function. Steady-state solution activity of each Sec-substituted enzyme (on a per-milligram basis) is lowered, although this may reflect the unquantified presence of recalcitrant inactive/immature/misfolded forms. Protein film electrochemistry, however, reveals detailed kinetic data that are independent of absolute activities. Like native Hyd-1, the variants have low apparent K_MH₂ values, do not produce H₂ at pH 6, and display the same onset overpotential for H₂ oxidation. Mechanistically important differences were identified for the C576U variant bearing the equivalent replacement found in native [NiFeSe]-hydrogenases, its extreme O₂ tolerance (apparent K_MH₂ and V_max [solution] values relative to native Hyd-1 of 0.13 and 0.04, respectively) implying the importance of a selenium atom in the position cis to the site where exogenous ligands (H⁻, H₂, O₂) bind. Observation of the same unusual electrocatalytic signature seen earlier for the proton transfer-defective E28Q variant highlights the direct role of the chalcogen atom (S/Se) at position 576 close to E28, with the caveat that Se is less effective than S in facilitating proton transfer away from the Ni during H₂ oxidation by this enzyme.

Hydrogenases catalyze highly efficient H₂ activation, providing a paradigm for renewable hydrogen technologies (1). In a small subgroup of [NiFe]-hydrogenases from sulfate-reducing bacteria and methanogens, natural substitution of cysteine (Cys; C) for selenocysteine (Sec; U) occurs in the active site (Fig. 1) (2). The [NiFeSe]-hydrogenases (group 1a) are reported to have higher activity than their [NiFe] counterparts—a feature seen in other enzymes where C and U are swapped (2–8). Escherichia coli produces [NiFe]-hydrogenase-1 (Hyd-1) (group 1d, O₂-tolerant) and Hyd-2 (group 1c, O₂-sensitive) membrane-bound [NiFe]-hydrogenases (3, 9–12). At neutral pH in vitro Hyd-1 performs H₂ oxidation only, whereas Hyd-2 can also produce H₂ (reduce H⁺) (13). Hydrogen oxidation activity in vivo is linked to reduction of different terminal electron acceptors depending on the bacterial species, availability of different oxidants, and their redox potential, for example fumarate or, in the case of the Knallgas bacterium Ralstonia eutropha, O₂. The production of Hyd-1 and Hyd-2 is maximal using fumarate as the terminal electron acceptor under anaerobic conditions (13). E. coli does not produce a [NiFeSe]-hydrogenase; Hyd-3 (group 4a) is U-containing in the formate dehydrogenase (FdhF) subunit only (3, 14).Open in a separate windowFig. 1.(A) Amino acid alignment of selected hydrogenases (Hyd-1 numbering) highlighting key residues (Cys/Sec, red; E28, green; D118, pink; R509, yellow) and differences (cyan). See also SI Appendix, Fig. S1 and Table S1. (B and D) The extended active site of Hyd-1 (Protein Data Bank [PDB] ID code 5A4M) (B) and Desulfomicrobium baculatum NiFeSe (PDB ID code 1CC1) (D) hydrogenases. (C) Representation of the active site, where “X” denotes the atom in a bridging position between the Ni and Fe atoms, the identity of which depends on the (in)active state of the enzyme (SI Appendix, Fig. S3).Minimally, Hyd-1 has two membrane-extrinsic subunits: HyaB containing the NiFe active site, and HyaA housing three FeS clusters to mediate long-range electron transfer (SI Appendix, Figs. S1 and S2) (15). The resulting complex, a (HyaA)₂(HyaB)₂ dimer, transfers electrons to a b-type cytochrome in a membrane-intrinsic HyaC subunit (10). For periplasmic [NiFeSe]-hydrogenases from sulfate-reducing bacteria, the normal redox partner is a soluble cytochrome c₃ (16). In protein film electrochemistry (PFE; see below), the FeS clusters connect the active site to an electrode, enabling catalysis to be controlled and recorded (1).The active-site metal atoms are coordinated by four conserved Cys residues (Fig. 1), two of which are terminal to the Ni and two of which are bridging (μ) between the Ni and Fe atoms. In [NiFeSe] homologs it is usually a terminal Cys (C576; Fig. 1B) that is replaced by Sec (Fig. 1D), although purported examples exist in which a bridging Cys residue C579 is substituted (2). Additionally, a nearby aspartate is substituted by serine in the active-site “canopy” of [NiFeSe]-hydrogenases (17). Other important residues include a strictly conserved arginine (R509) essential for fast and efficient H₂ oxidation in Hyd-1 (4, 17), and a glutamate (E28) adjacent to C576 which appears to be a universal proton gate (18, 19).The FeS relays are disparate (SI Appendix, Fig. S2): [NiFeSe]-hydrogenases coordinate three [4Fe-4S] clusters at all positions proximal, medial, and distal to the active site, whereas the medial site in [NiFe]-hydrogenases is a [3Fe-4S] cluster having a more positive reduction potential (20–23). In group 1d O₂-tolerant hydrogenases, such as Hyd-1, all the FeS clusters have more positive reduction potentials (24, 25); importantly, the proximal cluster is a unique [4Fe-3S]_6-Cys center that is essential for O₂ tolerance (9), a property requiring the invading O₂ molecule to be reduced to harmless water (1, 26). The [4Fe-3S]_6-Cys proximal cluster can perform two one-electron transfers back to the active site upon O₂ exposure during H₂ oxidation, a process requiring substantial conformational change to form the “superoxidized state” (11, 24, 25, 27, 28). A third electron is available from the high-potential medial [3Fe-4S] cluster and a fourth stems from oxidation of the Ni (SI Appendix, Fig. S3). A truly O₂-tolerant [NiFe]-hydrogenase is thus also an oxidase (26). Although [NiFeSe]-hydrogenases are considered “O₂-tolerant” (29), this property, requiring reductive destruction of O₂, may be limited to H₂ evolution (30).Selenocysteine, the versatile 21st amino acid, appears in proteins from all domains of life (31). Sec is structurally similar to Cys, except the thiol is replaced by a selenol (Fig. 2A). Selenium and sulfur are chalcogens; thus U and C share certain chemical properties, but the electronic structures of S and Se differ sufficiently to give selenoproteins distinctive catalytic efficacies (6). The much lower pKa of selenol compared with thiol renders it fully ionized at physiological pH (32), selenoproteins are more resistant to irreversible oxidation than their C-containing homologs (5), and diselenide bonds are more stable to reduction than disulfide bonds (7, 33). Most natural selenoenzymes are oxidoreductases having an essential (for efficient catalysis) U active-site residue, and many have Cys homologs from which they evolved (34).Open in a separate windowFig. 2.(A) Chemical structures of Cys and Sec show the selenol moiety (red) to be the only difference. (B) EF-Tu–driven site-specific incorporation of Sec at a UAG codon. mRNA, messenger RNA. (C) Coomassie blue-stained denaturing sodium dodecyl sulfate-polyacrylamide gel of Hyd-1 variants (C76U, C79U, C576U, and C579U) shows high purity in each case, comprising HyaA (37 kDa) and HyaB (66 kDa) only. (D) Tandem mass spectra of C76U and C576U show Sec incorporation at the desired position in the designated peptide. Red lines correlate with the cleavage products depicted in the peptide sequence with an accuracy of 5 ppm. See also SI Appendix, Figs. S4–S8.Advances in genetic code expansion have provided tools for effective, site-specific UAG-programmed Sec insertion into recombinant proteins in E. coli. These in vivo methods (35–37) rely on elongation factor EF-Tu, thus bypassing the natural complex U-specific selenoprotein synthesis machinery programmed by UGA (Fig. 2B). Recently, one of these methods was used to replace active-site Cys residues with Sec in ribonucleotide reductase (38). This encouraged us to produce Cys-to-Sec variants of Hyd-1 in E. coli by the same strategy (35, 37).A recent paper described the consequences of replacing the Sec residue of a natural [NiFeSe]-hydrogenase with Cys, thereby retroengineering it to resemble a [NiFe]-hydrogenase (39). Here we report the opposite and complementary study, substituting each active-site Cys residue for Sec at all four coordination positions in Hyd-1. The resulting data highlight why one particular position has special significance.     

From the Cover: Ancient noeggerathialean reveals the seed plant sister group diversified alongside the primary seed plant radiation

Noeggerathiales are enigmatic plants that existed during Carboniferous and Permian times, ∼323 to 252 Mya. Although their morphology, diversity, and distribution are well known, their systematic affinity remained enigmatic because their anatomy was unknown. Here, we report from a 298-My-old volcanic ash deposit, an in situ, complete, anatomically preserved noeggerathialean. The plant resolves the group’s affinity and places it in a key evolutionary position within the seed plant sister group. Paratingia wuhaia sp. nov. is a small tree producing gymnospermous wood with a crown of pinnate, compound megaphyllous leaves and fertile shoots each with Ω-shaped vascular bundles. The heterosporous (containing both microspores and megaspores), bisporangiate fertile shoots appear cylindrical and cone-like, but their bilateral vasculature demonstrates that they are complex, three-dimensional sporophylls, representing leaf homologs that are unique to Noeggerathiales. The combination of heterospory and gymnospermous wood confirms that Paratingia, and thus the Noeggerathiales, are progymnosperms. Progymnosperms constitute the seed plant stem group, and Paratingia extends their range 60 My, to the end of the Permian. Cladistic analysis resolves the position of the Noeggerathiales as the most derived members of a heterosporous progymnosperm clade that are the seed plant sister group, altering our understanding of the relationships within the seed plant stem lineage and the transition from pteridophytic spore-based reproduction to the seed. Permian Noeggerathiales show that the heterosporous progymnosperm sister group to seed plants diversified alongside the primary radiation of seed plants for ∼110 My, independently evolving sophisticated cone-like fertile organs from modified leaves.

The origin of the seed in the Late Devonian, ∼365 Mya, represents a key innovation in land plant evolution. Seeds provided a fundamentally new reproductive strategy that overcame the limitations of free-sporing, pteridophytic reproduction and enabled colonization of drier habitats (1‒3). Progymnosperms, the evolutionary stem group leading to seed plants, display a mosaic of evolutionary characters combining free-sporing reproduction with production of secondary xylem (wood) through a bifacial vascular cambium characteristic of seed plants (2‒4). Although not representing a monophyletic evolutionary group, progymnosperms are important for our present understanding of the origin of the seed and represent intermediates between pteridophytes and seed plants (1, 2, 4). Progymnosperms include the basally divergent homosporous Aneurophytales, as well as the more derived heterosporous Archaeopteridales, Protopityales (1, 4), and the enigmatic Carboniferous cone Cecropsis (5). The transition to seed plants requires multiple character state transitions from known progymnosperm sister groups comprising either Archaeopteris (6, 7) or Archaeopteris + Cecropsis (e.g., ref. 8). For each of these cases, viable intermediates are absent from the fossil record.In contrast, Noeggerathiales (9) have at times been proposed as progymnosperms (9–12), but this has been controversial. Comprising ∼20 genera and 50 species, Noeggerathiales are known from the late Carboniferous–Permian (323 to 251 Ma) tropical floras in North America, Europe, and East Asia (10), where they are identified as a related group based on shared features of heterospory, adaxial sporangial attachment to “sporophylls,” longitudinal sporangial dehiscence, plagiotropic pinnule attachment to the rachis, and once-pinnate compound megaphylls (9–12). However, their systematic position remained uncertain because details of their stem anatomy were unknown. In the absence of this anatomical information, they have been postulated as close relatives of the progymnosperms, leptosporangiate ferns, sphenopsids, the extant fern Tmesipteris, or as a distinct class of their own (see ref. 9 for summary). Recently, leaves of Plagiozamites oblongifolius that were interpreted as noeggerathialean (9) have also been interpreted as cycads (13). Wang et al. (9) considered the Noeggerathiales to be putative progymnosperms, but as their specimens lacked wood, a more confident assignment was not possible. This uncertainty meant that even when considered as progymnosperms, their relationship to other progymnosperms and seed plants has remained speculative.Here, we report a species that we name Paratingia wuhaia J. Wang et al. sp. nov. and ascribe to the Noeggerathiales. The fossils were collected from a single 66-cm-thick volcanic ash bed in the Chinese “vegetational Pompeii” from the Taiyuan Formation at Wuda open coalmine, Inner Mongolia (14, 15) (SI Appendix, Geological Information, Materials). The ash preserved in situ the morphology and anatomy of plants in exquisite detail (Fig. 1 A–N) and has been dated to 298.34 ± 0.09 Ma during the Asselian stage of the early Permian (16).Open in a separate windowFig. 1.P. wuhaia J. Wang et al. sp. nov. from the early Permian Taiyuan Formation of Wuda Coalfield, Inner Mongolia. (A) Holotype with an entire crown consisting of pseudostrobili and leaves. (B) Once-pinnate compound leaf with both large and small pinnules visible. Reprinted with permission from ref. 15. (C) Cross-section of a crown illustrating pseudostrobili around the stem. (D) Cross-section of pseudostrobilus with microsporangia around the axis with bilateral, inversed Ω-shaped vascular bundle. (E) Cross-section of a leaf rachis showing the same form of vascular bundle as that of pseudostrobili axes. (F–H) Partial cross, radial, and tangential sections of the stem showing the secondary xylem (wood). (I) Tangential section of pseudostrobilus showing sporangial arrangement with single line of megasporangia along with the axis (middle-right pseudostrobilus of the fragmental crown illustrated SI Appendix, Fig. S8; specimen PB22132; magnified in SI Appendix, Fig. S21C.) (J) Radial section of pseudostrobilus showing adaxial sporangia and axis lacking nodes (upmost-left pseudostrobilus of the fragmental crown in SI Appendix, Fig. S8; magnified in SI Appendix, Fig. S21A). (K) Tangential section showing adaxial sporangia and a single line of megasporangia along with the axis (SI Appendix, Fig. S7; specimen PB22131, second up-right pseudostrobilus). (L) Tangential section through same specimen as K showing megasporangial arrangement (magnified in SI Appendix, Figs. S21B and S23A). (M) Detail of the middle part of L showing the megasporangia and microsporangia. (N) Single spore macerated from the holotype. (Scale bars: A, 10 cm; B, 3 cm; C–E, 1 cm; F, 100 μm; G and H, 200 μm; I–L, 5 mm; M, 2 mm; N, 10 μm.)Isolated noeggerathialean fertile shoots (SI Appendix, Figs. S3–S8) and megaphyllous leaves (SI Appendix, Figs. S9‒S14) are common in the assemblage. The species P. wuhaia is based on >200 specimens, including >10 specimens that were preserved with intact crowns containing both leaves and fertile shoots in organic attachment to the stems, demonstrating that they belong to a single species (Fig. 1A and SI Appendix, Figs S2‒S4).     

The evolution and changing ecology of the African hominid oral microbiome

The oral microbiome plays key roles in human biology, health, and disease, but little is known about the global diversity, variation, or evolution of this microbial community. To better understand the evolution and changing ecology of the human oral microbiome, we analyzed 124 dental biofilm metagenomes from humans, including Neanderthals and Late Pleistocene to present-day modern humans, chimpanzees, and gorillas, as well as New World howler monkeys for comparison. We find that a core microbiome of primarily biofilm structural taxa has been maintained throughout African hominid evolution, and these microbial groups are also shared with howler monkeys, suggesting that they have been important oral members since before the catarrhine–platyrrhine split ca. 40 Mya. However, community structure and individual microbial phylogenies do not closely reflect host relationships, and the dental biofilms of Homo and chimpanzees are distinguished by major taxonomic and functional differences. Reconstructing oral metagenomes from up to 100 thousand years ago, we show that the microbial profiles of both Neanderthals and modern humans are highly similar, sharing functional adaptations in nutrient metabolism. These include an apparent Homo-specific acquisition of salivary amylase-binding capability by oral streptococci, suggesting microbial coadaptation with host diet. We additionally find evidence of shared genetic diversity in the oral bacteria of Neanderthal and Upper Paleolithic modern humans that is not observed in later modern human populations. Differences in the oral microbiomes of African hominids provide insights into human evolution, the ancestral state of the human microbiome, and a temporal framework for understanding microbial health and disease.

The oral cavity is colonized by one of the most diverse sets of microbial communities of the human body, currently estimated at over 600 prevalent taxa (1). Dental diseases, such as caries and periodontitis, remain health burdens in all human populations despite hygiene interventions (2, 3), and oral microbes are often implicated in extraoral inflammatory diseases (4, 5). To date, most oral microbiome research has focused on clinical samples obtained from industrialized populations that have daily oral hygiene routines and access to antibiotics (1, 6), but far less is known about the global diversity of the oral microbiome, especially from diverse past and present nonindustrialized societies (7). The oral cavity contains at least six distinct habitats, but dental biofilms, including both supra- and subgingival dental plaque, are among the most diverse and clinically important (1, 6, 8). During life, these dental biofilms naturally and repeatedly calcify, forming dental calculus (tooth tartar) (9), a robust, long-term record of the oral microbiome (10). Archaeological dental calculus has been shown to preserve authentic oral bacterial metagenomes in a wide range of historic and prehistoric populations and up to 50 thousand years ago (ka) (10–13). As such, dental calculus presents an opportunity to directly investigate the evolution of the hominid microbiome and to reconstruct ancestral states of the modern human oral microbiome. In addition, because research has shown that evolutionary traits, diet, and cultural behaviors shape modern human microbiome structure and function at other body sites, such as the gut and skin microbiomes (14–18), investigating ancient oral metagenomes has the potential to reveal valuable information about major events in modern human evolution and prehistory, such as predicted dietary changes during the speciation of Homo (19–21) and the direct interaction of Neanderthals and modern humans during the Late Pleistocene (22).To better understand the evolutionary ecology of the African hominid microbiome, we generated and analyzed 109 dental calculus metagenomes from present-day modern humans (n = 8), gorillas (Gorilla, n = 29), chimpanzees (Pan, n = 20), Neanderthals (n = 13), and two groups of archaeological modern humans associated with major lifestyle transitions (preagricultural, n = 20; preantibiotic, n = 14), as well as New World howler monkeys (n = 5) for comparison (SI Appendix, Fig. S1). To account for potential sampling biases, we analyzed multiple subspecies and populations of each African great ape genus, which were obtained from C20^th or C21^st-collected museum collections, and for modern humans we sampled multiple populations from both Africa and Europe. To this, we added previously published microbiome data from chimpanzees (n = 1) (13), Neanderthals (n = 4) (13), and present-day modern humans (n = 10) (23), for a total dataset of 124 individuals (Fig. 1A, SI Appendix, Table S1, and Dataset S1). We also generated eight new radiocarbon dates for archaeological individuals, for a total of 44 directly or indirectly dated ancient individuals in this study (Dataset S1).Open in a separate windowFig. 1.Sample locations and oral microbiome authentication of ancient dental calculus. (A) Sample locations. (B) PCoA comparing euclidean distances of microbial genera of well-preserved ancient and present-day dental calculus to environmental proxy controls (degraded archaeological bone) and present-day dental plaque and feces. Ancient dental calculus is distinct from gut and archaeological bone but overlaps with present-day dental plaque. (C) Representative DNA damage patterns for Neanderthals and ancient and present-day modern humans for four oral-specific bacterial species. The Neanderthal and upper Paleolithic modern human individuals show expected damage patterns consistent with authentic aDNA, whereas the present-day individual does not. See also SI Appendix, Fig. S4.Here, we investigate the structure, function, and core microbial members of the human oral microbiome within an evolutionary framework, seeking to determine whether a core microbiome can be defined for each African hominid group, whether the core is phylogenetically coherent, and whether some members of the core are specific to certain host groups. We test whether the oral microbiome of hominids reflects host phylogeny, finding that African hominid oral microbiota are distinguished by major taxonomic and functional differences that only weakly reflect host relationships and are likely influenced by other physiological, dietary, or behavioral factors. We compare the microbial profiles of Neanderthals and modern humans and, contrary to expectations (12, 13), find a high consistency of oral microbiome structure within Homo, regardless of geography, time period, or diet/lifestyle. We detect the persistence of shared genetic diversity in core taxa between Neanderthals and Upper Paleolithic humans prior to 14 ka, supporting a growing body of evidence for earlier admixture and interaction in Ice Age Europe (24, 25). Finally, we explore possible implications of our findings on Homo-associated encephalization (19, 26) and the role of dietary starch in human evolution (20, 21) by investigating the evolutionary history of amylase-binding capability by oral streptococci. We find that amylase binding is an apparent Homo-specific trait, suggestive of microbial coadaptation to starch-rich diets early in human evolution.     

A new sea-level record for the Neogene/Quaternary boundary reveals transition to a more stable East Antarctic Ice Sheet

Sea-level rise resulting from the instability of polar continental ice sheets represents a major socioeconomic hazard arising from anthropogenic warming, but the response of the largest component of Earth’s cryosphere, the East Antarctic Ice Sheet (EAIS), to global warming is poorly understood. Here we present a detailed record of North Atlantic deep-ocean temperature, global sea-level, and ice-volume change for ∼2.75 to 2.4 Ma ago, when atmospheric partial pressure of carbon dioxide (pCO₂) ranged from present-day (>400 parts per million volume, ppmv) to preindustrial (<280 ppmv) values. Our data reveal clear glacial–interglacial cycles in global ice volume and sea level largely driven by the growth and decay of ice sheets in the Northern Hemisphere. Yet, sea-level values during Marine Isotope Stage (MIS) 101 (∼2.55 Ma) also signal substantial melting of the EAIS, and peak sea levels during MIS G7 (∼2.75 Ma) and, perhaps, MIS G1 (∼2.63 Ma) are also suggestive of EAIS instability. During the succeeding glacial–interglacial cycles (MIS 100 to 95), sea levels were distinctly lower than before, strongly suggesting a link between greater stability of the EAIS and increased land-ice volumes in the Northern Hemisphere. We propose that lower sea levels driven by ice-sheet growth in the Northern Hemisphere decreased EAIS susceptibility to ocean melting. Our findings have implications for future EAIS vulnerability to a rapidly warming world.

The instability of polar continental ice sheets in a warmer future is an issue of major societal concern (1–5). Based on linear extrapolation of recent sea-level rise (2), mean global sea level could increase by 65 ± 12 cm by 2100 relative to the 2005 baseline, consistent with Intergovernmental Panel on Climate Change projections (1) of a ∼30- to 100-cm increase by 2100. Further, satellite observations (4) document substantial mass loss of both the Greenland Ice Sheet (GIS) and the West Antarctic Ice Sheet (WAIS) over the past decade—the two ice sheets that are most susceptible to global warming because of rapidly rising Arctic air temperatures (1) (GIS) and vulnerability to ocean-atmospheric warming (5, 6) (WAIS). The mass balance of the much larger EAIS and its contribution to ongoing sea-level change, however, remain poorly constrained (1).The role of atmospheric partial pressure of carbon dioxide (pCO₂) as a driver of long-term changes in ice volume and sea level over the Cenozoic Era (past ∼66 My) is widely documented (7–9) and there is compelling evidence (6, 10–12) of East Antarctic Ice Sheet (EAIS) retreat during warm intervals of the Pliocene epoch between ∼5.3 and 3.3 Ma when pCO₂ levels (13, 14) last reached values close to the present day (∼400 parts per million volume [ppmv]; Fig. 1 A and B and see SI Appendix, section S1). However, there is disagreement over EAIS behavior under pCO₂ levels (13) similar to those of preindustrial Quaternary times (<280 ppmv). A compilation of marine geochemical paleo-sea-level and pCO₂ records suggests that the EAIS was stable under these conditions (7). In contrast, while the amplitudes of change are controversial (15) (SI Appendix, section S2), sea-level reconstructions from paleoshorelines (16) and benthic geochemical data (9, 17, 18) (Fig. 2) imply EAIS melting during the Quaternary “super-interglacials” of Marine Isotope Stage (MIS) 11 (∼400 ka) and 31 (∼1.07 Ma) under relatively low pCO₂ conditions. Supporting evidence for EAIS retreat during the most recent “super-interglacial” MIS 11 comes from isotope measurements in mineral deposits recording past changes in subglacial East Antarctic waters (19), as well as records of ice-rafted debris (IRD) and detrital sediment neodymium isotopes from offshore the Wilkes Subglacial Basin (20). The latter records (20) also indicate EAIS retreat during the last interglacial MIS 5e (∼120 ka). Melting of the EAIS as inferred in the late Quaternary was likely driven by ocean–atmosphere warming around Antarctica and grounding-line retreat in response to ice–ocean interactions (19, 20).Open in a separate windowFig. 1.Neogene to Quaternary climate and sea-level evolution. (A) LR04 stack (21) for the past 5 My; arrow indicates the iNHG (∼3.6 to 2.4 Ma) and its culmination (thick-arrowed interval) (22); green line indicates the benthic δ¹⁸O level associated with MIS 101. (B) Atmospheric pCO₂ estimates of refs. 13 (blue) and 23 (purple) for the past 5 My; the late Quaternary glacial–interglacial pCO₂ range (1) is indicated as preindustrial pCO₂ band. Yellow shading in A and B highlights the study interval (∼2.75 to 2.4 Ma). (C and D) Site U1313 benthic δ¹⁸O and Mg/Ca raw data, respectively. (E) Site U1313 deep-sea temperature. (F) Site U1313 δ¹⁸O_sw-based sea level relative to present (black line); blue shading: 95% probability interval from Monte Carlo simulations (2σ); red line: threshold (11.6 msle) above which a smaller-than-present EAIS is signaled (24–26); m = marine part of EAIS, t = terrestrial part of EAIS. Glacials are highlighted in gray.Open in a separate windowFig. 2.Implication of different sea-level-δ¹⁸O_sw conversions for estimates of interglacial ice-volume loss. y axis shows lower-than-modern δ¹⁸O_sw values (∆δ¹⁸O_sw) and the x axis (log-scale) the corresponding sea-level increase for commonly used conversion factors (27–29) (0.011 [black], 0.010 [purple], and 0.008 ‰⋅m⁻¹ [red]) and those for Antarctica only (11) (0.014 ‰⋅m⁻¹) ignoring (yellow) and incorporating (brown) the impact of its marine-based ice sheets. Stars mark ∆δ¹⁸O_sw for interglacials of this study and corresponding sea-level equivalents in dependence of the conversion applied. Orange, blue, and purple diamonds show the same for MIS 31, 11 (18), and 5e (17), respectively. Vertical lines indicate the sea-level increase resulting from complete melting of the GIS (+7.3 m), WAIS (+4.3 m), and EAIS (+53.3 m) (24–26).To further investigate past EAIS response to climate forcing we studied the Neogene/Quaternary transition when mean pCO₂ (13, 23) fell from levels similar to the anthropogenically perturbed values of today into the Quaternary range, leading to progressive high-latitude cooling and the intensification of Northern Hemisphere Glaciation (21, 30–33) (iNHG; Fig. 1 A and B). Our approach is based on a simple approximation that, once estimated past global sea level exceeds 11.6 m sea-level equivalent (msle) above modern, which corresponds to the complete melting of the present-day GIS [7.3 msle (24, 25)] and the marine- and land-based WAIS [3.4 and 0.9 msle (25, 26), respectively], EAIS instability (i.e., a retreat from its present-day size) can be inferred (see SI Appendix, section S4.1 for details). We quantified sea-level and ice-volume changes for the interval ∼2.75 to 2.4 Ma (MIS G7 to 95) by measuring the oxygen-isotope composition (δ¹⁸O) and Mg/Ca ratio in well-preserved benthic foraminiferal calcite (Oridorsalis umbonatus) from Integrated Ocean Drilling Program (IODP) Site U1313 [41°0′N, 32°57′W; 3,426-m water depth (34)] in the North Atlantic Ocean (Fig. 1 C and D). Using this approach we reconstructed changes in seawater δ¹⁸O (δ¹⁸O_sw), a proxy for global sea level and continental ice volume (35). This was done by 1) calculating bottom-water temperatures (BWT) derived from Mg/Ca (36) (Fig. 1E), 2) combining Mg/Ca-derived BWTs with δ¹⁸O to determine δ¹⁸O_sw (37) (Fig. 1F), and 3) converting δ¹⁸O_sw to sea level using a relationship between changes in sea level and δ¹⁸O_sw of 0.011 ‰⋅m⁻¹ (27) (Materials and Methods and Fig. 1F). Ninety-five percent probability intervals calculated through Monte Carlo simulations for individual sea-level data points yield an average uncertainty for our sea-level estimates of ± 28 m (∼2σ [SD]) (Materials and Methods and Fig. 1F), roughly equivalent to the decay/growth of ice four times greater than the GIS. Our approach was validated by reconstructing δ¹⁸O_sw for the recent (∼0 to 7 ka) at IODP Site U1313 and for late Holocene core-top (multicorer) samples from a neighboring site (MSM58) which are indistinguishable from the observed modern-day values (see Materials and Methods and SI Appendix, section S4.2.8 for details).     

Genomic stability through time despite decades of exploitation in cod on both sides of the Atlantic

The mode and extent of rapid evolution and genomic change in response to human harvesting are key conservation issues. Although experiments and models have shown a high potential for both genetic and phenotypic change in response to fishing, empirical examples of genetic responses in wild populations are rare. Here, we compare whole-genome sequence data of Atlantic cod (Gadus morhua) that were collected before (early 20th century) and after (early 21st century) periods of intensive exploitation and rapid decline in the age of maturation from two geographically distinct populations in Newfoundland, Canada, and the northeast Arctic, Norway. Our temporal, genome-wide analyses of 346,290 loci show no substantial loss of genetic diversity and high effective population sizes. Moreover, we do not find distinct signals of strong selective sweeps anywhere in the genome, although we cannot rule out the possibility of highly polygenic evolution. Our observations suggest that phenotypic change in these populations is not constrained by irreversible loss of genomic variation and thus imply that former traits could be reestablished with demographic recovery.

As anthropogenic activities rapidly transform the environment, a fundamental question is whether wild populations have the capacity to adapt and evolve fast enough in response (1–3). Phenotypic change can result from phenotypic plasticity, but emerging examples of genomic change over only a few generations have made clear that rapid evolution is also possible (4–6). In the literature, one of the most dramatic and widely cited cases involves the declining age and size at maturation of Atlantic cod (Gadus morhua) following several generations of high fishing pressure (3, 7–10). Fisheries produce some of the fastest rates of phenotypic change ever observed in wild populations (2, 11), but the extent to which fisheries-induced evolution has occurred in the wild and the degree to which it is reversible remain strongly debated (12).The hypothesis that evolution underlies these phenotypic changes is supported by a range of observations. For example, theory on the selective nature of many fisheries reveals that higher rates of harvesting will—with only a few exceptions—favor earlier sexual maturation, greater investment in reproduction, and slower growth (13). In addition, experiments in the laboratory that selectively remove large or small individuals from a population reveal rapid evolution of body size and maturation time in only a few generations, as well as substantial impacts on fishery yields (14–16). Fisheries-induced evolution experiments in the laboratory also reveal selective sweeps through dramatic shifts in allele frequencies, loss of genetic diversity, and increases in linkage disequilibrium at specific locations in the genome (15, 17, 18).However, translating these findings to wild populations has been substantially more difficult. One concern is that phenotypic plasticity, gene flow, or spatial shifts in populations can also explain the substantial phenotypic and limited genotypic changes reported from the wild to date (10, 13, 19–23). The magnitude and rate of fisheries-induced evolution may also be quite small in the wild (19). While theory provides strong evidence that fishing can be a potent driver of evolutionary changes, a clear empirical demonstration of fisheries-induced evolution would require evidence that the observed change is genetic (13). Whether and to what extent the widespread genomic reorganization observed in experiments also occurs in wild-harvested populations therefore remain unknown.Genomic analyses of temporal samples before and after selective events have provided key opportunities to test for rapid adaptive evolution from standing genetic variation in wild populations by identifying unusually strong shifts in allele frequencies over time (4, 5). In addition, the history of genomic research with Atlantic cod (24, 25) provides a unique opportunity to test for genomic signatures of fisheries-induced evolution in particular. Archival samples collected by fisheries scientists decades or even centuries ago represent a valuable source of historical genomic material that can provide rare insight into the genetic patterns of the past (26). Here, we obtained whole-genome sequence data from well-preserved archives of Atlantic cod scales and otoliths (ear bones) that were originally collected from two populations on either side of the Atlantic Ocean: the northeast Arctic population sampled near Lofoten, Norway in 1907 and the Canadian northern cod population sampled near Twillingate, Newfoundland in 1940 (Fig. 1A and SI Appendix, Table S1). The Canadian northern population collapsed from overfishing in the early 1990s, while the northeast Arctic population experienced high fishing rates but smaller declines in biomass (10, 27, 28). Both populations have shown marked reductions in age at maturation, though with slight increases in maturation age in northeast Arctic cod after 2005 (Fig. 1B). We compared these historical genomes with modern data from the same locations (Fig. 1A and SI Appendix, Table S2). In total, we analyzed 113 individual genomes (Methods) from these two unique populations that had independently experienced intensive fishing during the last century (7, 10). We found a marked lack of large genomic changes or selective sweeps through time, suggesting instead that phenotypic plasticity or, potentially, highly polygenic evolution can explain the observed changes in phenotype.Open in a separate windowFig. 1.Spatiotemporal population structure based on genome-wide data in Atlantic cod. (A) In total, 113 modern and historical specimens were analyzed from northern cod collected in Newfoundland, Canada (1940, yellow; 2013, dark yellow) and from northeast Arctic cod collected in the Lofoten archipelago, Norway (1907, orange; modern: 2011, red; 2014, dark red). (B) Age at 50% maturity over time in each population. (C) PCA as implemented in PCAngsd. Velicier’s minimum average partial (MAP) test identified a single significant PC and only one PC is shown. Individuals are colored according to A. (D) Model-based ADMIXTURE ancestry components for historical (1907, 1940) and modern (2013, 2014) populations (k = 2; NGSadmix). Each individual is represented by a column colored to show the proportion of each ancestry component for Canada (dark yellow) and Norway (orange). Population differentiation based on pairwise weighted F_ST is also shown. (E) The correlation between the allele frequencies in historical and modern populations. Colors reflect the relative density of points, from darker (more density) to lighter (less density). R², coefficient of correlation.     

From the Cover: Invasiveness is linked to greater commercial success in the global pet trade

The pet trade has become a multibillion-dollar global business, with tens of millions of animals traded annually. Pets are sometimes released by their owners or escape, and can become introduced outside of their native range, threatening biodiversity, agriculture, and health. So far, a comprehensive analysis of invasive species traded as pets is lacking. Here, using a unique dataset of 7,522 traded vertebrate species, we show that invasive species are strongly overrepresented in trade across mammals, birds, reptiles, amphibians, and fish. However, it is unclear whether this occurs because, over time, pet species had more opportunities to become invasive, or because invasive species have a greater commercial success. To test this, we focused on the emergent pet trade in ants, which is too recent to be responsible for any invasions so far. Nevertheless, invasive ants were similarly overrepresented, demonstrating that the pet trade specifically favors invasive species. We show that ant species with the greatest commercial success tend to have larger spatial distributions and more generalist habitat requirements, both of which are also associated with invasiveness. Our findings call for an increased risk awareness regarding the international trade of wildlife species as pets.

The extraordinary movement of our own species through migration, colonization, and travel has driven the geographic expansion of countless other species since prehistoric times (1). Humans have deliberately introduced a diverse range of species, in particular domesticated crops and animals that have contributed to our success (1). Today, however, the trade in live organisms for nonutilitarian reasons has rocketed (2–4). In the last decade alone, billions of plants and animals comprising thousands of species were traded annually, fueling a multibillion-dollar global business (2, 3, 5, 6). In particular, the demand for nontraditional (also known as “exotic”) ornamentals and pets, i.e., organisms without a long history of domestication, has grown (2). These species are sometimes released into the wild or escape and may survive and reproduce (2, 7–9). Species with populations that have established outside of their native range are referred to as invasive species hereafter (see Table 1 for terminology). Some invasive species can have severe impacts on global biodiversity (10–13) and impose tremendous costs on society by damaging physical infrastructure, agriculture, forestry, and human health (14, 15). However, even though it is undisputed that the trade in pets and ornamentals contributes to the global movement of invasive animals (16–19), it is still unclear whether this trade specifically favors invasive species.Table 1.Glossary

Ecological adaptation in European eels is based on phenotypic plasticity

The relative role of genetic adaptation and phenotypic plasticity is of fundamental importance in evolutionary ecology [M. J. West-Eberhard, Proc. Natl. Acad. Sci. U.S.A. 102 (suppl. 1), 6543–6549 (2005)]. European eels have a complex life cycle, including transitions between life stages across ecological conditions in the Sargasso Sea, where spawning occurs, and those in brackish and freshwater bodies from northern Europe to northern Africa. Whether continental eel populations consist of locally adapted and genetically distinct populations or comprise a single panmictic population has received conflicting support. Here we use whole-genome sequencing and show that European eels belong to one panmictic population. A complete lack of geographical genetic differentiation is demonstrated. We postulate that this is possible because the most critical life stages—spawning and embryonic development—take place under near-identical conditions in the Sargasso Sea. We further show that within-generation selection, which has recently been proposed as a mechanism for genetic adaptation in eels, can only marginally change allele frequencies between cohorts of eels from different geographic regions. Our results strongly indicate plasticity as the predominant mechanism for how eels respond to diverse environmental conditions during postlarval stages, ultimately solving a long-standing question for a classically enigmatic species.

How species adapt to the diverse environmental conditions that they experience through their life is fundamental to understanding evolutionary processes. Species that occur across extreme environmental gradients must respond to a diverse range of conditions. This can be accomplished by individual-level phenotypic plasticity, meaning that individuals adjust their physiology to the prevailing environmental conditions (1), and by local genetic adaptation that may lead to reproductively isolated subpopulations.The European eel (Anguilla anguilla) provides a fascinating example. Adult eels are mostly found in freshwater bodies and brackish coastal areas from North Africa in the south to the North Cape in the north (well above the polar circle), from the Azores in the west to the Black Sea in the east. The spawning grounds long remained a mystery, but 100 y ago larvae were discovered in the Sargasso Sea, ∼7,000 km away from the mainland (2), and subsequent research has only recently begun to reveal the complexity of their spawning migration as maturing adults back to the Sargasso Sea (3–5). All European eels reproduce in the Sargasso Sea and the offspring drift passively, as leptocephali larvae, on oceanic currents toward the European continent. Consequently, which geographic region across Europe and North Africa eels inhabit as adults appears to be driven largely by a stochastic process. These observations are consistent with a single panmictic population and would preclude genetic adaptation to local conditions in Europe and North Africa. However, eel reproduction in the Sargasso Sea, which is 2,000 km wide, does not exclude the possibility of genetically differentiated subpopulations with distinct spawning areas or timing which may impact the likelihood of which geographic region the larvae reach after their trans-Atlantic migration. In fact, European eels and American eels (Anguilla rostrata), estimated to have split from a common ancestor around 3.75 million y ago (6), reproduce in parapatry in the Sargasso Sea but still maintain reproductive isolation with a low rate of hybridization (7, 8). Evaluating the hypothesis that the European eel consists of a single panmictic population is central to understanding eel ecology and evolution, and how eel populations may be affected by global change and other environmental threats (9, 10).Previous studies on the European eel using low-density marker sets (6, 11–15) or reduced representation sequencing (15) found little to no differentiation between geographic areas, consistent with a single panmictic population. However, low genetic differentiation at selectively neutral markers is a common observation in marine species with large geographical ranges and gene flow between subpopulations (16, 17), but does not necessarily capture patterns of local adaptation. For example, an early study based on 13 allozyme loci failed to identify genetic differentiation among Atlantic and Baltic herring (Clupea harengus) from diverse ecological conditions, and a single panmictic population could not be excluded (18). In sharp contrast, whole-genome sequencing revealed strong genetic differentiation between ecotypes of herring for a few percent of all genes (19–21). Thus, a lack of genetic differentiation even at many loci does not exclude the possibility that the European eel is divided into partially reproductively isolated subpopulations, genetically adapted to the diverse ecological conditions that individuals are exposed to during postlarval stages. Whole-genome data are required to evaluate the possibility that European eel populations are structured into subpopulations showing genetic differentiation.A high-quality reference genome for the European eel has recently been released by the Vertebrate Genomes Project (22). Here we used this assembly and low-coverage, individual whole-genome sequencing of 445 individuals from 10 geographic samples (median sequence coverage 1.4) covering most of the species range, stretching from Sweden to Ireland to Tunisia (Fig. 1A and SI Appendix, Table S1). We also include a sample of American eel (n = 49) as an outgroup. Our study scales considerably over existing genomics research on the European eel, both in number of individuals and number of loci sequenced (15). The sample design and sample size were selected to critically evaluate the presence of any genetic differentiation (due to drift or selection) between individuals from different geographic regions. We report a complete lack of genomic regions with significant differentiation between geographically separated samples.Open in a separate windowFig. 1.Sample overview and population structure analysis. (A) Sample locations (SI Appendix, Table S1) with NGSadmix (23) ancestry proportions for American (A. rostrata) and European (A. anguilla) eel samples (inset). K = 2 is best supported; see SI Appendix, Fig. S1 for additional analyses. (B) Principal-component analysis for whole-genome SNPs generated in PCAngsd (25), with points colored by sampling locality. PC 1 variance (39.9%) is driven by species divergence and A. rostrata ancestry in A. anguilla samples.     

Individual error correction drives responsive self-assembly of army ant scaffolds

An inherent strength of evolved collective systems is their ability to rapidly adapt to dynamic environmental conditions, offering resilience in the face of disruption. This is thought to arise when individual sensory inputs are filtered through local interactions, producing an adaptive response at the group level. To understand how simple rules encoded at the individual level can lead to the emergence of robust group-level (or distributed) control, we examined structures we call “scaffolds,” self-assembled by Eciton burchellii army ants on inclined surfaces that aid travel during foraging and migration. We conducted field experiments with wild E. burchellii colonies, manipulating the slope over which ants traversed, to examine the formation of scaffolds and their effects on foraging traffic. Our results show that scaffolds regularly form on inclined surfaces and that they reduce losses of foragers and prey, by reducing slipping and/or falling of ants, thus facilitating traffic flow. We describe the relative effects of environmental geometry and traffic on their growth and present a theoretical model to examine how the individual behaviors underlying scaffold formation drive group-level effects. Our model describes scaffold growth as a control response at the collective level that can emerge from individual error correction, requiring no complex communication among ants. We show that this model captures the dynamics observed in our experiments and is able to predict the growth—and final size—of scaffolds, and we show how the analytical solution allows for estimation of these dynamics.

Complex infrastructures of all kinds, whether technological, social, or biological, demand one thing above all else to function effectively: the property of resilience in the face of disruption. In this context, we consider a resilient system to be one “that returns to or exceeds its predisturbance level of performance following a perturbation,” as defined by Middleton and Latty in ref. 1. Understanding how resilience can be achieved in systems with many interacting components has emerged as a common goal across the disciplines of biology, engineering, and ecology (1, 2), and thus, it is important to identify and examine specific cases of infrastructural resilience in nature. Such examples abound among the social insects, in which difficult coordination problems are often solved through distributed control mechanisms relying on individual sensing and local interactions mediated by simple rules, without the need for complex communication (3, 4). Just as human societies rely on the organized flow of materials and information to function effectively, from the scale of the city to global trade networks (5–8), for many social insects, colony survival often depends on the effective coordination of foraging traffic to transport resources through the environment (9–11). These networks must be able to maintain their functionality when confronted with disruption and self-heal or otherwise rapidly respond to unpredictable conditions.To achieve resilience, biological systems at many scales, including these social insect infrastructures, utilize various mechanisms of feedback control. In a recent example, a core regulatory feedback mechanism common to many social insects was identified—the so-called common stomach—which achieves resilience through a saturation process that relies on negative feedback through integral control (12). Both positive and negative feedback have long been known to underlie self-organizing processes like the formation of ant foraging trails and collective decision making (13–15), and negative feedback has been shown to play a particularly important role among social insects in maintaining flexibility when environmental conditions change (16–18). Understanding how feedback control operates in distributed systems, like ant colonies, to facilitate resilience remains an important challenge, with implications across disciplines.The army ant Eciton burchellii, with large colonies of more than 500,000 workers and a distinctive nomadic life cycle (19), presents an ideal model for studying resilience and distributed control in natural systems. E. burchellii is considered a top predator in the neotropical forest leaf litter community (20) and a canonical example of self-organization among social insects (21). Unlike most ants that maintain a fixed nest location, colonies alternate every few weeks between a stationary and a nomadic phase, during which the entire colony moves nightly to a new temporary nest site, known as a bivouac (22, 23). Most workers spend the day foraging in massive swarm raids that start at dawn, combing the forest floor to flush out a diverse array of arthropod prey from the leaf litter (24, 25). At the front of the raid, fleeing prey items are attacked en masse, with those caught dismantled for transport along a branching trail network that can stretch for over 100 m back to the bivouac (26). Raids must be conducted at a rapid pace to maintain the high rates of prey delivery needed to support the ravenous appetite of the developing brood (27). This evolutionary pressure has led to a number of morphological and behavioral adaptations, including the spontaneous formation of traffic lanes that increase the efficiency of foraging trails (10), some of the highest running speeds among all ants (28), and a specialized porter caste morphology to facilitate the transport of bulky prey items, carried under the body (29).This suite of adaptations includes the remarkable capacity for self-assembly, by which ants can join together to create temporary structures that modify the environment (27, 30, 31). Various examples of self-assembly have been observed in other social insects (32, 33), and these structures have been described as intermediate-level parts of insect societies—adaptive units that function at a level between individual and colony (34). For E. burchellii and the closely related Eciton hamatum, self-assembled structures serve two primary functions: for shelter, in the case of the bivouac, and to facilitate the flow of traffic during raids and emigrations, with structures like plugs (27) and bridges (30, 31) that emerge along the trail network as needed. Our previous work on Eciton bridges has shown how these structures are robust to perturbations and responsive to traffic and environmental geometry, changing size and location to create shortcuts that benefit the colony.Here we reveal another type of self-assemblage in E. burchellii which we call “scaffolds,” given their function as temporary support structures. These structures have not been previously studied experimentally or clearly defined, although they have been described based on observations in previous literature without consistent terminology (35, 36). Scaffolds often form when the foraging trail crosses a sloped surface from which ants may slip and/or fall, such as a rock face, tree root, or even the wall of a building (as shown in SI Appendix, Movie S1). To form a scaffold, individual ants stop and grip the underlying surface with their tarsi, remaining stationary as traffic continues to pass. Scaffolds can be sparse, consisting of only a few dispersed ants, or extremely dense, with many ants overlapping one another in a continuous cluster (Fig. 1 and SI Appendix, Movie S2). To ascertain how scaffolds form and under what conditions, and to assess their effects on foraging traffic, we conducted field experiments with wild colonies of E. burchellii, using an apparatus that allowed for manipulation of the slope over which ants traversed (Fig. 1). We observed and quantified the growth of scaffold structures and measured traffic variables to assess their influence on scaffold formation.Open in a separate windowFig. 1.Overview of experimental apparatus and video analysis procedure. (A) Photo of apparatus setup in the field during an experiment, with platform set to 90°. (B) Diagram of 3D-printed apparatus assembly, with removable spacer to adjust platform angle shown in red. (C) Frames extracted from video of experiment with platform set to 90°, at experiment start (Left), after 120 s (Middle), and after 300 s (Right). (D) Diagram of image subtraction algorithm for video analysis. Moving ants are shown in gray, with stationary ants shown in black (indicating scaffold area over time) and the final scaffold area shown in red. (E) Dimensions of adjustable platform overlaid on photo of experiment in progress, with platform set to 80° (removable spacer to adjust angle visible at bottom left of platform).Our experimental results show that scaffolds repeatedly and predictably form on inclined surfaces, and we describe the relative influences of environmental and traffic variables on their growth. Informed by these observations, we propose a theoretical model that describes scaffold growth as a control response at the group level that emerges through a simple process of error correction at the individual level. Our model describes a hypothetical mechanism underlying scaffold formation that links individual-level sensing of disruption to the adaptive collective response of scaffold formation. Comparing the model predictions with our experimental observations, we show that the model captures well the observed dynamics of scaffold growth, and we derive an analytical solution that allows for straightforward prediction of these dynamics.     

Avian mud nest architecture by self-secreted saliva

Mud nests built by swallows (Hirundinidae) and phoebes (Sayornis) are stable granular piles attached to cliffs, walls, or ceilings. Although these birds have been observed to mix saliva with incohesive mud granules, how such biopolymer solutions provide the nest with sufficient strength to support the weight of the residents as well as its own remains elusive. Here, we elucidate the mechanism of strong granular cohesion by the viscoelastic paste of bird saliva through a combination of theoretical analysis and experimental measurements in both natural and artificial nests. Our mathematical model considering the mechanics of mud nest construction allows us to explain the biological observation that all mud-nesting bird species should be lightweight.

Bird nests come in a variety of forms made from diverse building materials (1, 2). Each type of bird nest is subjected to mechanical constraints imposed by material characteristics. To overcome these constraints, birds have devised brilliant architectural technologies, which provide inspiration for a novel materials processing scheme and help us to better understand animal behavior.For instance, some birds including storks (Cicioniidae) and eagles (Accipitidae) build nests by piling up hard filamentary materials such as twigs, harnessing their friction as the cohesion mechanism (3). Weaverbirds (Ploceidae) weave soft filamentary materials such as grass and fine leaves into a woven nest tied to a tree branch. Some bird species use their own saliva in nest building, which Darwin considered an example of natural selection (4). An extreme case is the Edible-nest Swiftlets, which build their nest purely of self-secreted saliva so that it can be attached to cliff walls and cave ceilings where the above twig piles and tied leaves are not allowed (5).Swallows (Hirundinidae), phoebes (Sayornis), and other mud nesters have developed a unique building material, a mixture of mud and their own saliva, in contrast to those made of purely collected or self-secreted materials (6) (Fig. 1). During construction, mud nesters repeatedly pile a beakful of wet mud on the nest, and liquid bridges are formed in the nest due to evaporation. While building a nest usually takes several weeks, a transition from wet to dry structures can occur within a few hours. Hence, the capillary forces of liquid bridges temporarily provide cohesion such as those in sandcastles. However, unlike sandcastles, dehydrated saliva comes into play for permanent cohesion after complete evaporation (SI Appendix, Supplementary Note 1).Open in a separate windowFig. 1.A nest of the barn swallow (H. rustica). (A) Photograph of a barn swallow nest, taken from under the ceiling of a house in Suwon-si, Gyunggi-do, South Korea (37°16′13.5″N 126°59′01.0″E). (B) SEM image of the nest surface. (C) Chemical composition analysis of the surface shown in B by EDS. The red area indicates a region containing mostly carbon atoms, which may originate from bird saliva. The green area indicates a region containing mostly the silicon atoms of clay particles.Mud itself cannot confer sufficient cohesion and adhesion in mud nests. The ability of mud nests to bear tensile loads originates from the gluing agent in the bird''s saliva, which permeates into granules as a liquid and binds them as a solid after solvent evaporation (6–8) (SI Appendix, Supplementary Note 2). The gluing agent is called mucin, a family of large glycoproteins that are ubiquitous in animal organs and form a mucus gel with versatile functionality (9). Fig. 1B shows the scanning electron microscopy (SEM) image of a barn swallow’s mud nest consisting of platelet clay particles and larger grains. Energy-dispersive spectroscopy (EDS) mapping image of Fig. 1C clearly shows regions corresponding to organic material which is presumed to be from bird’s saliva.Of particular interest and worth biophysical investigation are the tensile strength of the mud nest with hardened saliva, design principles associated with the saliva-originated strength, and the resulting effects on the evolution of these mud-nesting birds. Principles behind cohesion in granular materials, such as wet sands (10), cemented powder aggregates (11), construction materials (12), and pharmaceutical tablets (13), have been studied to date, exploring the stress transmission, elasticity, and failure (14–18), and the formation of solidified bridges (19–21). However, little attention has been paid to the cohesion effects of self-secreted polymer materials upon evaporation and the biologically constructed granular architecture like birds’ mud nests. Here we devised experimental techniques to measure the strength of the relatively small and fragile nest specimens in order to mechanically characterize birds’ mud nests. We elucidate how solutes from bird saliva generate solid bridges that give rise to macroscopic tensile strength, which has long awaited physicochemical explanation since its first observation (4). To characterize the design principle of bird''s mud nests, we investigated natural and three-dimensional (3D)-printed artificial nests with various tools for visualization and mechanical testing. Along with the experimental studies, we theoretically investigated the effects of biopolymer concentration on nest strength. This combination of theory and experiment suggests that there is a size limit for mud-nesting birds, which is supported by biological data.     

Three genomes in the algal genus Volvox reveal the fate of a haploid sex-determining region after a transition to homothallism

Transitions between separate sexes (dioecy) and other mating systems are common across eukaryotes. Here, we study a change in a haploid dioecious green algal species with male- and female-determining chromosomes (U and V). The genus Volvox is an oogamous (with large, immotile female gametes and small, motile male gametes) and includes both heterothallic species (with distinct male and female genotypes, associated with a mating-type system that prevents fusion of gametes of the same sex) and homothallic species (bisexual, with the ability to self-fertilize). We date the origin of an expanded sex-determining region (SDR) in Volvox to at least 75 Mya, suggesting that homothallism represents a breakdown of dioecy (heterothallism). We investigated the involvement of the SDR of the U and V chromosomes in this transition. Using de novo whole-genome sequences, we identified a heteromorphic SDR of ca 1 Mbp in male and female genotypes of the heterothallic species Volvox reticuliferus and a homologous region (SDLR) in the closely related homothallic species Volvox africanus, which retained several different hallmark features of an SDR. The V. africanus SDLR includes a large region resembling the female SDR of the presumptive heterothallic ancestor, whereas most genes from the male SDR are absent. However, we found a multicopy array of the male-determining gene, MID, in a different genomic location from the SDLR. Thus, in V. africanus, an ancestrally female genotype may have acquired MID and thereby gained male traits.

As first noted by Darwin when studying plant sexuality, self-fertilization can lead to inbreeding depression, though this potential disadvantage relative to outcrossing can be offset by higher probability of fertilization success (1, 2). Thus, transitions between inbreeding and outbreeding mating systems attract the attention of evolutionary biologists and have been documented in sexual systems across a broad range of taxa including animals, land plants, algae (Fig. 1), protists, and fungi (3–9). Hermaphroditic mating systems in diploid animals and plants are fairly common, and the evolution of hermaphroditism in animals has been extensively studied (10–12). The molecular genetic bases of hermaphroditism have also been recently studied in several flowering plants (13–15).Open in a separate windowFig. 1.Schematic representation of life cycles of two closely related species of Volvox sect. Merrillosphaera (SI Appendix, Fig. S1) showing different sexual systems (heterothallism and homothallism), based on Nozaki et al. (29). Note gonidia (g) in asexual spheroids and sperm packets (sp) and eggs (e) in sexual spheroids. (A) NIES-3781. (B) NIES-3783. (C and D) NIES-3782. (E and H) NIES-3784. (F and G) NIES-3780. For materials and methods of light micrographs, refer to Nozaki et al. (29). All photographs are original.In haploid organisms, two basic types of mating systems are recognized: heterothallic (with two or more self-incompatible mating types in isogamous species or with males and females in anisogamous/oogamous species) and homothallic (self-compatible with isogamy or bisexual with anisogamy/oogamy) (SI Appendix, Table S1). In heterothallic species, gamete compatibility in isogamy or maleness versus femaleness is usually determined by a single complex locus on a chromosome (9). The process by which this genetic system breaks down to allow a single genotype to acquire functions of both sexes and mating-types in the homothallic species is poorly understood outside of fungi (8).Organisms with a heterothallic haploid generation such as algae and early diverging land plants like mosses and liverworts have sex-determining regions (SDRs) that are differentiated between the two sexes and located on male and female determining (U and V) chromosomes (9). SDRs exhibit suppressed recombination and harbor fully sex-linked genes including sex-specific genes (found in only one of the two sexes) and gametologs (genes with alleles in both haplotypes of the SDR) (9). A recent study has characterized genes involved in a change from sexual reproduction to parthenogenesis in a heterothallic brown alga with UV chromosome (16). However, genomics of transitions between heterothallism and homothallism have not been previously studied in the context of haploid SDRs, and the fates of sex-determining and sex-related genes present on ancestral SDRs in UV chromosomes after transitions to homothallism are unknown.The volvocine green algal lineage is an especially well-studied evolutionary model for investigating the origins of sexes and transitions in sexuality, because it includes extant organisms with a graded range of sexual or mating phenotypes from unicellular isogamous Chlamydomonas through genera with increasing degrees of sexual dimorphism, such as multicellular, isogamous Gonium, and oogamous Volvox (Fig. 2) (17–22). In heterothallic species of the isogamous volvocine genera such as Chlamydomonas, there are two mating-types, plus and minus, both of which are needed for sexual reproduction (9). Ferris and Goodenough (23) characterized the first mating-type locus (MT) in Chlamydomonas reinhardtii, which was discovered to have features similar to those of a typical SDR, including large size and suppressed recombination. Subsequent studies of heterothallic species of multicellular volvocine algae revealed conservation of orthologs of the mating-type determining gene MID in the mating type minus genotype of the isogamous genera Gonium and Yamagishiella and in the SDR of males of anisogamous Eudorina and the oogamous Volvox (19, 21, 24). This indicates that the evolution of anisogamy has involved genetic changes that are closely associated with the MT, presumably because this ensures the ability of male and female gametes to fuse with each other as proposed on theoretical grounds (25). The MT or SDR haplotypes in heterothallic volvocine algae range from 7 kbp to 1 Mbp and are structurally heteromorphic with chromosome rearrangements that distinguish the haplotypes and various degrees of genetic differentiation between them (21) (Fig. 2). The ca 1 Mbp, highly differentiated Volvox carteri SDR has expanded in size about fivefold relative to that in Chlamydomonas (19), and this was likely to have occurred after the transition from isogamy to anisogamy/oogamy (21). However, the timing of the SDR expansion in the genus Volvox and its significance with respect to the evolution of sexual dimorphism remain unknown.Open in a separate windowFig. 2.Volvocine green algal phylogeny and SDR or MT evolution. Asexual or vegetative phase, sexual phase, MTM and MTF, or SDLR and phylogenetic positions of V. reticuliferus and V. africanus are illustrated with those of five other species previously studied [C. reinhardtii, Gonium pectorale, Yamagishiella unicocca, Eudorina sp., and V. carteri (19, 21, 23, 24)]. Note that all three of these Volvox species belongs to the section Merrillosphaera (SI Appendix, Fig. S1).To date, genomic studies of volvocine green algae have focused on heterothallic species (Fig. 2), but most genera within this lineage have members that also underwent transitions from heterothallism to homothallism (7). These algae are found in freshwater habitats and form a dormant diploid zygote with a resistant wall that allows it to survive under unfavorable conditions; diploidy in the zygote could be of selective value for survival when chromosomal damage and mutation may be occurring at increased frequencies (26). Thus, transitions from heterothallism to homothallism in the volvocine algae may be favored for production of the resistant diploid cells (zygotes) by self-fertilization. We previously noted the presence of a MID gene in two homothallic species of the genus Volvox (27), and an artificial homothallic phenotype was shown in the male genotype of V. carteri (28). However, it remains unknown how a naturally occurring homothallic mating system could arise from an ancestral SDR in UV chromosome system (or vice versa).In most members of the genus Volvox, the transition from vegetative to sexual reproduction involves modified embryogenesis and the production of morphologically distinct male and/or female sexual spheroids (7). Recently, we characterized two species of Volvox that are closely related but have different mating systems: heterothallic Volvox reticuliferus, which upon sexual induction makes differentiated male or female, and homothallic Volvox africanus, which produces both male spheroids and bisexual spheroids from single clonal cultures (27, 29) (Figs. 1 and ​and22 and SI Appendix, Table S1). These two species diverged ca 11 MYA; together with the more distantly related V. carteri, they represent a monophyletic group, the “section Merrillosphaera,” which originated ca 75 MYA and is an infrageneric taxon of the polyphyletic genus Volvox (SI Appendix, Fig. S1). While ancestral heterothallism in Merrillosphaera is more likely than homothallism, statistical support for this inference is not strong, and the directionality of transitions between these two types of sexuality within this clade remains somewhat inconclusive (7).Here, we performed de novo whole-genome sequencing of male and female genotypes of heterothallic V. reticuliferus and of homothallic V. africanus in order to characterize their SDR and SD-like regions (SDLRs), respectively. We used this sequence information to infer the likely ancestral state of heterothallism in this clade and to reconstruct the SDR changes which occurred during the evolution of homothallism in V. africanus, including tracking the fates of male- and female-specific genes and male versus female gametologs derived from the putative heterothallic ancestral species.     

Capillary equilibrium of bubbles in porous media

In geologic, biologic, and engineering porous media, bubbles (or droplets, ganglia) emerge in the aftermath of flow, phase change, or chemical reactions, where capillary equilibrium of bubbles significantly impacts the hydraulic, transport, and reactive processes. There has previously been great progress in general understanding of capillarity in porous media, but specific investigation into bubbles is lacking. Here, we propose a conceptual model of a bubble’s capillary equilibrium associated with free energy inside a porous medium. We quantify the multistability and hysteretic behaviors of a bubble induced by multiple state variables and study the impacts of pore geometry and wettability. Surprisingly, our model provides a compact explanation of counterintuitive observations that bubble populations within porous media can be thermodynamically stable despite their large specific area by analyzing the relationship between free energy and bubble volume. This work provides a perspective for understanding dispersed fluids in porous media that is relevant to CO₂ sequestration, petroleum recovery, and fuel cells, among other applications.

Bubbles are generated, trapped, and mobilized within porous media as a consequence of incomplete fluid–fluid displacements (1, 2), phase changes (3, 4), chemical and biochemical reactions (5, 6), or injection of emulsified fluids and foams (7, 8). Compared to continuously connected phases, the behavior of dispersed bubbles, or ganglia, are far less understood. In particular, the thermodynamic stability of bubbles, despite their large specific surface area, remains a puzzle. The difficulty comes from the fact that each bubble can attain a volume (V), topology, and capillary pressure (P_c) that is distinct from other bubbles in the medium (9). The variability poses challenges to understanding the transport and trapping mechanisms of bubbles in geologic CO₂ sequestration (10, 11), hydrocarbon recovery (12, 13), fuel cell water management (14, 15), and vadose zone oxygen supply (16, 17).The dominant factor controlling a bubble’s behavior in a porous medium is capillarity, which is typically much larger than either viscous, gravitational, or inertial forces (18, 19). Capillary pressure, P_c, allows a closure relationship for two-phase Darcy Eqs. (20–22) and influences thermodynamic properties like phase partition (23). Capillary pressure is derived from the Young–Laplace equation P_c = γκ, where γ is the interfacial tension and κ is the surface curvature. In an open space without obstacles, a bubble spontaneously evolves into a sphere to minimize its total interfacial energy. Thus, P_c is a continuous and monotonically decreasing function of V (Fig. 1A). However, in a porous medium, bubble’s P_c–V relation is more complicated due to the geometric confinement imposed by the porous structure and topological evolution (24). A bubble can no longer remain spherical as it grows in size but must conform to the geometry of the pore(s) it occupies. Therefore, a bubble’s P_c is a function of not only its volume and interfacial tension but also its topology as dictated by the confining porous medium, as confirmed by recent laboratory experiments and numerical simulations (25–29). The mere presence of confinement therefore engenders a host of phenomena that would otherwise be absent, such as capillary trapping (30, 31), anticoarsening of bubble populations (32, 33), and complex ganglion dynamics (11, 18). Furthermore, theoretical studies in mathematical topology (28, 34, 35) prove that immiscible fluids can be fully characterized by d+1 Minkowski functionals, where d is the problem dimension. Such characterizations remove the path-dependent (or hysteretic) behavior common to these systems (34, 35).Open in a separate windowFig. 1.(A) Spherical bubbles inside a bulk fluid. (B) Micromodel observations show that bubbles are nonspherical in porous media and may occupy multiple pores. This image is from SI Appendix, Movie S1. (C) A 2D porous medium comprised of an ordered array of identical circular grains. A bubble occupying multiple pores including a zoom-in to a portion of it. (D) Illustration of the full state. (E) Illustration of the critical state. (F) Decomposition of a bubble into four distinct parts: minor arc menisci shown by dark blue cap-shaped regions, throats shown by light blue diamond-shaped regions, inner bulk bodies shown by red star-shaped regions, and major arc menisci shown by dark green cap-shaped regions.Recent developments in microfluidics and micro computed tomography imaging allow detailed pore-scale visualizations of fluids inside porous media, including the morphology of bubbles and ganglia (25, 36–39). Garing et al. (25) experimentally measured the equilibrium capillary pressure of trapped air bubbles inside sandstone and bead-pack samples. They found that, unlike bubbles within a bulk fluid, the P_c of trapped bubbles shows no clear dependence on V and seems to fall within a bounded interval, except for vanishingly small V. Xu et al. (40) proposed an empirical correlation for the P_c trapped bubbles based on microfluidic observations. In this correlation, as V increases, P_c decreases until a minimum is reached and then increases linearly. In the first stage, the bubble is unconfined, whereas in the second, it is reshaped by the surrounding solid walls. The proposed correlation, however, is only valid for bubbles in a single pore and not bubbles that span multiple pores. The latter seems to be rather common in nature as evidenced by recent direct observations (Fig. 1B) (2, 25).Here, we propose a simple conceptual model to describe the equilibrium states of a bubble with arbitrary size trapped inside a porous medium. The model accounts for the bubble’s morphology, the geometry of the solid matrix, and the wettability between the two. We derive all metastable configurations of the bubble analytically and highlight the thermodynamic states the bubble assumes when it is static, growing, or shrinking. We also show that the relationship between surface free energy (F) and volume (V) of large bubbles is approximately linear, which explains the previously counterintuitive observation that such bubbles are thermodynamically stable despite having large surface areas. Our work provides a step toward understanding the capillary state, stability, and evolution of dispersed immiscible fluids in porous media.