Archean komatiite volcanism controlled by the evolution of early continents

The generation and evolution of Earth’s continental crust has played a fundamental role in the development of the planet. Its formation modified the composition of the mantle, contributed to the establishment of the atmosphere, and led to the creation of ecological niches important for early life. Here we show that in the Archean, the formation and stabilization of continents also controlled the location, geochemistry, and volcanology of the hottest preserved lavas on Earth: komatiites. These magmas typically represent 50–30% partial melting of the mantle and subsequently record important information on the thermal and chemical evolution of the Archean–Proterozoic Earth. As a result, it is vital to constrain and understand the processes that govern their localization and emplacement. Here, we combined Lu-Hf isotopes and U-Pb geochronology to map the four-dimensional evolution of the Yilgarn Craton, Western Australia, and reveal the progressive development of an Archean microcontinent. Our results show that in the early Earth, relatively small crustal blocks, analogous to modern microplates, progressively amalgamated to form larger continental masses, and eventually the first cratons. This cratonization process drove the hottest and most voluminous komatiite eruptions to the edge of established continental blocks. The dynamic evolution of the early continents thus directly influenced the addition of deep mantle material to the Archean crust, oceans, and atmosphere, while also providing a fundamental control on the distribution of major magmatic ore deposits.Volcanism on Earth is the dynamic surface expression of our planet’s thermal cycle, with heat created from radioactive decay and lost through mantle convection (1). In the Archean eon (>2.5 bya), Earth’s heat flux was significantly higher than that observed today (1, 2) due to the combined effects of a more radioactive mantle (1, 3) and residual heat from planetary accretion (4). This resulted in the eruption of komatiites: ultra-high temperature, low-viscosity lavas with MgO >18% and eruption temperatures >1,600 °C (5) formed from mantle plumes (1, 2, 6). These rare, ancient magmas are dominantly restricted to the early history of the planet (3.5–1.5 Ga; ref. 7) and represent the remnants of huge volcanic flow fields (8) consisting of the hottest lavas preserved on Earth (5, 9, 10). These now-extinct volcanic systems and flow complexes had the potential to cover significant portions of the early continents, and were likely analogous to large igneous provinces in size and magma volume (11, 12). Komatiites are vital to our understanding of Earth’s thermal evolution (1–3, 7, 13–16), and represent a window into the dynamic secular development of the mantle throughout the early history of our planet (5). Subsequently, understanding the physical and chemical processes that govern their localization, volcanology, and geochemistry is vital in deciphering this information.In the Yilgarn Craton of Western Australia (Fig. 1), two major pulses of komatiite activity occurred at ∼2.9 Ga (southern Youanmi Terrane; refs. 17–19) and 2.7 Ga (Kalgoorlie Terrane, Eastern Goldfields Superterrane; refs. 5, 10). These represent two separate plume events that impinged onto preexisting continental crust (20–23), with the resulting magmas heterogeneously distributed across the craton (8, 10, 17–19, 23, 24). In this study, we provide the first evidence of a fundamental relationship between the spatiotemporal variation in komatiite abundance, geochemistry, and volcanology and the evolution of an Archean microcontinent, reflected in the changing isotopic composition of the crust.Open in a separate windowFig. 1.Map of the Archean Yilgarn Craton showing the basic granite-greenstone bedrock geology and location of the ∼2.9 and 2.7 Ga komatiite localities. Individual terranes/domains (39, 40) are labeled. Greenstone belts are labeled as follows: MD, Marda–Diemals; SC, Southern Cross; FO, Forrestania; LJ, Lake Johnston; RAV, Ravensthorpe; AW, Agnew–Wiluna; and KAL, Kalgoorlie/Kambalda. Komatiite localities are from Barnes and Fiorentini (10) (Table S4).We used Lu-Hf and U-Pb isotopic techniques on multiple magmatic and inherited zircon populations from granitoid rocks and felsic volcanic units, which represent the exposed Archean crust of the Yilgarn Craton. All zircon grains were dated using the sensitive high-resolution ion microprobe (SHRIMP), before in situ laser ablation inductively coupled plasma mass-spectrometry (LA-ICP-MS) analysis for Lu-Hf isotopes. The Lu-Hf isotopic data are expressed as εHf, which denotes the derivation of the ¹⁷⁶Hf/¹⁷⁷Hf ratio of the sample from the contemporaneous ratio of the chondritic uniform reservoir (CHUR), multiplied by 10⁴. The term “juvenile” refers to crustal material that plots on or close to the depleted mantle evolution line, suggesting derivation from a depleted mantle source. In contrast, “reworked” refers to the remobilization of preexisting crust by partial melting and/or erosion and sedimentation (25, 26). Complete sample information, methodology, and geochemical datasets (U-Pb, Lu-Hf, and komatiite) are available in the Supporting Information, Figs. S1–S3, and Tables S1–S4.The Lu-Hf data are displayed as a series of time-slice contour maps, which show “snapshots” of the changing source and age of the crust at 3,050–2,820; 2,820–2,720; and 2,720–2,600 Ma [Figs. 2–4; intervals based on the work of Mole et al. (21); full Hf dataset displayed in Fig. S4]. In these maps, point data representing the median εHf value of granites and felsic volcanics are plotted as contour maps that show the spatial extent of “blocks” of specific Lu-Hf isotopic character and their evolution through time. This method is based on previous isotopic mapping of the Yilgarn Craton using the analogous Sm-Nd system (27). Importantly, the Lu-Hf data presented here replicate the features of the Sm-Nd work (27, 28), with the added ability to look further back in time due to the in situ analysis of abundant inherited zircons (21).Open in a separate windowFig. 2.Lu-Hf (εHf) map of the Yilgarn Craton at 3,050–2,820 Ma. (A) Hf isotope map with the location of sample sites and komatiite localities. (B) Interpretive map of the area, showing the individual crustal blocks identified from the Hf isotope map and corresponding probability density plots. The blue curve represents the median εHf for discrete temporal groups (n_g), whereas the red curve represents all of the individual grain analyses (n_a). Dark gray polygons shown in the background of all maps represent supracrustal belts (Fig. 1).Open in a separate windowFig. 4.Lu-Hf (εHf) map of the Yilgarn Craton at 2,720–2,600 Ma. (A) Hf isotope map with the location of sample sites and komatiite localities. (B) Interpretive map of the area, showing the individual crustal blocks identified from the Hf isotope map and corresponding probability density plots. The blue curve represents the median εHf for discrete temporal groups (n_g), whereas the red curve represents all of the individual grain analyses (n_a). Dark gray polygons shown in the background of all maps represent supracrustal belts (Fig. 1).The variable isotopic signatures of the crust (Figs. 2–4) can be interpreted as proxies for lithospheric thickness (Figs. 5 and ​and6;6; ref. 29), where young, juvenile εHf values (εHf > 0) indicate relatively thin lithosphere and old, reworked values (εHf < 0) reflect thicker lithosphere (29, 30); a pattern observed in the modern-day western United States (29–32). Here, this information is combined to document the four-dimensional lithospheric architecture of the Yilgarn Craton and development of an Archean microcontinent.Open in a separate windowFig. 5.Isotopic cross-section and interpreted lithospheric architecture during the emplacement of ∼2.9 Ga komatiites in the southern Youanmi Terrane. (A) εHf map showing the isotopic architecture at the time of the ∼2.9 Ga plume emplacement, with the approximate extent of the plume head (red) and tail (yellow) shown for scale; (B) isotopic cross-section (A–A′) documenting the changing εHf of the crust from east to west (circles represent median; squares represent individual analyses) together with the occurrence and MgO content (Table S4) of ultramafic–mafic magmatism; and (C) interpreted lithospheric architecture based on the changing isotopic properties of the crust. The white ellipses represent the types of magma available in a particular area and the dashed lines show the approximate limits of their source regions. TiO₂ vs. Al₂O₃ data (17, 19) are shown for komatiites of the relevant greenstone belts, demonstrating the progressive eastward homogenization of Barberton-type melts. The location of the continent core and continent edge are shown based on the Lu-Hf data. The eruption of komatiite would likely have been facilitated by plume-related extension at this interface. Approximate thickness values for developed Archean lithosphere (∼250–150 km) were taken from Boyd et al. (41), Artemieva and Mooney (42), and Begg et al. (43). The approximate scale of the plume head (∼1,600 km), tail (200–100 km), and thickness (∼150–100 km) were taken from Campbell et al. (15). Note that the plume-tail material moves above the plume head, despite impacting the lithosphere later, as it is hotter, more buoyant, and subsequently emplaced at higher flux (15).Open in a separate windowFig. 6.Isotopic cross-sections and interpreted lithospheric architecture during the emplacement of the ∼2.7 Ga komatiites in the Eastern Goldfields (Kalgoorlie Terrane). (A) εHf map showing the isotopic architecture at the time of the ∼2.7 Ga plume emplacement, with the approximate extent of the plume head (red) and tail (yellow) shown for scale; (B) isotopic cross-section (B–B′) documenting the changing εHf of the crust from east to west (circles represent median; squares represent individual analyses) together with the occurrence and MgO content (Table S4) of ultramafic–mafic magmatism; and (C) interpreted lithospheric architecture based on the changing isotopic properties of the crust. The white ellipses represent the types of magma available in a particular area and the dashed lines show the approximate limits of their source regions. TiO₂ vs. Al₂O₃ data (9) are shown for komatiites of the Eastern Goldfields, demonstrating the eastward dilution and removal of Barberton-type melts. The location of the continent core and continent edge are shown based on the Lu-Hf dataset. The eruption of komatiite would likely have been facilitated by plume-related extension at this interface. The dashed red lines shown in C account for the potential variation in lithospheric architecture within the Lake Johnston block based on the Lu-Hf data. External data used to construct these diagrams are the same as for Fig. 5. Note that the plume-tail material moves above the plume head, despite impacting the lithosphere later, as it is hotter, more buoyant, and subsequently emplaced at higher flux (15).     

Dietary options and behavior suggested by plant biomarker evidence in an early human habitat

The availability of plants and freshwater shapes the diets and social behavior of chimpanzees, our closest living relative. However, limited evidence about the spatial relationships shared between ancestral human (hominin) remains, edible resources, refuge, and freshwater leaves the influence of local resources on our species’ evolution open to debate. Exceptionally well-preserved organic geochemical fossils—biomarkers—preserved in a soil horizon resolve different plant communities at meter scales across a contiguous 25,000 m² archaeological land surface at Olduvai Gorge from about 2 Ma. Biomarkers reveal hominins had access to aquatic plants and protective woods in a patchwork landscape, which included a spring-fed wetland near a woodland that both were surrounded by open grassland. Numerous cut-marked animal bones are located within the wooded area, and within meters of wetland vegetation delineated by biomarkers for ferns and sedges. Taken together, plant biomarkers, clustered bone debris, and hominin remains define a clear spatial pattern that places animal butchery amid the refuge of an isolated forest patch and near freshwater with diverse edible resources.Spatial patterns in archaeological remains provide a glimpse into the lives of our ancestors (1–5). Although many early hominin environments are interpreted as grassy or open woodlands (6–8), fossil bones and plant remains are rarely preserved together in the same settings. As a result, associated landscape reconstructions commonly lack coexisting fossil evidence for hominins and local-scale habitat (microhabitat) that defined the distribution of plant foods, refuge, and water (7). This problem is exacerbated by the discontinuous nature and low time resolution often available across ancient soil (paleosol) horizons, including hominin archaeological localities. One notable exception is well-time-correlated 1.8-million-y-old paleosol horizons exposed at Olduvai Gorge. Associated horizons contain exceptionally preserved plant biomarkers along with many artifacts and fossilized bones. Plant biomarkers, which previously revealed temporal patterns in vegetation and water (8), are well preserved in the paleosol horizon and document plant-type spatial distributions that provide an ecosystem context (9, 10) for resources that likely affected the diets and behavior of hominin inhabitants.Plant biomarkers are delivered by litter to soils and can distinguish plant functional type differences in standing biomass over scales of 1–1,000 m² (11). Trees, grasses, and other terrestrial plants produce leaf waxes that include long-chain n-alkanes such as hentriacontane (nC₃₁), whereas aquatic plants and phytoplankton produce midchain homologs (e.g., nC₂₃) (12, 13). The ratio of shorter- versus long-chain n-alkane abundances distinguish relative organic matter inputs from aquatic versus terrestrial plants to sediments (13):P_aq = (nC₂₃ + nC₂₅)/(nC₂₃ + nC₂₅ + nC₂₉ + nC₃₁).Sedges and ferns are prolific in many tropical ecosystems (14). These plants both have variable and therefore nondiagnostic n-alkane profiles. However, sedges produce distinctive phenolic compounds [e.g., 5-n-tricosylresorcinol (ⁿR₂₃)] and ferns produce distinctive midchain diols [e.g., 1,13-dotriacontanediol (C₃₂-diol)] (SI Discussion).Lignin monomers provide evidence for woody and nonwoody plants. This refractory biopolymer occurs in both leaves and wood, serves as a structural tissue, and accounts for up to half of the total organic carbon in modern vegetation (11). Lignin is composed of three phenolic monomer types that show distinctive distributions in woody and herbaceous plant tissues. Woody tissues from dicotyledonous trees and shrubs contain syringyl (S) and vanillyl (V) phenols (12), whereas cinnamyl (C) phenols are exclusively found in herbaceous tissues (12). The relative abundance of C versus V phenols (C/V) is widely used to distinguish between woody and herbaceous inputs to sedimentary and soil organic matter (15).Plant biomarker ¹³C/¹²C ratios (expressed as δ¹³C values) are sensitive indicators of community composition, ecosystem structure, and climate conditions (8). Most woody plants and forbs in eastern Africa use C₃ photosynthesis (6), whereas arid-adapted grasses use C₄ photosynthesis (8, 14). These two pathways discriminate differently against ¹³C during photosynthesis, resulting in characteristic δ¹³C values for leaf waxes derived from C₃ (about –36.0‰) and C₄ (–21.0‰) plants (16). Carbon isotopic abundances of phenolic monomers of lignin amplify the C₃–C₄ difference and range between ca. –34.0‰ (C₃) and –14.0‰ (C₄) in tropical ecosystems (15). Terrestrial C₃ plant δ¹³C values decrease with increased exposure to water, respired CO₂, and shade (8), with lowest values observed in moist regions with dense canopy (17). Although concentration and δ¹³C values of atmospheric CO₂ can affect C₃ plant δ¹³C values (17), this influence is not relevant to our work here, which focuses on a single time window (SI Discussion). The large differences in leaf-wax δ¹³C values between closed C₃ forest to open C₄ grassland are consistent with soil organic carbon isotope gradients across canopy-shaded ground surfaces (6) and serve as a quantitative proxy for woody cover (f_woody) in savannas (8).As is observed for nonhuman primates, hominin dietary choices were likely shaped by ecosystem characteristics over habitat scales of 1–1,000 m² (3–5). To evaluate plant distributions at this small spatial scale (9), we excavated 71 paleosol samples from close-correlated trenches across a ∼25,000-m² area that included FLK Zinjanthropus Level 22 (FLK Zinj) at Olduvai Gorge (Fig. 1). Recent excavations (18–21) at multiple trenches at four sites (FLK^NN, FLK^N, FLK, and FLK^S, Fig. 1D) exposed a traceable thin (5–50 cm), waxy green to olive-brown clay horizon developed by pedogenic alterations of playa lake margin alluvium (22). Weak stratification and irregular redox stains suggest initial soil development occurred during playa lake regression (18, 22), around 1.848 Ma (ref. 23 and SI Discussion). To date, craniodental remains from at least three hominin individuals (18–20), including preadolescent early Homo and Paranthropus boisei, were recovered from FLK Zinj. Fossils and artifacts embedded in the paleosol horizon often protrude into an overlying airfall tuff (18, 19), which suggests fossil remains were catastrophically buried in situ under volcanic ash. Rapid burial likely fostered the exceptional preservation of both macrofossils (10) and plant biomarkers across the FLK Zinj land surface.Open in a separate windowFig. 1.Location and map of FLK Zinj paleosol excavations. (A and B) Location of FLK Zinj as referenced to reconstructed depositional environments at Olduvai Gorge during the early Pleistocene (18, 22) and the modern gorge walls. The perennial lake contained shallow saline–alkaline waters that frequently flooded the surrounding playa margin (i.e., floodplain) flats. (C) Outline of FLK Zinj paleosol excavation sites used for our spatial biomarker reconstructions. (D) Concentric (5 m) gridded distribution map of FLK Zinj paleosol excavations relative to previous archaeological trenches (18–21). Major aggregate complexes (FLK^NN, FLK^N, FLK, and FLK^S) are color-coded to show excavation-site associations.Plant biomarker signatures reveal distinct types of vegetation juxtaposed across the FLK Zinj land surface (Figs. 2–4 and Fig. S1). In the northwest, FLK^NN trenches show high nC₂₃ δ¹³C values (Fig. 2B) as well as high C/V and P_aq values (Figs. 3 and ​and4A).4A). They indicate floating or submerged aquatic plants (macrophytes) in standing freshwater (13), a finding that is consistent with nearby low-temperature freshwater carbonates (tufa), interpreted to be deposited from spring waters (22). Adjacent FLK^N trenches have lower P_aq values (Fig. 4A) with occurrences of fern-derived C₃₂-diol and sedge-derived ⁿR₂₃ (Fig. 2 C and D). These biomarker distributions indicate an abrupt (around 10 m) transition from aquatic to wetland vegetation. Less than 100 m away (Fig. 1C), low nC₃₁ δ¹³C values (Fig. 2A) and low C/V and very low P_aq values (Figs. 3 and ​and4A)4A) collectively indicate dense woody cover (Fig. 4B). In the farthest southeastern (FLK^S) trenches, high C/V values and high δ¹³C values for C lignin phenols (Fig. 3) indicate open C₄ grassland.Open in a separate windowFig. 2.Spatial distributions and δ¹³C values for plant biomarkers across FLK Zinj. Measured and modeled δ¹³C values (large and smaller circles, respectively) are shown for (A) nC₃₁ from terrestrial plants, (B) nC₂₃ from (semi)aquatic plants, (C) C₃₂-diol from ferns, and (D) ⁿR₂₃ from sedges (see refs. 12 and 13 and SI Discussion). Modeled values [inverse distance-weighted (9)] account for spatial autocorrelation (15-m radius) in standing biomass (35) over scales of soil organic matter accumulation (11). Black dots represent paleosols with insufficient plant biomarker concentrations for isotopic analysis.Open in a separate windowFig. 3.Molecular and isotopic signatures for lignin phenols across FLK Zinj. Bivariate plots are shown for diagnostic lignin compositional parameters (see refs. 12 and 15 and Fig. 1C). Symbols are colored according to respective δ¹³C values for the C lignin phenol, p-coumaric acid. FLK symbols are uncolored due to insufficient p-coumaric acid concentrations for isotopic analysis. Representative lignin compositional parameters (12, 15) are shown for monocotyledonous herbaceous tissues (G), dicotyledonous herbaceous tissues (H), cryptogams (N), and dicotyledonous woody tissues (W).Open in a separate windowFig. 4.Spatial relationships shared between local plant resources and hominin remains. Measured and modeled values (large and smaller circles, respectively) are shown for (A) P_aq (13) and (B) f_woody (8). Modeled values [inverse distance-weighted (9)] account for spatial autocorrelation (15-m radius) in standing biomass (35) over scales of soil organic matter accumulation (11). (C) Kernel density map of cut-marked bones (18–21) across the FLK Zinj land surface (Fig. S4). High estimator values indicate hotspots of hominin butchery (Fig. S5). A shaded rectangle captures the area (ca. 0.68 probability mass) with highest cut-marked bone densities and is shown in A and B for reference.Open in a separate windowFig. S1.Total ion chromatograms for saturated hydrocarbons in representative paleosols at (A) FLK^NN, (B) FLK^N, (C) FLK, and (D) FLK^S. C₂₃, tricosane; C₂₅, pentacosane; C₂₉ nonacosane; C₃₁, hentriacontane.Biomarkers define a heterogeneous landscape at Olduvai and suggest an influence of local resources on hominin diets and behavior. It is recognized (2, 24–26) that early Homo species and P. boisei had similar physiological characteristics. These similarities in physical attributes suggest behavioral differences were what allowed for overlapping ranges and local coexistence (sympatry) of both hominins. For instance, differences in seasonal subsistence strategies or different behavior during periods of drought and limited food could have reduced local hominin competition and fostered diversification via niche specialization (27–29).Physical and isotopic properties of fossil teeth indicate P. boisei was more water-dependent [low enamel δ¹⁸O values (24)] and consumed larger quantities of abrasive, ¹³C-enriched foodstuffs [flat-worn surfaces (25) and high enamel δ¹³C values (26)] than coexisting early Homo species. Although ¹³C-enriched enamel is commonly attributed to consumption of C₄ grasses or meat from grazers (14), this was not likely, because P. boisei craniodental features are inconsistent with contemporary gramnivores (24, 25) or extensive uncooked flesh mastication (26). Numerous scholars have proposed the nutritious underground storage organs (USOs) of C₄ sedges were a staple of hominin diets (14, 24, 26, 27). Consistent with this suggestion, occurrences of ⁿR₂₃ attest to the presence of sedges at FLK^NN and FLK^N (Fig. 2D). However, the low δ¹³C values measured for ⁿR₂₃ at these same sites (Fig. 2D and Fig. S2) indicate C₃ photosynthesis (12, 16), a trait common in modern sedges that grow in alkaline wetlands and lakes (30) (Fig. S3). Thus, biomarker signatures support the presence of C₃ sedges in the wetland area of FLK Zinj.Open in a separate windowFig. S2.Total ion chromatogram [TIC (A)] and selected ion chromatograms for derivatized 5-n-alkylresorcinols [m/z 268 (●)] and midchain diols [m/z 369 (○)] from a representative paleosol at FLK^N. Also shown are δ¹³C values for homologous (B) 5-n-alkylresorcinols and (C) midchain diols. C₃₂-diol, dotriacontanediol; ⁿR₂₃, tricosylresorcinol.Open in a separate windowFig. S3.Summary phyogenetic consensus tree of Cyperaceae (sedges) based on nucleotide (rcbL and ETS1f) sequence data (50–54, 95, 96). Important taxonomic distinctions discussed in SI Discussion, Fern Alkyldiols are shown explicitly. Triangle-enclosed digits represent the number of additional branches at different levels of taxonomic classification. CEFA, Cypereae Eleocharideae Fuireneae Abildgaardieae; CSD, Cariceae Scirpeae Dulichieae.Alternative foodstuffs with abrasive, ¹³C-enriched biomass include seedless vascular plants (cryptogams), such as ferns and lycophytes [e.g., quillworts (27, 30)]. Ferns are widely distributed throughout eastern Africa in moist and shaded microhabitats (31) and are often found near dependable sources of drinking water (32). Today, ferns serve as a dietary resource for humans and nonhuman primates alike (27), and fiddlehead consumption is consistent with the inferred digestive physiology [salivary proteins (33)] and the microwear on molars (34) of P. boisei in eastern Africa (25, 26). Ferns were present at FLK^N, based on measurements of C₃₂-diol (Fig. 2D). Further, the high δ¹³C values measured for these compounds are consistent with significant fern consumption by P. boisei at Olduvai Gorge.Ferns and grasses were not the only plant foods present during the time window documented by FLK Zinj. Further, the exclusive reliance on a couple of dietary resources was improbable for P. boisei, because its fossils occur in diverse localities (24–26). Aquatic plants are an additional candidate substrate, as evidenced by high P_aq values at FLK^NN and FLK^N (Fig. 4A). Floating and submerged plants proliferate in wetlands throughout eastern Africa today (13, 14), and many produce nutritious leaves and rootstock all year long (27, 28). Although C₄ photosynthesis is rare among modern macrophytes (30), they can assimilate bicarbonate under alkaline conditions, which results in C₄-like isotope signatures in their biomass (30). Their leaf waxes, such as nC₂₃ (13), are both present and carry ¹³C-enriched signatures at FLK^NN and FLK^N (Fig. 2B). It is also likely that aquatic macrophytes sustained invertebrates and fish with comparably ¹³C-enriched biomass, as they do in modern systems (14), and we suggest aquatic animal foods could have been important in P. boisei diets (27, 28).Biomarkers across the FLK Zinj soil horizon resolve clear patterns in the distribution of plants and water and suggest critical resources that shaped hominin existence at Olduvai Gorge. The behavioral implications of local conditions require understanding of regional climate and biogeography (3–5, 7), because hominin species likely had home ranges much larger than the extent of excavated sites at FLK Zinj. Lake sediments at Olduvai Gorge include numerous stacked tuffs with precise radiometric age constraints (23). These tephrostratigraphic correlations (21) tie the FLK Zinj landscape horizon to published records of plant biomarkers in lake sediments that record climate cycles and catchment-scale variations in ecology. Correlative lake sediment data indicate the wet and wooded microhabitats of FLK Zinj sat within a catchment dominated by arid C₄ grassland (8). Under similarly arid conditions today, only a small fraction of landscape area (ca. 0.05) occurs within 5 km of either forest or standing freshwater (35). Given a paucity of shaded refuge and potable water in the catchment, the concentration of hominin butchery debris (18–21) exclusively within the forest microhabitat and adjacent to a freshwater wetland (Fig. 4) is notable. We suggest the spatial patterns defined by both macro- and molecular fossils reflect hominins engaged in social transport of resources (1–5), such as bringing animal carcasses and freshwater-sourced foods from surrounding grassy or wetland habitats to a wooded patch that provided both physical protection and access to water.     

Insect nicotinic receptor interactions in vivo with neonicotinoid,organophosphorus, and methylcarbamate insecticides and a synergist

The nicotinic acetylcholine (ACh) receptor (nAChR) is the principal insecticide target. Nearly half of the insecticides by number and world market value are neonicotinoids acting as nAChR agonists or organophosphorus (OP) and methylcarbamate (MC) acetylcholinesterase (AChE) inhibitors. There was no previous evidence for in vivo interactions of the nAChR agonists and AChE inhibitors. The nitromethyleneimidazole (NMI) analog of imidacloprid, a highly potent neonicotinoid, was used here as a radioligand, uniquely allowing for direct measurements of house fly (Musca domestica) head nAChR in vivo interactions with various nicotinic agents. Nine neonicotinoids inhibited house fly brain nAChR [³H]NMI binding in vivo, corresponding to their in vitro potency and the poisoning signs or toxicity they produced in intrathoracically treated house flies. Interestingly, nine topically applied OP or MC insecticides or analogs also gave similar results relative to in vivo nAChR binding inhibition and toxicity, but now also correlating with in vivo brain AChE inhibition, indicating that ACh is the ultimate OP- or MC-induced nAChR active agent. These findings on [³H]NMI binding in house fly brain membranes validate the nAChR in vivo target for the neonicotinoids, OPs and MCs. As an exception, the remarkably potent OP neonicotinoid synergist, O-propyl O-(2-propynyl) phenylphosphonate, inhibited nAChR in vivo without the corresponding AChE inhibition, possibly via a reactive ketene metabolite reacting with a critical nucleophile in the cytochrome P450 active site and the nAChR NMI binding site.The nicotinic nervous system has two principal sites of insecticide action, the nicotinic receptor (nAChR) activated by acetylcholine (ACh) and neonicotinoid agonists (1–6), and acetylcholinesterase (AChE) inhibited by organophosphorus (OP) and methylcarbamate (MC) compounds to generate and maintain localized toxic ACh levels (Fig. 1) (7). The nAChR and AChE targets have been identified in insects by multiple techniques but not by direct assays of the ACh binding site in the brain of poisoned insects. Here we use the outstanding insecticidal potency of the nitromethyleneimidazole (NMI) analog of imidacloprid (IMI) (8) as a radioligand (9), designated [³H]NMI, to directly measure the house fly (Musca domestica) nAChR not only in vitro but also in vivo, allowing us to validate by a previously undescribed method the neonicotinoid direct and OP/MC indirect nAChR targets (Fig. 2). This approach also helped solve the intriguing mechanism by which an O-(2-propynyl) phosphorus compound strongly synergizes neonicotinoid insecticidal activity (10) by dual inhibition of cytochrome P450 (CYP) (11–13) and the nAChR agonist site (described herein). Insecticide disruption at the insect nAChR can now be readily studied in vitro and in vivo with a single radioligand allowing better understanding of the action of several principal insecticide chemotypes (Fig. 3).Open in a separate windowFig. 1.The insect nicotinic receptor is the direct or indirect target for neonicotinoids, organophosphorus compounds and methylcarbamates, which make up about 45% of the insecticides by number and world market value (2, 7).Open in a separate windowFig. 2.In this study, Musca nicotinic receptor in vivo interactions with major insecticide chemotypes are revealed by a [³H]NMI radioligand reporter assay. *Position of tritium label.Open in a separate windowFig. 3.Two neonicotinoid nicotinic agonists and two anticholinesterase insecticides.     

Integration of high-content screening and untargeted metabolomics for comprehensive functional annotation of natural product libraries

Traditional natural products discovery using a combination of live/dead screening followed by iterative bioassay-guided fractionation affords no information about compound structure or mode of action until late in the discovery process. This leads to high rates of rediscovery and low probabilities of finding compounds with unique biological and/or chemical properties. By integrating image-based phenotypic screening in HeLa cells with high-resolution untargeted metabolomics analysis, we have developed a new platform, termed Compound Activity Mapping, that is capable of directly predicting the identities and modes of action of bioactive constituents for any complex natural product extract library. This new tool can be used to rapidly identify novel bioactive constituents and provide predictions of compound modes of action directly from primary screening data. This approach inverts the natural products discovery process from the existing ‟grind and find” model to a targeted, hypothesis-driven discovery model where the chemical features and biological function of bioactive metabolites are known early in the screening workflow, and lead compounds can be rationally selected based on biological and/or chemical novelty. We demonstrate the utility of the Compound Activity Mapping platform by combining 10,977 mass spectral features and 58,032 biological measurements from a library of 234 natural products extracts and integrating these two datasets to identify 13 clusters of fractions containing 11 known compound families and four new compounds. Using Compound Activity Mapping we discovered the quinocinnolinomycins, a new family of natural products with a unique carbon skeleton that cause endoplasmic reticulum stress.Notwithstanding the historical importance of natural products in drug discovery (1) the field continues to face a number of challenges that affect the relevance of natural products research in modern biomedical science (2). Among these are the increasing rates of rediscovery of known classes of natural products (3–6) and the high rates of attrition of bioactive natural products in secondary assays due to limited information about compound modes of action in primary whole-cell assays (7). Although pharmaceutical companies recognize that natural products are an important component of drug discovery programs because of the different pharmacologies of natural products and synthetic compounds (8), there is a reluctance to return to “grind and find” discovery methods (9). Therefore, there is a strong need for technologies that address these issues and provide new strategies for the prioritization of lead compounds with unique structural and/or biological properties (10).Natural product drug discovery is challenging in any assay system because extract libraries are typically complex mixtures of small molecules in varying titers, making it difficult to distinguish biological outcomes (11). This is compounded by issues of additive effects of multiple bioactive compounds and the presence of nuisance compounds that cause false positives in assay systems (12). To address these issues, our laboratory has recently developed several image-based screening platforms that are optimized for natural product discovery (13–16). The cytological profiling platform optimized by Schulze and coworkers characterizes the biological activities of extracts using untargeted phenotypic profiling. These phenotypic profiles are compared with natural products extracts and a training set of compounds with known modes of action to characterize the bioactivity landscape of the screening library (17, 18). This cytological profiling tool forms the basis of the biological characterization component of the Compound Activity Mapping platform, as described below.In the area of chemical characterization of natural product libraries, untargeted metabolomics is gaining attention as a method for evaluating chemical constitution (3, 19–22). Modern “genes-to-molecules” and untargeted metabolomics approaches taking advantage of principal component analysis and MS² spectral comparisons have also been developed to quickly dereplicate complex extracts and distinguish noise and nuisance compounds from new molecules (23–27). Unfortunately, although these techniques are well suited to the discovery of new chemical scaffolds, they are unable to describe the function or biological activities of the compounds they identify. Therefore, there is still a need for new approaches to systematically identify novel bioactive scaffolds from complex mixtures.To overcome some of these outstanding challenges we have developed the Compound Activity Mapping platform to integrate phenotypic screening information from the cytological profiling assay with untargeted metabolomics data from the extract library (Fig. 1). By correlating individual mass signals with specific phenotypes from the high-content cell-based screen (Fig. 2), Compound Activity Mapping allows the prediction of the identities and modes of action of biologically active molecules directly from complex mixtures, providing a mechanism for rational lead selection based on desirable biological and/or chemical properties. To evaluate this platform for natural products discovery we examined a 234-member extract library, from which we derived 58,032 biological measurements (Fig. 1C) and 10,977 mass spectral features (Fig. 1A). By integrating and visualizing these data we created a Compound Activity Map for this library composed of 13 clusters containing 16 compounds from 11 compound classes (Fig. 3). This integrated data network enabled the discovery of four new compounds, quinocinnolinomycins A–D (1–4, Fig. 4), which are the first examples to our knowledge of microbial natural products containing the unusual cinnoline core (Fig. 5). Clustering the cytological profiles of the quinocinnolinomycins with those of the Enzo library training set suggests that these compounds induce endoplasmic reticulum (ER) stress and the protein unfolding response.Open in a separate windowFig. 1.Overview of Compound Activity Mapping. (A) Representation of the chemical space in the tested extract library. The network displays extracts (light blue) connected by edges to all m/z features (red) observed from the metabolomics analysis, illustrating the chemical complexity of even small natural product libraries. (B) Histograms of activity and cluster scores for all m/z features with cutoffs indicated as red lines (for full-size histograms see SI Appendix, Fig. S5). (C) Compound Activity Map, with the network displaying only the m/z features predicted to be associated with consistent bioactivity, and their connectivity to extracts within the library. (D) Expansion of the staurosporine cluster (dotted box in C) with extract numbers and relevant m/z features labeled.Open in a separate windowFig. 2.Determination of synthetic fingerprints and cluster and activity scores. (A) Table of Pearson correlations for the cytological profiles between all extracts containing a specific m/z feature (m/z of 489.1896, rt of 1.59). In each cytological profile, yellow stripes correspond to positive perturbations in the observed cytological attribute and blue stripes correspond to negatively perturbed attributes. The cluster score is determined by calculating the average of the Pearson correlation scores for all relevant extracts. (B) Calculated synthetic fingerprint and activity score for feature (m/z of 489.1896, rt of 1.59). Synthetic fingerprints are calculated as the averages of the values for each cytological attribute to give a predicted cytological profile for each bioactive m/z feature in the screening set.Open in a separate windowFig. 3.Annotated Compound Activity Map. An expanded view of the Compound Activity Map from Fig. 1C, with the extracts and m/z features separated into subclusters and colored coded using the Gephi modularity function. Each bioactive subcluster is composed of extracts containing a family of compounds with a defined biological activity. The Compound Activity Map is annotated with a representative molecule from each of the families of compounds that have been independently confirmed by purification and chemical analysis.Open in a separate windowFig. 4.The prioritization, isolation, and confirmation of the quinocinnolinomycins A–D (1–4). (A) Bioactive m/z features plotted on a graph of activity score vs. cluster score. The color of the dot corresponds to the retention time of the m/z feature with the color bar and scale below in minutes. (B) Isolated cluster from Fig. 1C and Fig. 3 containing both the relevant extracts (blue) and bioactive m/z features (red). (C) HPLC trace of extract RLPA-2003E and the isolation of quinocinnolinomycins A–D (highlighted with blue boxes on HPLC trace). (D) Cell images of pure compounds screened as a twofold dilution series for quinocinnolinomycins A and B in both stain sets compared with images of vehicle (DMSO) wells. (E) Comparison of the synthetic and actual cytological fingerprints of the pure compounds is presented below the relevant images, demonstrating the relationship between experimental and calculated cytological profiles for these two metabolites.Open in a separate windowFig. 5.Structure elucidation of quinocinnolinomycins A–D (1–4). (A) Structures of quinocinnolinomycins A–D. (B) Key NMR correlations used in the structure elucidation of quinocinnolinomycin A. COSY correlations are indicated by bold lines. Heteronuclear multiple-bond correlations are indicated by curved arrows. (C) ∆δ^SR values for the Mosher’s α-methoxy-α-trifluoromethylphenylacetic acid (MTPA) ester analysis of the secondary alcohol in quinocinnolinomycin A (1) to assign the absolute configuration at position C11.     

Pictet–Spengler reaction-based biosynthetic machinery in fungi

The Pictet–Spengler (PS) reaction constructs plant alkaloids such as morphine and camptothecin, but it has not yet been noticed in the fungal kingdom. Here, a silent fungal Pictet–Spenglerase (FPS) gene of Chaetomium globosum 1C51 residing in Epinephelus drummondhayi guts is described and ascertained to be activable by 1-methyl-l-tryptophan (1-MT). The activated FPS expression enables the PS reaction between 1-MT and flavipin (fungal aldehyde) to form “unnatural” natural products with unprecedented skeletons, of which chaetoglines B and F are potently antibacterial with the latter inhibiting acetylcholinesterase. A gene-implied enzyme inhibition (GIEI) strategy has been introduced to address the key steps for PS product diversifications. In aggregation, the work designs and validates an innovative approach that can activate the PS reaction-based fungal biosynthetic machinery to produce unpredictable compounds of unusual and novel structure valuable for new biology and biomedicine.Microbes and plants produce a multitude of unpredictably structured organic molecules known as secondary metabolites (natural products), from which more than half of globally marketed drugs have been developed (1–3). Large-scale genomic mining has indicated that microbial secondary metabolites are substantially underestimated because many biosynthetic genes remain silent or less active in the laboratory cultivation conditions (4, 5). Accordingly, there has long been an urgent need to develop a new strategy that enables microorganisms to produce more unforeseeable bioactive compounds, which are important to the drug discovery efforts to combat life-threatening diseases (6, 7), and to the complexity-based driving force for synthetic and material chemistry (8–10).Characterized by forming a piperidine ring through a condensation of β-arylethylamine with an aldehyde, Pictet–Spengler (PS) reaction contributes greatly to the framework diversification of important alkaloidal phytochemicals such as morphine, camptothecin, and reserpine (Fig. 1A), with plant-derived Pictet–Spenglerase (called strictosidine synthetase, STR) mechanistically addressed (11). The PS mechanism has been presumed to involve in the tetrahydroisoquinoline antibiotic biosynthesis in the bacterium Streptomyces lavendulae (12, 13), and likely in the biosynthetic pathway of hyrtioreticulin F in the marine sponge Hyrtios reticulatus (14). However, surprisingly, nothing is known concerning the PS reaction in the fungal kingdom.Open in a separate windowFig. 1.Alkaloids derived from the PS reaction. (A) Representatives for PS reaction-based phytochemicals. (B) 1-MT has been found to be a potent up-regulator for the FPS gene expression, and a suitable FPS substrate for its PS condensation with flavipin to yield unnatural natural products (1−8) with unprecedented skeletons.Most if not all Chaetomium fungi in the Chaetomiaceae family produce l-tryptophan–derived alkaloids, but “refuse” to generate any PS reaction-based secondary metabolite (15–18). However, a comparative genomic analysis has clarified that C. globosum 1C51 does have an FPS gene (CHGG_06703, STR-like) (SI Appendix, Fig. S25), but remains silent or poorly activated in the laboratory cultivations because no PS-derived secondary metabolite has been detected in the fungal culture. Therefore, this C. globosum 1C51 strain was adopted here to test for the activation of its “unworking” PS reaction-based biosynthetic machinery. As a result, 1-methyl-l-tryptophan (1-MT) was demonstrated to be able to up-regulate the FPS expression and condense with the fungal aldehyde flavipin (3,4,5-trihydroxy-6-methyl phthalaldehyde) to form unexpectedly a family of skeletally unprecedented alkaloids, trivially named chaetoglines A−H (1−8) (Fig. 1B). A gene-implied enzyme inhibition (GIEI) strategy, derived from the hypothesis-based enzyme modulation described elsewhere (19, 20), was introduced to identify the key diversification steps for the PS reaction-derived compounds (Figs. 2–4). Chaetoglines B (2) and F (6) have been found to be more antibacterial than tinidazole (a coassayed drug prescribed in clinic for bacterial infections) against pathogenic anaerobes Veillonella parvula, Bacteroides vulgatus, Streptococcus sp., and Peptostreptococcus sp. Moreover, alkaloid 6 is potently inhibitory on acetylcholinesterase (AChE), an effective target enzyme exploited for the treatment of Alzheimer’s disease (21, 22).Open in a separate windowFig. 2.LC-MS profile-based comparisons for the chaetogline production in monooxygenase inhibitor exposed fungal cultures: (A) for 3–6; (B) for 1–2 and 7–8. The ESI-MS spectra (C) of 1−8 (①) displayed the corresponding protonated and Na⁺-liganded molecular ions. Samples were ethyl acetate extracts derived from 1-MT supplemented C. globosum 1C51 cultures without (②) and with the separate exposure to PB, PR, and MMI at 0.1 (③∼⑤) and 1.0 mM (⑥∼⑧), respectively.Open in a separate windowFig. 4.Proposed generation of the fungal PS-derived products. Catalyzed by FPS, 1-MT and flavipin (a fungal aldehyde) undergo PS reaction to form chaetoglines A−H (1−8) in concert with tailing reactions including oxidation, decarboxylation, and Aldol reaction. The Schiff base intermediate 9 tends to tautomerize via 10 and 11 to give chaetogline C (3) that can be methyl-esterified into chaetogline D (4). Intramolecular cyclization of 10 gives 12, which is oxidizable into chaetoglines A (1) and E (5), the latter yielding chaetogline F (6) after the oxidative and decarboxylative aromatization. Chaetogline E (5) can also be oxidized to intermediate 13, which gives 14 after condensing presumably via the decarboxylative Aldol reaction (34, 35) with 1-M-IAA derived from 1-MT by the fungi (SI Appendix, Fig. S24). Intermediate 14 undergoes intramolecular cyclization, monooxygenation, and isomerization to form chaetogline G (7), which after decarboxylation gives chaetogline H (8), a precursor of chaetogline B (2).     

Rare crested rat subfossils unveil Afro–Eurasian ecological corridors synchronous with early human dispersals

Biotic interactions between Africa and Eurasia across the Levant have invoked particular attention among scientists aiming to unravel early human dispersals. However, it remains unclear whether behavioral capacities enabled early modern humans to surpass the Saharo–Arabian deserts or if climatic changes triggered punctuated dispersals out of Africa. Here, we report an unusual subfossil assemblage discovered in a Judean Desert’s cliff cave near the Dead Sea and dated to between ∼42,000 and at least 103,000 y ago. Paleogenomic and morphological comparisons indicate that the specimens belong to an extinct subspecies of the eastern African crested rat, Lophiomys imhausi maremortum subspecies nova, which diverged from the modern eastern African populations in the late Middle Pleistocene ∼226,000 to 165,000 y ago. The reported paleomitogenome is the oldest so far in the Levant, opening the door for future paleoDNA analyses in the region. Species distribution modeling points to the presence of continuous habitat corridors connecting eastern Africa with the Levant during the Last Interglacial ∼129,000 to 116,000 y ago, providing further evidence of the northern ingression of African biomes into Eurasia and reinforcing previous suggestions of the critical role of climate change in Late Pleistocene intercontinental biogeography. Furthermore, our study complements other paleoenvironmental proxies with local—instead of interregional—paleoenvironmental data, opening an unprecedented window into the Dead Sea rift paleolandscape.

Situated at the gateway of Africa, the Levant witnessed major Afro–Eurasian biotic exchanges during the Neogene-Quaternary (1–3), including multiple hominin dispersal events (4–11). The extent and timing of potential paleoenvironmental connections between these regions during the Late Pleistocene have received considerable attention due to their role in the global expansion of Homo sapiens out of Africa (9, 12–17). Due to multiple dating uncertainties, the paucity of the fossil record, and the disparate scale and resolution of paleoenvironmental proxies (17), it remains unclear whether technological and behavioral capacities enabled early modern humans to surpass the Saharo–Arabian biogeographic barrier (6, 18–20) or if climatic changes triggered punctuated dispersals out of Africa through the creation of ecological corridors (14–16, 21, 22). The southern part of the Levant is a geographical bottleneck, longitudinally divided by mountain ranges into a wooded Mediterranean zone and the more arid, inland regions of the Dead Sea rift. Under monsoon-dominated climatic conditions, these arid regions could have supported savanna-like environments that funnelled the dispersal of African faunas by linking the Arabian and Sinai regions with the eastern Mediterranean (16, 23). While the paleoenvironmental interpretation of the well-known Mediterranean Levantine fossil record continues to be scrutinized and debated (16, 17, 24–26), there is no record of faunal movement in the Dead Sea arid areas during the Late Pleistocene (1, 27).Here, we report evidence that the southeastern Levant, today blocked by a rain shadow desert, had more densely vegetated habitats in the Late Pleistocene that supported African faunal immigrants synchronously with early human dispersals. Recent excavations (28) and surveys (29) in the Cave of the Skulls (CoS), in the southern Judean Desert, have yielded exceptionally well-preserved fossils of a subspecies of the eastern African crested rat, Lophiomys imhausi maremortum subspecies nova (subsp. nov.) (Figs. 1–3), dated to between ∼42,000 y ago (ka) and at least 103 ka. The subspecific status is sustained by paleogenomic data, which indicate low molecular divergence (<2%) between the Levantine population and the extant populations from eastern Africa. A few specimens dated to >44 and ∼112 ka from Sodmein Cave (SOD), in northeastern Egypt [Fig. 1A (30)], are also ascribed to this subspecies.Open in a separate windowFig. 1.Geographical setting of this study. (A) Map of Africa and the Middle East showing fossil sites with Lophiomys: 1, Irhoud Ocre (Jebel Irhoud); 2, Khemis; 3, Lissasfa; 4, Salobreña; 5, Amama; 6, Oued Ahtmenia; 7, Sheikh Abdallah; 8, SOD; 9, Alayla, Middle Awash; and 10, CoS. (B) Paleoartistic reconstruction of L. i. maremortum subsp. nov. Artwork by Aya Marck, used with permission. (C) Elevation map showing the location of the CoS. (D) Planar map of the CoS, showing the different sectors excavated (sectors A to Q) and the areas with higher density of L. i. maremortum subsp. nov. (E) The CoS main entrance, opening to a vertical cliff in Nahal Tze’elim. During excavations, the cave had to be approached with the use of ropes. Photograph by Guy Fitoussi (Israel Antiquities Authority), used with permission. (F) A complete skull of L. i. maremortum subsp. nov. (SMNH-TAU {"type":"entrez-nucleotide","attrs":{"text":"M17121","term_id":"1015602695"}}M17121) found in situ in sector Q during a survey of the cave. Note that the hemimandibles are in anatomical connection. Photograph by I.A.L.Open in a separate windowFig. 3.Dental metrics and first upper molar geometric morphometrics of Lophiomys imhausi maremortum subsp. nov in comparison to extant populations of L. imhausi and extinct relatives. (A) Dental metrics of upper (M1 to M3) and lower (m1 to m3) molars. Note that the specimens of L. i. maremortum subsp.nov. (CoS, in blue) are consistently larger than the extant populations (northeastern, NE, in orange; southwestern, SW, in green). (B) Geometric morphometric of the first upper molar (M1) visualized by means of a between-group principal component analysis. Centroid size is calculated as the sqrt of the sum of squared distances of all landmarks from the center of the shape before Procrustes analysis. Independently of size, the three populations are well distinguished by M1 shape.The crested rat is a large (>800 g) eastern African rodent that shows some of the most extraordinary adaptive features among living mammals, including a granulated helmet-like skull, a poisonous pelt, and a three-chambered stomach with hindgut fermentation (31, 32). The crested rat has the ability to become toxic to predators by chewing up the roots and bark of the poison-arrow tree (Acokanthera ssp.) and using its mouth to spread the toxin—a cardiac glycoside called ouabain—on strips of specialized hairs on its body flanks (33, 34). In general, the crested rat inhabits steep, rocky valleys within woodlands and montane forests under variable but typically seasonal rainfall, and its diet includes leaves, fruits, and shoots (32). Populations are apparently stable in relatively wet and densely vegetated habitats in the montane forests of Kenya and the highlands of Ethiopia, referred to here as southwest populations (SW; Fig. 4A), where the mean annual precipitation can reach 2,500 mm (35). It is also known from drier and more sparsely vegetated mountain ranges or lowland woodlands in Somalia, Djibouti, Eritrea, and southeastern Sudan, referred here as northeast populations (NE; Fig. 4A), where the mean annual precipitation can be as little as 350 mm (36). Despite its apparently broad climatic tolerance, the crested rat is rare in the wild, and it usually prefers riparian vegetation, where it can hide and den in holes under trees (34). Therefore, it is likely that the potential distribution and dispersal capability of Lophiomys is limited by the availability of trees and by vegetation density.Open in a separate windowFig. 4.Estimated habitat corridors derived from SDMs of Lophiomys in relation to paleoclimatic, paleobotanical, and archaeological evidence in the Late Pleistocene. (A) Map of the study region colored by NDVI, with areas in green indicating more vegetation cover and areas in orange showing sparsely vegetated areas. The circles are confirmed occurrences of extant L. imhausi. The southwest (SW) population is colored in green, and the northeast (NE) population is colored in orange. (B) Climatic niche envelope of extant L. imhausi SW and NE populations and at the location of the CoS with the present climate. The variables are the same used in the SDMs. Units have been scaled to fit the chart. For an explanation of variables, reference SI Appendix, Text S13 and SI Appendix, Table S17. (C–E) Maps of the study region showing potential habitat suitability areas for Lophiomys based on Maxent-based species distribution modeling during (C) current time, (D) LGM, and (E) the LIG. The red colors indicate higher habitat suitability estimates. (F) Diagram showing δ¹⁸O curves in the Southern Levant and pollen records of the Dead Sea. From Left to Right, chronology and δ¹⁸O record of Soreq, Peqi’in, and Tzavoa Caves (data from refs. 47 and 64), arboreal and sclerophyllous pollen concentration, Pistacia pollen concentration, and Artemisia and Amaranthaceae pollen concentration (data from refs. 23 and 46). (G) A photograph of Nahal Tze’elim taken from the CoS main entrance, with the Dead Sea on the top of the image. Photograph by I.A.L. (H) Photo of Iran’s Zagros Mountains. Photograph by Marijn van den Brink, distributed under a CC-BY 3.0 license. (I) A photo of Djibouti’s National Day Park, showing how the Judean Desert could have looked like during part of the Late Pleistocene. Photograph by Steven Dessein, distributed under a CC-BY 3.0 license. (J) Map of the southern Levant showing current mean annual precipitation (in blue) and Late Pleistocene archaeological sites (yellow squares). The spatial arrangement of sites suggests that the nowadays-arid regions—including the Dead Sea basin—could have funneled early human dispersals into the Levant. Precipitation curves are redrawn from ref. 65.The evolutionary and biogeographic history of crested rats is poorly documented. Nowadays, the Lophiomyinae is represented by only one eastern African species. However, the group was more diverse in the past, and its distribution spread over areas of eastern and northern Africa and Spain (Fig. 1; see SI Appendix, Table S1, for a list of all sites and references). The subfossils described here represent the most extensive paleontological collection (n > 250) of Lophiomys ever found. Based on the known habitats occupied by extant Lophiomys and paleoclimatic data, this finding suggests that continuous habitat corridors connected eastern Africa with the Levant during the Last Interglacial (LIG, ∼129 to 116 ka) along the western Red Sea mountains. In contrast to the present-day hyper-arid conditions, the Late Pleistocene Judean Desert was likely wetter and included a variable but relatively persistent arboreal component. Altogether, our study provides further evidence of the expansion of African biomes into Eurasia and the role of habitat corridors in early hominin dispersals out of Africa.     

High thermoelectric performance by resonant dopant indium in nanostructured SnTe

From an environmental perspective, lead-free SnTe would be preferable for solid-state waste heat recovery if its thermoelectric figure-of-merit could be brought close to that of the lead-containing chalcogenides. In this work, we studied the thermoelectric properties of nanostructured SnTe with different dopants, and found indium-doped SnTe showed extraordinarily large Seebeck coefficients that cannot be explained properly by the conventional two-valence band model. We attributed this enhancement of Seebeck coefficients to resonant levels created by the indium impurities inside the valence band, supported by the first-principles simulations. This, together with the lower thermal conductivity resulting from the decreased grain size by ball milling and hot pressing, improved both the peak and average nondimensional figure-of-merit (ZT) significantly. A peak ZT of ∼1.1 was obtained in 0.25 atom % In-doped SnTe at about 873 K.Good thermoelectric (TE) materials should not only have high figure-of-merit (Z), but also be environmentally friendly and cost-effective (1–5). The nondimensional figure-of-merit (ZT) is defined as ZT = [S²σ/(κ_L+κ_e)]T, where S is the Seebeck coefficient, σ the electrical conductivity, κ_L the lattice thermal conductivity, κ_e the electronic thermal conductivity, and T the absolute temperature. Lead chalcogenides and their alloys can be engineered to exhibit high ZTs; however, environmental concern regarding Pb prevents their deployment in large-scale applications (6–10). Tin telluride (SnTe), a lead-free IV–VI narrow band-gap semiconductor has not been considered favorably as a good thermoelectric material because of its low ZT due to the relatively low Seebeck coefficient and high electronic thermal conductivity caused by intrinsic Sn vacancies (11–13), although SnTe has been used to alloy with other tellurides for better TE properties (14–26). Even though there has been no real success in achieving good TE properties of lead-free SnTe, the similarity between the electronic band structure of SnTe and that of PbTe and PbSe (27–31) suggests it has the potential to be a good TE material, especially given the two valence bands (light-hole and heavy-hole bands) that contribute to the hole density of states. The main difficulty here, however, is the fact that the separation between the light-hole and heavy-hole band edges in SnTe is estimated to be in the range of ∼0.3 to ∼0.4 eV (27, 29), larger than those of PbTe or PbSe (9), rendering the benefit of the heavier mass for the Seebeck coefficient less significant.In this paper, we prepared In-doped SnTe by high-energy ball milling and hot pressing and measured the samples up to 873 K without experiencing any mechanical strength issues. We show, based on both experiments and first-principles simulation, that a small amount of In-doping helps create resonant states around the Fermi level inside the valence band, which increases the Seebeck coefficient, especially at room temperature, leading to improvements in both average ZT and peak ZT, combined with the decreased lattice thermal conductivity due to the increased density of grain boundaries (32–34). Peak ZT value reaches ∼1.1 at about 873 K for SnTe doped with 0.25 atom % In.Single-phased In-doped SnTe was obtained by ball milling and hot pressing. Fig. 1 presents the X-ray diffraction (XRD) patterns of In_xSn_1-xTe (x = 0, 0.0025, 0.005, and 0.01). All the peaks can be indexed to the face-centered structure (space group Fmm). No impurity phase was found, despite the increasing content of In. First-principles calculations (Table S1) indicated it is energetically favorable for In to substitute for Sn, which is consistent with the case in In-doped PbTe and PbSe. In previous work, we found In substitutes for Pb in PbTe and PbSe, which is the same with In-doped SnTe, but it is n-type doping in In_xPb_1-xTe and In_xPb_1-xSe, which is different from p-type doping by In in SnTe, as we are reporting in this work (35, 36).Open in a separate windowFig. 1.XRD patterns for In_xSn_1-xTe (x = 0, 0.0025, 0.005, and 0.01) prepared by ball milling and hot pressing.The electrical conductivities decrease with increasing temperature, as shown in Fig. 2A, showing the typical behavior of degenerate semiconductors. With increasing content of In, the electrical conductivity decreases, especially at room temperature, from ∼7 × 10⁵ S⋅m⁻¹ to ∼2 ×10⁵ S⋅m⁻¹. The hole concentration indicated by the Hall measurement, however, changes in an interesting way with increasing In content: it drops below the intrinsic value at the beginning and starts to rise after x ≥ 0.0025 (as shown in Fig. 3A). Based on this observation, we conclude In atoms should be p-type dopants and explain the change of the carrier concentration as follows. The intrinsic SnTe is p-type because of the Sn vacancies (19). Those vacancies create empty electronic states and behave like p-type dopants. If we dope SnTe with In, In atoms first fill the Sn vacancies. Despite being p-type dopants, they are not as “strong” as the vacancies, in the sense that they induce fewer holes (examined by the simulation shown in Table S1); thus, at low doping levels, the p-type charge concentration decreases. However, as the doping level is increased, at some point all the Sn vacancies are filled with In, and beyond that point, excessive In atoms substitute for Sn, and the p-type charge concentration increases again (Fig. 3A). However, when In is more than the solubility limit in SnTe, the extra In atoms act as donors, which decreases the hole carrier concentration (x = 0.01) (37). The fact that the electrical conductivity decreases all the way indicates that the In dopants affected the hole mobility significantly (shown in Fig. 3B), as the result of both increased effective mass and impurity scattering. The Seebeck coefficients increase with temperature in the whole temperature range and also increase with In content, as shown in Fig. 2B. No bipolar effect is evident, even up to 873 K, in all the compositions despite the small band gap ∼0.18 eV for SnTe (29, 31). All the measured Seebeck coefficients are positive, consistent with the density of states (DOS) calculation presented in Fig. 4 and different from In-doped PbTe and PbSe, in which In turned out to be an n-type dopant (36, 38). Fig. 2C shows the power factors for undoped and In-doped SnTe. The highest power factor reaches ∼2.0 × 10⁻³ W⋅m⁻¹⋅K⁻² at about 873 K, higher than all the reported power factors of doped PbTe and PbSe at this temperature (9, 39–41). Most importantly, the average power factor is increased a great deal by In doping. Compared with the undoped SnTe prepared by melting and hand milling (M+HM) (broken line), the electrical properties of the ball-milled samples are not different.Open in a separate windowFig. 2.Temperature dependence of (A) electrical conductivity, (B) the Seebeck coefficient, and (C) the power factor for In_xSn_1-xTe (x = 0, 0.0025, 0.005, and 0.01). The undoped SnTe prepared by melting, hand milling, and hot pressing (M+HM) is shown for comparison (broken line).Open in a separate windowFig. 3.Hall carrier concentration (A) and Hall mobility (B) at room temperature with respect to the doping content x. ○, undoped SnTe; ●, In-doped SnTe.Open in a separate windowFig. 4.Comparison of DOS for undoped SnTe (broken line), Bi-doped SnTe (solid line), and In-doped SnTe (bold solid line). Sharp features are observed in the DOS of In-doped SnTe near the band edge, to which the abnormal Seebeck coefficient might be attributed. The simulated supercell configuration corresponds to 3 atom % In concentration, which is higher than that achieved in the experiment. The Fermi level in the simulation resides at 6.207 eV, slightly below the DOS hump. With the experimental In concentration, the Fermi level is expected to sit around the DOS peak.Fig. 5 shows variation of the Seebeck coefficient vs. carrier concentration for both pure SnTe and In-doped SnTe. The Seebeck coefficients of undoped SnTe with different hole concentrations (2 × 10²⁰ to 1.8 × 10²¹ cm⁻³) were obtained previously by annealing under different conditions (open circles) (27). The carrier concentration obtained in this work is ∼2.35 × 10²⁰ cm⁻³ (filled circle). Unlike PbTe and PbSe (7, 9, 36, 39, 40), the Seebeck coefficient of SnTe shows abnormal variation with increasing carrier concentration, which was qualitatively explained previously by two parabolic band models (27) and density functional theory (DFT) calculations (31). The valence band model (VBM), which takes into account the nonparabolicity of the light-hole band (solid line), provides a quantitative fit to all the Seebeck coefficient data, except for those of In-doped samples, and thus is expected to best depict the contribution from the intrinsic band structure of SnTe (29). The model details for TE transport of p-type SnTe may be found in SI Text. Compared with the same model we used for PbTe and PbSe (9, 36), two major differences should be stated. The L point energy gap, Eg, is smaller for SnTe, making the nonparabolicity larger. This makes the Seebeck coefficient drop faster with increasing concentration, as seen in Fig. S1. The light-hole–heavy-hole band edge energy difference is 0.12 eV for PbTe, 0.26 eV for PbSe, and 0.35 eV for SnTe (9, 29, 36); thus, the heavy-hole contribution is relatively weaker for SnTe. This may be seen from the fact that there is not much difference between the predictions of VBM and those of the two-band Kane model (which ignores the heavy-hole band contribution) at room temperature for SnTe, until 10 × 10¹⁹ cm⁻³. However, the contribution from the heavy-hole band gradually increases at higher temperatures (Fig. S2) as for PbSe (9, 36), helping improve the Seebeck coefficient at high temperature and suppress the bipolar effect. Although the Seebeck coefficients of bismuth- (Bi-) and Cu-doped samples agree well with the VBM model, as shown in Fig. 5, indicating pure doping effects, the deviation of the In-doped samples from the VBM model implies that there must be mechanisms through which In dopants significantly alter the band structure of pure SnTe near the band edge. One of the possible mechanisms is the introduction of resonant levels (6, 42–44) into the valence band. Fig. 4 shows the DOS of pure SnTe, Bi-doped SnTe, and In-doped SnTe near the top of the valence band. A well-defined peak is observed in the DOS of In-doped SnTe that may contribute to the large deviation of the Seebeck coefficient from the VBM model. One may question whether the observed features are a result of the limited size of the supercell and thus the artificial interactions between In atoms. Similar features, however, are not observed in Bi-doped SnTe with the same supercell size. Therefore, we believe the added feature originates from the interactions of the In atoms with the host atoms. Because of the limitation of computing resources, a sufficiently dense k-mesh for calculating transport properties for the supercells is not possible at this stage; also, the simulated supercells are too small to represent a realistic doping concentration. [The simulated supercell corresponds to 3% In concentration, with a Fermi level located slightly below the DOS “hump.” With the doping concentration achieved in the experiments, the Fermi level is expected to reside close to the DOS peak. An alternative simulation method, such as a Korringa–Kohn–Rostoker coherent-potential-approximation (KKR-CPA) calculation (44), is required in cases of more dilute doping concentrations.] Thus, a direct evaluation of the effect of the features in DOS on the Seebeck coefficient is not available for now. However, the rich features introduced by In atoms are speculated to play an important role in the enhanced TE properties.Open in a separate windowFig. 5.Room temperature Pisarenko plot for ball-milled and hot-pressed In_xSn_1-xTe (x = 0, shown by ●; x = 0.001, 0.0015, 0.0025, 0.005, 0.0075, and 0.01, shown by ▲) in comparison with reported data on undoped SnTe (○), Bi-doped SnTe (□), and Cu-doped SnTe (♢) by Brebrick and Strauss (27). The solid curve is based on the VBM (light nonparabolic band and heavy parabolic band) with the heavy-hole effective mass of SnTe m*/m_e = 1.92.The other problem we should resolve is the high thermal conductivity induced by intrinsic Sn vacancies, causing very high electrical conductivity. By In doping, the decreased electrical conductivity results in a reduced electronic part of the thermal conductivity determined by the Wiedemann–Franz law (κ_e = LσT), where L is the Lorenz number. The Lorenz number is calculated using the VBM in a way similar to that of the Seebeck coefficient, including contributions from both nonparabolic light-hole and parabolic heavy-hole bands. The detailed expressions used are included in SI Text. Fig. 6 A–C gives the temperature dependences of the thermal diffusivity, specific heat, total thermal conductivity, and lattice thermal conductivity (obtained by subtracting the electronic contribution from the total thermal conductivity) of the undoped and In-doped SnTe, respectively. With increasing temperature, the total thermal conductivity decreases rapidly without showing any bipolar effect, consistent with the behavior of the Seebeck coefficient in Fig. 2B. The total thermal conductivities of all In-doped SnTe are lower than the undoped sample. Compared with the undoped SnTe prepared by melting and hot pressing (dotted line), the samples prepared by ball milling and hot pressing exhibit lower lattice thermal conductivity, which may be attributed to the increased density of grain boundaries by ball milling. In Fig. 7, the representative microstructure of ball-milled and hot-pressed In-doped SnTe is presented. Scanning electron microscopic (SEM) images shown in Fig. 7A indicate that the In_0.0025Sn_0.9975Te samples consist of both big grains with diameters of several tens of microns and small grains. The observed small cavities may contribute to the lower lattice thermal conductivity. The densities of all the samples are listed in Table S2. The size of the small grains is about 1 μm, as shown in Fig. 7B, less than one tenth that of the big grains. Nanograins in the samples also are observed via transmission electron microscopy (TEM). Fig. 7C shows a typical bright-field TEM image of the nanograins, with sizes around 100 nm. As a result, the lattice thermal conductivity of the samples is greatly reduced by significantly enhanced boundary scatterings of the phonons, as shown in Fig. 6C. Selected area electron diffraction and high-resolution TEM (HRTEM) images show that all the grains, whether in microns or nanometers, are single crystals with clean boundaries and good crystallinity, as shown in Fig. 7D. The crystalline grains and boundaries would benefit the transport of charge carriers, as observed in nanograined Bi_xSb_2-xTe₃ bulks (45), without degrading the electronic properties (Fig. 2).Open in a separate windowFig. 6.Temperature dependence of (A) thermal diffusivity (the undoped SnTe prepared by melting and hot pressing is shown by the broken line), (B) specific heat (the specific heat of sample x = 0 is used for the undoped SnTe prepared by melting and hot pressing), and (C) total thermal conductivity and lattice thermal conductivity for In_xSn_1-xTe (x = 0, 0.0025, 0.005, and 0.01) (the undoped SnTe prepared by melting and hot pressing is shown by the broken line).Open in a separate windowFig. 7.Representative SEM (A and B), TEM (C), and HRTEM (D) images for as-prepared In_0.0025Sn_0.9975Te samples by ball milling and hot pressing.Fig. 8 summarizes the ZT values of different samples. The two intrinsic valence bands contribute to the peak ZT value ∼0.7 at about 873 K for the undoped SnTe. The decreased lattice thermal conductivity by ball milling further boosts the peak ZT value to ∼0.8. However, the ZT values in both cases are quite low, below 600 K, resulting in low average ZTs. The enhanced Seebeck coefficient by resonant states increased both the peak and average ZTs in the In-doped nanostructured SnTe. A peak ZT ∼1.1 is obtained at about 873 K in In_0.0025Sn_0.9975Te.Open in a separate windowFig. 8.Temperature dependence of ZT for In_xSn_1-xTe (x = 0, 0.0025, 0.005, and 0.01) compared with the reported data on undoped SnTe (0.5–2.0 atom % Te excess) (■) by Vedeneev et al. (25) and codoped SnTe (0.5–2.0 atom % Te excess) with 1 atom % In and 1 atom % Ag (●) by Vedeneev et al. (25). The undoped SnTe prepared by melting and hot pressing is included for comparison (broken line).In summary, nanostructured In-doped SnTe with a ZT >1 has been prepared by ball milling and hot pressing. The improved ZT (peaked around 1.1 at about 873 K in 0.25 atom % In-doped SnTe) incorporates both the high Seebeck coefficient resulting from the two valence bands and the local resonant states around Fermi level created by In-doping and the lowered lattice thermal conductivity owing to the increased phonon interface scattering. This lead-free TE material is a potential candidate to replace lead chalcogenides used at medium to high temperatures for waste heat recovery applications. Further improvement is expected by adding suitable nanoinclusions or alloying with SnSe and SnS to decrease the thermal conductivity and increase the Seebeck coefficient.     

From the Cover: Environmental drivers of annual population fluctuations in a trans-Saharan insect migrant

Many latitudinal insect migrants including agricultural pests, disease vectors, and beneficial species show huge fluctuations in the year-to-year abundance of spring immigrants reaching temperate zones. It is widely believed that this variation is driven by climatic conditions in the winter-breeding regions, but evidence is lacking. We identified the environmental drivers of the annual population dynamics of a cosmopolitan migrant butterfly (the painted lady Vanessa cardui) using a combination of long-term monitoring and climate and atmospheric data within the western part of its Afro-Palearctic migratory range. Our population models show that a combination of high winter NDVI (normalized difference vegetation index) in the Savanna/Sahel of sub-Saharan Africa, high spring NDVI in the Maghreb of North Africa, and frequent favorably directed tailwinds during migration periods are the three most important drivers of the size of the immigration to western Europe, while our atmospheric trajectory simulations demonstrate regular opportunities for wind-borne trans-Saharan movements. The effects of sub-Saharan vegetative productivity and wind conditions confirm that painted lady populations on either side of the Sahara are linked by regular mass migrations, making this the longest annual insect migration circuit so far known. Our results provide a quantification of the environmental drivers of large annual population fluctuations of an insect migrant and hold much promise for predicting invasions of migrant insect pests, disease vectors, and beneficial species.

Insect migration occurs on an enormous scale (1), with billions of individuals undertaking multigenerational migrations between seasonally favorable climatic zones around the globe (2–6). These long-range migration cycles profoundly influence terrestrial ecosystems via the large-scale transfer of biomass, energy, and nutrients (4–8), the provision of ecosystem services (8–10), impacts on agricultural productivity (11), and spread of disease (5, 12); thus, it is imperative that we better understand insect movement patterns. Recently, there has been a step change in our knowledge of the year-round spatial distribution and migratory routes of a few well-studied species (13, 14), particularly the monarch butterfly (Danaus plexippus) (15, 16) and (to a lesser extent) the painted lady butterfly (Vanessa cardui) (17, 18). However, interannual population dynamics of such insect migrants remain poorly known. One of the characteristic features is the interannual variation in the abundance of the first wave of immigrants to reach the temperate zone, which can vary by several orders of magnitude between successive years (3, 11, 17, 19). It is generally believed that this variation is driven by the effect of winter climate on breeding success and survival in the tropical and subtropical winter-breeding regions, particularly when these regions are arid or semiarid (1, 3, 20).Here, we study the painted lady butterfly, a cosmopolitan, continuously breeding migrant that undertakes seasonally predictable, long-range movements between tropical/subtropical winter-breeding regions and temperate zone summer-breeding regions (17–20). We focus on the western portion of its Afro-Palearctic migration system (from the Gulf of Guinea to Fennoscandia) due to the unparalleled monitoring data available on the spring and summer generations in parts of this range (Butterfly Monitoring Scheme [BMS] data from western Europe; Fig. 1 A and B) and in order to quantify the environmental drivers of interannual variation in abundance. In this western section, winter breeding was traditionally considered to occur predominantly in the Maghreb region of northwestern (NW) Africa (21, 22). However, recent studies suggest that it can occur over a much larger latitudinal range, from the Gulf of Guinea coast to the north Mediterranean coast, with two-way movement across the Sahel and Sahara Desert linking European and sub-Saharan African populations (23–26). The colonization of Western Europe consists of a northward progression of successive generations throughout spring and summer. The European component starts when butterflies that had emerged in the Maghreb (Morocco, Algeria, and Tunisia) a few days previously (17, 27) arrive in the Mediterranean region during March and April and immediately produce the next generation there. What is unclear is just how important the winter generations produced south of the Sahara are in seeding or reinforcing the early-spring generation in the Maghreb. Subsequent late-spring and summer generations reach as far north as Fennoscandia, and then the autumn generation undertakes an extremely long migration [often high above the ground, utilizing fast tailwinds (17, 28)] back to NW Africa and sub-Saharan West Africa (17, 29). Here, the annual cycle, comprising six or more generations per year, resumes (17–20, 24).Open in a separate windowFig. 1.Painted lady population data in western Europe. (A) Phenology of painted ladies in Europe showing peaks that correspond to either migrants or local generations. In the Mediterranean region (NE Spain), the light-blue period corresponds to the spring immigration and the dark-blue period to the summer emergence of a locally bred generation. In NW Europe (NL: the Netherlands; Eng & Wal: England and Wales), the light-pink period corresponds to the early-summer immigration and the dark-pink period to the late-summer emergence of a locally bred generation. (B) Log-collated annual index (across all sites in each country) for NE Spain in spring (1 March to 30 May) and summer (1 June to 31 July) and for NW Europe in early summer (15 May to 15 July) and late summer (16 July to 30 September). Abundance indices are expressed on a log scale, with zero reflecting the average for that region and season across all years. See Fig. 4 for the factors explaining years of peak abundance (e.g., 1996, 2003, 2006, 2009, and 2015).Extreme interannual variation in the abundance of the spring immigrants (and subsequent summer population) is a feature of painted lady population dynamics in both Europe (17, 19, 20, 30) and North America (29, 31, 32). Some painted ladies arrive in western Europe every spring; however, the pattern of abundance is one of irregular spectacular mass arrivals interspersed with years of much-reduced immigration (Fig. 1B). Here, we determine the key environmental conditions, and when/where they act during the migratory cycle, that drive this extreme annual variability in the European population dynamics each summer. In particular, we tackle the question of whether sub-Saharan, North African, and/or southern Iberian environmental conditions during the previous winter or spring are the primary drivers of the size of the spring and early-summer immigrations, initially to the Mediterranean and ultimately to northern Europe. To identify links between the generations monitored in Europe and the African breeding cycles both north and south of the Sahara, we use the following: 1) winter and spring environmental data (normalized difference vegetation index [NDVI], precipitation, temperature, and frequency of favorable tailwinds) covering the critical regions and periods (Fig. 2) in which large populations could potentially originate; 2) 21 y of BMS records from the Mediterranean (northeastern [NE] Spain) and NW Europe (the United Kingdom and the Netherlands); and 3) atmospheric trajectory simulations along the migratory route from sub-Saharan Africa to Europe.Open in a separate windowFig. 2.Correlations between spring painted lady counts in NE Spain with the NDVI, precipitation, and temperature. Red areas on the maps indicate regions that have positive significant correlations between the variable plotted and spring painted lady counts in NE Spain, while blue areas are negative correlations. See also Fig. 3 and SI Appendix, Fig. S2 for plots of painted lady spring numbers against the winter NDVI. These correlation plots were used to identify the ecoregions that were likely to be important (see delineation of these ecoregions in Fig. 3 and SI Appendix, Fig. S1) and to select the most important variables for the modeling (SI Appendix, Tables S1 and S2).     

Iceberg discharges of the last glacial period driven by oceanic circulation changes

Proxy data reveal the existence of episodes of increased deposition of ice-rafted detritus in the North Atlantic Ocean during the last glacial period interpreted as massive iceberg discharges from the Laurentide Ice Sheet. Although these have long been attributed to self-sustained ice sheet oscillations, growing evidence of the crucial role that the ocean plays both for past and future behavior of the cryosphere suggests a climatic control of these ice surges. Here, we present simulations of the last glacial period carried out with a hybrid ice sheet–ice shelf model forced by an oceanic warming index derived from proxy data that accounts for the impact of past ocean circulation changes on ocean temperatures. The model generates a time series of iceberg discharge that closely agrees with ice-rafted debris records over the past 80 ka, indicating that oceanic circulation variations were responsible for the enigmatic ice purges of the last ice age.Compared with the present interglacial period, the last glacial period (LGP) (∼110–10 ka before the present), and almost certainly previous ones (1), were characterized by substantial climatic variability on millennial timescales. This variability is mainly manifested in two types of events. Dansgaard–Oeschger (D/O) events are most notable in Greenland ice core records and involve decadal-scale warming of more than 10 K (interstadials) followed by slow cooling lasting several centuries and a final more rapid fall to cold background (stadial) conditions (2). Heinrich (H) events consist of massive iceberg discharges from the Laurentide Ice Sheet at intervals of ∼7 ka during peak glacial conditions throughout the LGP (3). Both D/O and H events are associated with widespread centennial- to millennial-scale climatic changes, including a synchronous temperature response over the North Atlantic and an antiphase temperature relationship over Antarctica and most of the Southern Ocean, as revealed by a wealth of deep-sea sediments, ice core, and terrestrial records (4). The Atlantic meridional overturning circulation (AMOC) is thought to play a central role in these abrupt glacial climatic changes. Although the paleoceanographic evidence on this link is scarce and mostly restricted to a few high-resolution deep-sea sediment records of the last deglaciation (5, 6), both modeling studies and reconstructions provide strong support for the hypothesis that D/O events were caused by reorganizations of the AMOC (7, 8). H events, identified as enhanced ice-rafted detritus (IRD) in North Atlantic deep-sea sediments (3, 9), occur during climatic minima of the Northern Hemisphere. They have classically been attributed to internal oscillations of the Laurentide (10) and assumed to lead to important disruptions of the Atlantic Ocean circulation (11). However, paleoclimate data have revealed that most H events likely occurred about a thousand years after North Atlantic Deep Water (NADW) formation had already slowed down or largely collapsed (12, 13), implying that the initial AMOC reduction could not have been caused by the H events themselves. This evidence directly conflicts with the common interpretation that freshwater fluxes representing the iceberg discharges caused the shift into cold (i.e., stadial) conditions. This furthermore highlights the need for a new paradigm through which to understand the triggering mechanism of H events. As already advanced one decade ago (14), any new theory should be able to account for the fact that the cold periods in which H events appear are not caused by the iceberg discharges and that the latter occur systematically several centuries after the North Atlantic cooling. More recently, the interaction between ocean circulation and ice sheet dynamics has been suggested to play a major role in triggering H events (15–17). This hypothesis has been assessed in particular for the first H event, H1, with both models and data showing that reduced NADW formation and a weakened AMOC lead to subsurface warming in the Nordic and Labrador Seas. This results in rapid melting of the Labrador ice shelves causing substantial ice stream acceleration and enhanced iceberg discharge (18–20).Here, we investigate the effects of oceanic circulation changes associated with millennial-scale climate variability on the Laurentide Ice Sheet dynamics within a more realistic modeling framework. To this end, we drive a hybrid ice sheet–ice shelf model (21) with time-varying oceanic subsurface temperature fields for the LGP (Materials and Methods) obtained by combining glacial climate simulations and information from proxy data. Climatic boundary conditions are otherwise fixed to glacial conditions, so that the only external forcing felt by the ice sheet model is the change in subsurface ocean temperatures. These are translated into basal melting rates via a linear equation dependent on a single tunable parameter (see SI Text for details and sensitivity tests). Climate simulations are performed with a global atmosphere–ocean model for glacial stadial and interstadial states (with weak and strong AMOC states, respectively) (22). These provide the range and spatial distribution of oceanic temperatures felt by the ice sheet. The temporal millennial-scale variability is based on a proxy-derived index used to interpolate in time between the stadial and interstadial ocean temperature fields. To produce this index, we assume that millennial-scale variability registered in the Greenland Ice Core Project (GRIP) ice core record (2) reflects variations in the North Atlantic oceanic state (Fig. 1A). To characterize the latter, following previous work (1), we use a threshold in the derivative of the GRIP signal to determine the timing of stadial to interstadial transitions (Fig. 1B). This allows for an objective classification of climatic states into stadials or interstadials (i.e., cold and warm surface periods). We furthermore assume millennial-scale variability as registered in the GRIP record reflects variations in NADW formation that have an imprint on subsurface temperatures in antiphase with respect to the surface state. Stadials are thus associated with periods of reduced NADW formation and weakened AMOC and warm subsurface temperatures, whereas during interstadials a stronger AMOC with active NADW formation cools the subsurface, in agreement with previous studies (15, 16, 19, 20, 23, 24). Considering a fast relaxation time of the subsurface temperature when convection resumes, and slow relaxation when convection is weak (20), we generate a subsurface warming index that slowly peaks during stadial climatic periods and more abruptly collapses when entering interstadial climatic periods (Fig. 1C). This index is thus directly derived from the GRIP time series and represents the only external forcing to the ice sheet model (see Materials and Methods and SI Text for details).Open in a separate windowFig. 1.Derivation of the Labrador Sea subsurface oceanic index. (A) GRIP d18O ice core record. (B) Smooth derivative of GRIP d18O record (red) with positive and negative thresholds which define transitions between stadial and interstadial states (grey). (C) Subsurface warming index. The cold subsurface state corresponds to an interstadial state (i.e., warm surface “climatic” state) with a mean subsurface (700- to 1,100-m depth) temperature of −0.9 °C. The warm subsurface state corresponds to a stadial state (i.e., cold surface “climatic” state) with a mean subsurface temperature of 1.1 °C.The application of subsurface oceanic forcing to the ice sheet model induces significant millennial-scale variability in the otherwise stable Laurentide Ice Sheet, as reflected in the velocity at the Hudson Strait outlet and iceberg discharge into the ocean (Fig. 2). For almost every peak in subsurface warming, there is a corresponding large and abrupt acceleration of the ice flow. Transitions between slow and fast states of the Hudson Strait ice stream occur several times during the LGP, with velocities varying between ∼1,000 m⋅a⁻¹ during buttressing periods and ∼4,000 m⋅a⁻¹ during periods of ice shelf breakup. The magnitude of the velocity does not directly correlate with the magnitude and duration of the subsurface warming, however, because of the competing timescales of ice sheet growth, ice advection from inland, and ice shelf breakup and growth. These three mechanisms lead to a nonlinear response of the system that appears to modulate the dynamics of the floating and inland ice in this region. When interstadial subsurface (i.e., cold) oceanic conditions are applied, the Labrador Sea ice shelf experiences low melt rates and can extend far enough to reach the western coast of Greenland (Fig. 3). In this way, significant backforce is felt by the Hudson Strait ice stream and velocities are greatly reduced. This allows the main Laurentide ice dome to grow and subsequently advect ice from inland toward the margin because of a permanently active Hudson Strait ice stream, preconditioning the ice sheet for more ice discharge into the ocean. When stadial subsurface (i.e., warm) oceanic conditions are applied, the ice shelf melts away from Greenland and no longer buttresses the ice stream that feeds it (Fig. 3). This allows a surge of velocity at the mouth of the ice stream, which propagates inland over several centuries and results in a significant increase in ice discharge into the Labrador Sea (Fig. 4). The magnitude of such a discharge event depends on the state of the ice sheet before the ice shelf collapse.Open in a separate windowFig. 2.(A) Labrador Sea subsurface oceanic index; (B) simulated Hudson Strait ice velocity (in kilometers per year); (C) simulated Labrador Sea calving rate (in Sverdrups); (D) magnetic susceptibility from core MD95-2024 (45.7°W, 50.2°N) (25); (E) lithic fraction from core JPC-13 (33.5°W, 53.1°N) (26). For the comparison, the timescales of the proxy data were converted to the SS09 timescale of the GRIP record (2).Open in a separate windowFig. 3.Laurentide ice stream velocities (in kilometers per year) before (Left) and during (Right) H event 2, along with locations of the cores of the IRD proxies shown in Fig. 2. The dashed line in the right panel indicates the location of the profiles shown in Fig. 4.Open in a separate windowFig. 4.Ice sheet profiles of the Laurentide (as indicated in Fig. 3) before and during HE2. Time series show ice-shelf thickness (gray; in meters), basal stress (10⁵ Pa), velocity (in kilometers per year), and thickness (in meters) for the upstream (magenta) and downstream (dark blue) sections of the Hudson Strait ice stream. The background shading in the right panel represents buttressed (light blue), transition (light red) and unbuttressed (white) periods.The simulated time series of calving into the Labrador Sea compares very well with proxies of calving obtained from marine sediment cores from the North Atlantic (Figs. 2 and ​and3).3). Both a high-resolution record of magnetic susceptibility from core MD95-2024 (45.7°W, 50.2°N) (25) and a record of lithic fraction from core JPC-13 (33.5°W, 53.1°N) (26), i.e., IRD proxies, show the same timing of peaks corresponding to major discharge events. In some isolated cases, such as between H4 and H3, or between H2 and H1, the simulated time series agrees better with the latter core. However, no spurious discharge events are simulated that are not apparent in at least one core. A comparison of the time series of the prescribed subsurface warming, calving, and IRD proxies highlights the fact that, for every strong peak in calving (i.e., H event), there is a necessary peak in subsurface warming. Several sensitivity tests with Grenoble ice shelf and land ice model (GRISLI) show that the amount of calving produced is a result of the nonlinear preconditioning of the ice sheet (SI Text). However, in our simulations the triggering mechanism for large ice discharges is always an ice shelf breakup precipitated by subsurface warming in the Labrador Sea. The strongest calving events are furthermore found to take place during the longest subsurface warming periods. We would expect the same relationship if subsurface temperature reconstructions were available, reflecting positive feedbacks operating between NADW formation and ice discharge (16). A more persistent reduction in NADW formation results in longer periods of subsurface warming. This in turn has a larger impact on ice shelves, which tends to increase ice discharge and suppress NADW formation further.H events are among the most dramatic examples of millennial-scale variability of the Quaternary climate and their interpretation has remained elusive for decades. In recent years, the increasing availability of observations of the present-day ice sheets has confirmed the unexpected and crucial role that the ocean exerts on the dynamics of the ice sheets (27). Ice shelves represent the necessary interface for this coupled system. Whereas little information exists for a constrained reconstruction of the floating parts of the Laurentide, the maximum extent of the ice shelf simulated here is restricted to the continental shelf area between Greenland and the Hudson Strait. Such a configuration does not appear to contradict the relatively sparse proxy data available in this region (28), and is glaciologically consistent. Furthermore, the existence of a persistent ice stream through the Hudson Strait, as simulated here, is supported by geological evidence and modeling (29, 30).Combined with the simulations presented here, the fact that the subsurface warming index generated from GRIP data aligns so well with the IRD proxies lends strong support to the hypothesis that millennial-scale glacial ice discharges are the result of a response to oceanic forcing. A characteristic time longer than the forcing timescale is the result of the nonlinearities of the ice sheet/ice shelf system. These arise from the different characteristic times of the ice shelf breakup and regrowth and by the time needed by the ice sheet to propagate the signal from its oceanic perturbation across the ice streams (Fig. 4). These phenomena favor the occurrence of resonance in the system and finally determine the observed pacing of ∼7 ka.Our simulations provide a physically based framework through which to understand the coupled ice sheet–ocean system. Open questions remain concerning the relationship between IRD proxies and actual calving rates, which can result from outburst floods, iceberg melting, and ocean circulation changes (31). One important related aspect concerns the fact that it is very difficult to constrain the melting rates that icebergs experience during their trip across the North Atlantic. This allows for alternative explanations considering the observed IRD belt as mainly the reflection of colder oceanic temperatures when Heinrich layers were formed (32). Under this interpretation, however, the amount of IRDs in marine cores close to the ice sheet source would reflect a signal absent of Heinrich-like events. This seems not to be the case, because Heinrich peaks can be observed in cores of the Labrador Sea (33). However, the explanation for the ultimate causes behind the underlying glacial oceanic variability remains elusive. Nonetheless, the work presented here shows that proxies and modeling reveal a consistent picture of the origin of the massive iceberg discharges of the last glacial cycle, including the enigmatic H events.     

Conjugated polymers based on metalla-aromatic building blocks

Conjugated polymers usually require strategies to expand the range of wavelengths absorbed and increase solubility. Developing effective strategies to enhance both properties remains challenging. Herein, we report syntheses of conjugated polymers based on a family of metalla-aromatic building blocks via a polymerization method involving consecutive carbyne shuttling processes. The involvement of metal d orbitals in aromatic systems efficiently reduces band gaps and enriches the electron transition pathways of the chromogenic repeat unit. These enable metalla-aromatic conjugated polymers to exhibit broad and strong ultraviolet–visible (UV–Vis) absorption bands. Bulky ligands on the metal suppress π–π stacking of polymer chains and thus increase solubility. These conjugated polymers show robust stability toward light, heat, water, and air. Kinetic studies using NMR experiments and UV–Vis spectroscopy, coupled with the isolation of well-defined model oligomers, revealed the polymerization mechanism.

Conjugated polymers are macromolecules usually featuring a backbone chain with alternating double and single bonds (1–3). These characteristics allow the overlapping p-orbitals to form a system with highly delocalized π-electrons, thereby giving rise to intriguing chemical and physical properties (4–6). They have exhibited many applications in organic light-emitting diodes, organic thin film transistors, organic photovoltaic cells, chemical sensors, bioimaging and therapies, photocatalysis, and other technologies (7–10). To facilitate the use of solar energy, tremendous efforts have been devoted in recent decades to developing previously unidentified conjugated polymers exhibiting broad and strong absorption bands (11–13). The common strategies for increasing absorption involve extending π-conjugation by incorporating conjugated cyclic moieties, especially fused rings; modulating the strength of intramolecular charge transfer between donor and acceptor units (D–A effect); increasing the coplanarity of π conjugation through weak intramolecular interactions (e.g., hydrogen bonds); and introducing heteroatoms or heavy atoms into the repeat units of conjugated polymers (11–16). Additionally, appropriate solubility is a prerequisite for processing and using polymers and is usually achieved with the aid of long alkyl or alkoxy side chains (12, 17).Aromatic rings are among the most important building blocks for conjugated polymers. In addition to aromatic hydrocarbons, a variety of aromatic heterocycles composed of main-group elements have been used as fundamental components. These heteroatom-containing conjugated polymers show unique optical and electronic properties (4–10). However, while metalla-aromatic systems bearing a transition metal have been known since 1979 due to the pioneering work by Thorn and Hoffmann (18), none of them have been used as building blocks for conjugated polymers. The HOMO–LUMO gaps (E_g) of metalla-aromatics are generally narrower (Fig. 1) than those of their organic counterparts (19–22). We reasoned that this feature should broaden the absorption window if polymers stemming from metalla-aromatics are achievable.Open in a separate windowFig. 1.Comparison of traditional organic skeletons with metalla-aromatic building blocks (the computed energies are in eV). (A) HOMO–LUMO gaps of classic aromatic skeletons. (B) Carbolong frameworks as potential building blocks for novel conjugated polymers with broad absorption bands and improved solubility.In recent years, we have reported a series of readily accessible metal-bridged bicyclic/polycyclic aromatics, namely carbolong complexes, which are stable in air and moisture (23–25). The addition of osmium carbynes (in carbolong complexes) and alkynes gave rise to an intriguing family of d_π–p_π conjugated systems, which function as excellent electron transport layer materials in organic solar cells (26, 27). These observations raised the following question: Can this efficient addition reaction be used to access metalla-aromatic conjugated polymers? It is noteworthy that incorporation of metalla-aromatic units into conjugated polymers is hitherto unknown. In this contribution, we disclose a polymerization reaction involving M≡C analogs of C≡C bonds, which involves a unique carbyne shuttling strategy (Fig. 2A). This led to examples of metalla-aromatic conjugated polymers (polycarbolongs) featuring metal carbyne units in the main chain. On the other hand, the development of polymerization reactions plays a crucial role in involving certain building blocks in conjugated polymers (28–32). These efficient, specific, and feasible polymerizations could open an avenue for the synthesis of conjugated polymers.Open in a separate windowFig. 2.Design of polymers and synthesis of monomers. (A) Schematic illustration of the polymerization strategy. (B) Preparation of carbolong monomers. Insert: X-ray molecular structure for the cations of complex 3. Ellipsoids are shown at the 50% probability level; phenyl groups in PPh₃ are omitted for clarity.     

Universal features in the photoemission spectroscopy of high-temperature superconductors

The energy gap for electronic excitations is one of the most important characteristics of the superconducting state, as it directly reflects the pairing of electrons. In the copper–oxide high-temperature superconductors (HTSCs), a strongly anisotropic energy gap, which vanishes along high-symmetry directions, is a clear manifestation of the d-wave symmetry of the pairing. There is, however, a dramatic change in the form of the gap anisotropy with reduced carrier concentration (underdoping). Although the vanishing of the gap along the diagonal to the square Cu–O bond directions is robust, the doping dependence of the large gap along the Cu–O directions suggests that its origin might be different from pairing. It is thus tempting to associate the large gap with a second-order parameter distinct from superconductivity. We use angle-resolved photoemission spectroscopy to show that the two-gap behavior and the destruction of well-defined electronic excitations are not universal features of HTSCs, and depend sensitively on how the underdoped materials are prepared. Depending on cation substitution, underdoped samples either show two-gap behavior or not. In contrast, many other characteristics of HTSCs, such as the dome-like dependence of on doping, long-lived excitations along the diagonals to the Cu–O bonds, and an energy gap at the Brillouin zone boundary that decreases monotonically with doping while persisting above (the pseudogap), are present in all samples, irrespective of whether they exhibit two-gap behavior or not. Our results imply that universal aspects of high- superconductivity are relatively insensitive to differences in the electronic states along the Cu–O bond directions.Elucidating the mechanism of high-temperature superconductivity in the copper–oxide materials remains one of the most challenging open problems in physics. It has attracted the attention of scientists working in fields as diverse as materials science, condensed matter physics, cold atoms, and string theory. To clearly define the problem of high-temperature superconductors (HTSCs), it is essential to establish which of the plethora of observed features are universal, namely, qualitatively unaffected by material-specific details.An important early result concerns the universality of the symmetry of the order parameter for superconductivity. The order parameter was found to change sign under a 90° rotation (1, 2), which implies that the energy gap must vanish along the diagonal to the Cu–O bonds, i.e., the Brillouin zone diagonal. This sign change is consistent with early spectroscopic studies of near-optimally-doped samples (those with the highest in a given family), where a energy gap (3, 4) was observed (ϕ being the angle from the Cu–O bond direction), the simplest functional form consistent with d-wave pairing. More recently, there is considerable evidence (5–8) that, with underdoping, the anisotropy of the energy gap deviates markedly from the simple form. Although the gap node at is observed at all dopings, the gap near the antinode (near and 90°) is significantly larger than that expected from the simplest d-wave form. Further, the large gap continues to persist in underdoped (UD) materials as the normal-state pseudogap (9–11) above . This suggests that the small (near-nodal) and large (antinodal) gaps are of completely different origin, the former related to superconductivity and the latter to some other competing order parameter.This two-gap picture has attracted much attention (8), raising the possibility that multiple energy scales are involved in the HTSC problem. There is mounting evidence for additional broken symmetries (12–14) in UD cuprates, once superconductivity is weakened upon approaching the Mott insulating state. The central issue is the role of these additional order parameters in impacting the universal properties of high- superconductivity.In this paper we use angle-resolved photoemission (ARPES) to examine the universality of the two-gap scenario in HTSCs by addressing the following questions. To what extent are the observed deviations from a simple d-wave energy gap independent of material details? How does the observed gap anisotropy correlate, as a function of doping, with other spectroscopic features such as the size of the antinodal gap, and the spectral weights of the nodal and antinodal quasiparticle excitations?We systematically examine the electronic spectra of various families of cation-substituted Bi₂Sr₂CaCu₂O_8+δ single crystals as a function of carrier concentration to elucidate which properties are universal and which are not. We present ARPES data on four families of float-zone-grown Bi₂Sr₂CaCu₂O_8+δ single crystals, where was adjusted by both oxygen content and cation doping. As-grown samples, labeled Bi2212, have an optimal of 91 K. These crystals were UD to by varying the oxygen content. Ca-rich crystals (grown from material with a starting composition Bi_2.1Sr_1.4Ca_1.5Cu₂O_8+δ) with an optimal of 82 K are labeled Ca. Two Dy-doped families grown with starting compositions Bi_2.1Sr_1.9Ca_{1 − x}Dy_xCu₂O_8+δ with x = 0.1 and 0.3 are labeled Dy1 and Dy2, respectively. A full list of the samples used and their determined from magnetization measurements are shown in SI Text, where we also show high-resolution X-ray data that give evidence for the excellent structural quality of our samples.Our main result is that the Dy1 and Dy2 samples show clear evidence of a two-gap behavior in the UD regime , with loss of coherent quasiparticles in the antinodal region of k space where the gap deviates from a simple d-wave form. In marked contrast, the UD Bi2212 samples and the Ca samples show a simple d-wave gap in the superconducting state and sharp quasiparticles over the entire Fermi surface in a similar range of the UD regime. We conclude by discussing the implications of the nonuniversality of the two-gap behavior for the phenomenon of high superconductivity.We begin our comparison of the various families of samples by focusing in Fig. 1 on the superconducting state antinodal spectra as a function of underdoping. The antinode is the Fermi momentum k_F on the Brillouin zone boundary, where the energy gap is a maximum and, as we shall see, the differences between the various samples are the most striking. We show data at optimal doping, corresponding to the highest in each family, in Fig. 1A. Increasing Dy leads to a small suppression of the optimal compared with Bi2212, together with an increase in the antinodal gap and a significant reduction of the quasiparticle weight. This trend continues down to moderate underdoping, as seen in Fig. 1B, where we show UD Bi2212 and Dy2 samples with very similar . For more severely UD samples, with , spectral changes in the Dy-substituted samples are far more dramatic. In Fig. 1C, we see that quasiparticle peaks in the Dy samples are no longer visible, even well below , consistent with earlier work on Y-doped Bi2212 and also Bi2201 and La_1.85Sr_0.15CuO₄ (5, 15–18). In contrast, Bi2212 and Ca-doped samples with comparable continue to exhibit quasiparticle peaks. In this respect the latter two are similar to epitaxially grown thin-film samples that exhibit quasiparticle peaks all of the way down to the lowest (19).Open in a separate windowFig. 1.Superconducting state antinodal ARPES spectra. We use the label “Bi2212” for samples without cation doping, “Dy1” for 10% Dy, “Dy2” for 30% Dy, and “Ca” for Ca-doped samples. The temperature is indicated along with . OP denotes optimal doped, UD underdoped, and OD overdoped samples. (A) Antinodal spectra for OP samples of three different families: Bi2212 (blue), Dy1 (green), and Dy2 (red), showing an increase in gap and a decrease in quasiparticle weight with increasing Dy content. (B) Antinodal spectra for UD samples with similar (≃66 K) for Bi2212 (blue) and Dy2 (red). As in A, there is a larger gap and smaller coherent weight in the Dy-substituted sample. (C) Same as in B, but for four UD samples with near 55 K for Bi2212 (dark blue), Ca (light blue), Dy1 (green), and Dy2 (red). The Bi2212 and Ca spectra are very similar to each other and quite different from those of the Dy1 and Dy2 materials. (D) Doping evolution of the antinodal spectra of four Dy1 samples from OP to UD . (E) Doping evolution of the antinodal spectra of four Dy2 samples from OP to UD . We see in D and E the sudden loss of quasiparticle weight for below 60 K. (F) Doping evolution of the antinodal spectra of three Bi2212 samples and three Ca samples, showing well-defined quasiparticle peaks in all cases.A significant feature of the highly UD Dy samples in Fig. 1C is that, in addition to the strong suppression of the quasiparticle peak, there is severe loss of low-energy spectral weight. To clearly highlight this, we show the doping evolution of antinodal spectra for the Dy1 (Fig. 1D) and Dy2 (Fig. 1E) samples. These observations are in striking contrast with the Bi2212 and Ca-doped data in Fig. 1F, where we do see a systematic reduction of the quasiparticle peak with underdoping, but not a complete wipeout of the low-energy spectral weight. To the extent that the superconducting state peak–dip–hump line shape (20, 21) originates from one broad normal-state spectral peak, the changes in spectra of the Dy materials are not simply due to a loss of coherence, but more likely a loss of the entire spectral weight near the chemical potential.The doping evolution of the k-dependent gap is illustrated in Figs. 2 and​and 3. 3. In Fig. 2 we contrast the optimally doped Dy1 (T_c = 86 K) sample (Fig. 2 A and B) with a severely UD Dy1 (T_c = 38 K) sample (Fig. 2 C and D ), the spectra being particle–hole-symmetrized to better illustrate the gap. The OP 86 K sample shows a well-defined quasiparticle peak over the entire Fermi surface (Fig. 2A) with a simple d-wave gap of the form (blue curve in Fig. 2B). For the UD 38 K sample, we see in Fig. 2C well-defined quasiparticles near the node (red spectra), but not near the antinode (blue spectra). The near-nodal gaps (red triangles in Fig. 2D) are obtained from the energy of quasiparticle peaks and continue to follow a d-wave gap (blue curve in Fig. 2D). However, once the quasiparticle peak is lost closer to the antinode, one has to use some other definition of the gap scale. We identify a break in the slope of the spectrum, by locating the energy scale at which it deviates from the black straight lines (Fig. 2C), which leads to the gap estimates (blue squares) in Fig. 2D.Open in a separate windowFig. 2.Superconducting state spectra and energy gap for OD and highly UD Dy1 samples. (A) Symmetrized spectra at k_F, from the antinode (Upper) to the node (Lower) for an OP 86 K Dy1 sample. (B) Gap as a function of Fermi surface angle (0° is the antinode and 45° the node). The blue curve is a d-wave fit to the data. (C) Same as A for an UD 38 K Dy1 sample. Curves, near the node, with discernible quasiparticle peaks are shown in red; those near the antinode are shown in blue. (D) Gap along the Fermi surface from data of C.Open in a separate windowFig. 3.Energy gap anisotropies of various samples. (A) OD 79 K Ca (where ); (B) UD 54 K Ca; (C) OP 81 K Dy2; and (D) UD 59 K Dy2. The two near-optimal samples in A and C both show a simple d-wave gap. This behavior persists in the UD Ca sample of B, but the UD Dy2 sample of D has a two-gap behavior despite having a similar to the UD Ca sample.Despite the larger error bar associated with gap scale extraction in the absence of quasiparticles, it is nevertheless clear (Fig. 2D) that the UD 38 K Dy1 sample has an energy gap that deviates markedly from the simple d-wave form. This observation is called two-gap in the UD regime, in contrast with a single gap near optimality (Fig. 2B). It is easy to observe from Fig. 2 that the Fermi surface angle at which the energy gap starts to deviate from the form matches the one at which the spectral peak gets washed out. This is very similar to the two-gap behavior demonstrated in refs. 5, 15–18. From this, one might conclude that two-gap behavior is directly correlated with a loss of well-defined quasiparticle excitations in the antinodal region. However, we point to recent ARPES data on Y-doped Bi2212 (6, 7), where two-gap behavior has been observed despite the presence of small antinodal quasiparticle peaks.We next show that the two-gap behavior is not a universal feature of all UD samples. To make this point, we compare in Fig. 3 the gap anisotropies of the Ca-doped samples (Fig. 3 A and B) with the Dy2 samples (Fig. 3 C and D) with essentially identical , where both families have the same optimal . The near-optimal samples, OD 79 K Ca (Fig. 3A) and OP 81 K Dy2 (Fig. 3C) samples, both have a simple d-wave anisotropy (although different maximum gap values at the antinode). However, upon underdoping to similar values, the two have markedly different gap anisotropies. The UD 59 K Dy2 sample (Fig. 3D) shows two-gap behavior, and an absence of quasiparticles near the antinode (similar to the discussion in connection with Fig. 2 above). However, the UD 54 K Ca sample (Fig. 3B) continues to exhibit sharp spectral peaks and a single-gap, despite a very similar as the UD 59 K Dy2.Having established the qualitative differences in the gap anisotropies for various samples as a function of underdoping, we next summarize in Fig. 4 the doping evolution of various spectroscopic features. Instead of estimating the carrier concentration in our samples using an empirical equation (22) (that may or may not be valid for various cation substitutions), we prefer to use the measured to label the doping. In Fig. 4A we show the doping evolution of the antinodal energy gap, which is consistent with the known increase in the gap with underdoping.Open in a separate windowFig. 4.Antinodal gaps and quasiparticle weights. (A) Antinodal energy gap as a function of doping for various samples is seen to grow monotonically with underdoping. Here, and in B and C, the doping is characterized by the measured quantity , with UD samples shown to the left of and OD samples to the right. All results are at temperatures well below . (B) Coherent spectral weight for antinodal quasiparticles as a function of doping. Dy-doped samples exhibit a rapid suppression of this weight to zero for UD , whereas the Ca-doped samples show robust antinodal peaks even for . (C) Coherent spectral weight for nodal quasiparticles as a function of doping, which is seen to be much more robust than the antinodal one.The coherent spectral weight Z for antinodal quasiparticles is plotted in Fig. 4B (for details on the procedure used to estimate this weight, from a ratio of spectral areas, see SI Text). The Dy1 and Dy2 samples both show a sudden and complete loss of Z with underdoping (23), which coincides with the appearance of two-gap behavior. In marked contrast with the Dy samples, the Bi2212 and Ca samples that exhibit a single d-wave gap show a gradual drop in the antinodal Z. On the other hand, we find that the nodal excitations are much less sensitive to how the sample is UD compared with the antinodal ones. Similar sharp nodal excitations have been observed in Dy-doped Bi2212 samples in ref. 7 as well. The nodal quasiparticle weight Z in Fig. 4C decreases smoothly with underdoping for all families of samples, as expected for a doped Mott insulator (24).The two-gap behavior and the attendant loss of quasiparticle weight near the antinode imply a nodal–antinodal dichotomy, aspects of which have been recognized in k space (25–27) and in real space (28–30). Two possible, not mutually exclusive, causes of this behavior are disorder and competing orders.It is known that antinodal states are much more susceptible to impurity scattering, whereas near-nodal excitations are protected (31). However, it is not a priori clear why certain cation substitutions (Dy) should lead to more electronic disorder than others (Ca). As shown by our X-ray studies in SI Text, there is no difference in the structural disorder in Dy and Ca samples. One possibility is that Dy has a local moment, but there is no direct experimental evidence for this.The two-gap behavior in UD materials, with a large antinodal gap that persists above , is suggestive of an order parameter, distinct from d-wave superconductivity, which sets in at the pseudogap temperature . There are several experiments (12–14) that find evidence for a broken symmetry at . However, it is not understood how the observed small, and often subtle, order parameter(s) could lead to large antinodal gaps of , with a loss of spectral weight over a much larger energy range (Fig. 1 D and E).We now discuss the pertinence of competing order parameters based on our measurements. First, in our ARPES data, we have not found any direct evidence for density wave ordering (say, from zone folding). Second, our X-ray data did not provide any signature for additional diffraction peaks expected for long-range density wave ordering. However, none of these null results provide definitive evidence for the absence of a density wave ordering, particularly if it were short range. In contrast, in previously published work (5, 15–18), two-gap behavior has been conjectured to be a direct consequence of phase competition between d-wave superconductivity and some type of density wave ordering. As we have demonstrated, two-gap behavior in and of itself is a sample-specific issue and hence, even if we assume a linkage between competing order and two-gap behavior, it cannot be central to the question of superconductivity in HTSC systems.Whatever the mechanism leading to qualitatively different gap anisotropies for the UD Dy and Ca samples, it only produces relatively small, quantitative changes in key aspects of these materials, such as the dependence of on doping, the presence of sharp nodal quasiparticles, and the pseudogap. We thus conclude that antinodal states do not make a substantial contribution to the universal features of HTSCs. Clearly, two gaps are not necessary for high-temperature superconductivity.     

Mandibular and dental characteristics of Late Triassic mammaliaform Haramiyavia and their ramifications for basal mammal evolution

As one of the earliest-known mammaliaforms, Haramiyavia clemmenseni from the Rhaetic (Late Triassic) of East Greenland has held an important place in understanding the timing of the earliest radiation of the group. Reanalysis of the type specimen using high-resolution computed tomography (CT) has revealed new details, such as the presence of the dentary condyle of the mammalian jaw hinge and the postdentary trough for mandibular attachment of the middle ear—a transitional condition of the predecessors to crown Mammalia. Our tests of competing phylogenetic hypotheses with these new data show that Late Triassic haramiyids are a separate clade from multituberculate mammals and are excluded from the Mammalia. Consequently, hypotheses of a Late Triassic diversification of the Mammalia that depend on multituberculate affinities of haramiyidans are rejected. Scanning electron microscopy study of tooth-wear facets and kinematic functional simulation of occlusion with virtual 3D models from CT scans confirm that Haramiyavia had a major orthal occlusion with the tallest lingual cusp of the lower molars occluding into the lingual embrasure of the upper molars, followed by a short palinal movement along the cusp rows alternating between upper and lower molars. This movement differs from the minimal orthal but extensive palinal occlusal movement of multituberculate mammals, which previously were regarded as relatives of haramiyidans. The disparity of tooth morphology and the diversity of dental functions of haramiyids and their contemporary mammaliaforms suggest that dietary diversification is a major factor in the earliest mammaliaform evolution.Haramiyidans are among the first mammaliaforms to appear during the Late Triassic in the evolutionary transition from premammalian cynodonts. Their fossils have a cosmopolitan distribution during the Late Triassic to the Jurassic (1–8), tentatively with the youngest record in the Late Cretaceous of India (9). Most of these occurrences are of isolated teeth. For this reason, Haramiyavia clemmenseni (1) holds a special place in mammaliaform phylogeny: It is the best-preserved Late Triassic haramiyid with intact molars, nearly complete mandibles, and also postcranial skeletal elements (Figs. 1 and ​and22 and SI Appendix, Figs. S1–S4) (1). By its stratigraphic provenance from the Tait Bjerg Beds of the Fleming Fjord Formation, East Greenland (Norian-Rhaetic Age) (7), Haramiyavia is also the oldest known haramiyid (5, 7). Haramiyids, morganucodonts, and kuehneotheriids are the three earliest mammaliaform groups that are distinctive from each other in dental morphology and masticatory functions (10–12).Open in a separate windowFig. 1.(A and B) Composite reconstruction of Haramiyavia clemmenseni right mandible in lateral (A) and medial (B) views. Dark red: original bone with intact periosteal surface; brown: broken surface of preserved bone or remnant of bone; light blue: morphologies preserved in mold outlines or clear impression. (C) Morganucodon mandible in medial view.Open in a separate windowFig. 2.Molar features of Haramiyavia. (A) Right M1–M3 in occlusal view (medio-lateral orientation by the zygomatic root and the palate). (B) Occlusal facets of upper molars. (C) Lingual view of M1–M3. (D) Root structures of upper molars (M1 and M2 show three partially divided anterior roots connected by dentine and two posterior roots connected by dentine; M3 has two anterior and two posterior roots connected respectively by dentine). These roots have separate root canals. (E) Buccal view of M1–M3. All roots are bent posteriorly, suggesting that crowns shifted mesially, relative to the roots, during the tooth eruption, also known as mesial drift of teeth (arrowhead), typical of successive eruption of multirooted postcanines. (F) SEM photograph of lower m3 in a posterior occlusal view. (G) Approximate extent of wear facets by orthal occlusion (a1 cusp in embrasure of upper molars) (blue) and palinal movement of b2–b4 cusps sliding across the median furrow of upper molars (green). (H) Lingual view of m3. There are no wear facets on lingual side of cusps a1–a4. (I) Buccal view of m3 showing wear facets on the buccal sides of cusps b1–b4 and on apices.Haramiyids are characterized by their complex molars with longitudinal rows of multiple cusps. The cusp rows occlude alternately between the upper and lower molars. Primarily because of similarities in molar morphology, haramiyids are considered to be related to poorly known theroteinids of the Late Triassic (5, 13) and eleutherodontids of the Middle to Late Jurassic (14–17). Collectively haramiyids and eleutherodontids are referred to as “haramiyidans” (10, 14, 15, 18, 19). Recent discoveries of diverse eleutherodontids or eleutherodontid-related forms with skeletons from the Tiaojishan Formation (Middle to Late Jurassic) of China (18–20) have greatly augmented the fossil record of haramiyidans, ranking them among the most diverse mammaliaform clades of the Late Triassic and Jurassic.Historically, it has been a contentious issue whether haramiyidans (later expanded to include theroteinids and eleutherodontids) are closely related to the more derived multituberculates from the Middle Jurassic to Eocene (13) or represent a stem clade of mammaliaforms excluded from crown mammals (21, 22). The conflicting placement of haramiyidans was attributable in part to the uncertainties in interpreting the isolated teeth of most Late Triassic haramiyids (21, 22). More recent phylogenetic disagreements have resulted from different interpretations of mandibular characters in Haramiyavia (17–20, 23–25), which has not been fully described (figure 2 in ref. 1).Here we present a detailed study of the mandibles and teeth of Haramiyavia from the exhaustive documentation during initial fossil preparation (Fig. 1 and SI Appendix, Figs. S1–S4), from scanning electron microscopy (SEM) images, and from computed tomography (CT) scans and 3D image analyses of the two fossil slabs with mandibles (MCZ7/95A and B), plus a referred specimen of upper molars in a maxilla (MCZ10/G95) (Figs. 2 and ​and3,3, SI Appendix, Figs. S5–S8 and Tables S2 and S3, and Movie S1). These new data are informative for testing alternative mammaliaform phylogenies (Fig. 4 and SI Appendix) and are useful for reconstructing evolutionary patterns of feeding function in the earliest mammaliaforms.Open in a separate windowFig. 3.Molar occlusion of haramiyids. (A) In Haramiyavia the upper and lower molars form an en echelon pattern, a series of parallel and step-like occlusal surfaces in lingual and buccal views (based on 3D scaled models from CT scans of MCZ7/G95 and MCZ10/G95). (B–E) In Haramiyavia are shown the occlusal paths of cusps a1–a4 of the lingual row (B), cusps b2–b4 of the buccal row (shown with the lingual half of the tooth cut away) (C), and tooth orientation and the cut-away plane (D). During the orthal occlusion phase, the tallest lingual cusp a1 occludes into the embrasure of the preceding and the opposing upper molars (B and E), and the tallest buccal cusp b2 occludes into the upper furrow and behind the A1–B1 saddle of the upper molar (C). During the palinal occlusion phase, cusps b1–b4 of the buccal row slide posteriorly in the upper furrow, and in the upper row cusps B5–B1 slide in the lower furrow (lower molars with blue and green shading, superpositioned by flipped upper molar in clear outlines) (E). (F) Extent of wear on molars during the orthal phase (blue) and the palinal movement (green) produced by OFA simulation (Movie S1). (G–I) In Thomasia, reconstruction of upper and lower molar series on the basis of wear surfaces and tooth crown morphology (revised from refs. 2, 3, and 10). (G) The en echelon occlusal surfaces in lingual view. (H) Orthal occlusion (blue) is followed by palinal occlusal movement (green). (I) Occlusal wear facets of molars. Facets worn by orthal occlusion are shown in blue, and facets worn by palinal occlusion are shown in gray hatching. Cusp and facet designations are after refs 3 and 6.Open in a separate windowFig. 4.Hypotheses concerning the phylogenetic relationship of Haramiyavia and timing estimates of the basal diversification of crown mammals. (A) Haramiyavia is a close relative of multituberculates, both nested in the crown Mammalia. This hypothesis (haramiyidan node position 1) was based on a misinterpretation of a previous illustration of a fragment of the mandible (17, 18). (B) Haramiyavia is a stem mammaliaform, as determined by incorporating the features preserved on both mandibles into phylogenetic estimates (haramiyidan node position 2). (C) Placement of Haramiyavia and other haramiyidans among mammaliaforms according to this study. Many mandibular features were treated as unknown by studies favoring a Late Triassic diversification of mammals (18, 23). A more complete sampling of informative features revealed by this study now has overturned the previous placement. Clades: crown Mammalia (node a); Mammaliaformes (node b); haramiyidans (node 1 or 2, alternative positions); Eleutherodontida (node c). The rescored datasets and analyses are presented in SI Appendix.     

Little evidence for enhanced phenotypic evolution in early teleosts relative to their living fossil sister group

Since Darwin, biologists have been struck by the extraordinary diversity of teleost fishes, particularly in contrast to their closest “living fossil” holostean relatives. Hypothesized drivers of teleost success include innovations in jaw mechanics, reproductive biology and, particularly at present, genomic architecture, yet all scenarios presuppose enhanced phenotypic diversification in teleosts. We test this key assumption by quantifying evolutionary rate and capacity for innovation in size and shape for the first 160 million y (Permian–Early Cretaceous) of evolution in neopterygian fishes (the more extensive clade containing teleosts and holosteans). We find that early teleosts do not show enhanced phenotypic evolution relative to holosteans. Instead, holostean rates and innovation often match or can even exceed those of stem-, crown-, and total-group teleosts, belying the living fossil reputation of their extant representatives. In addition, we find some evidence for heterogeneity within the teleost lineage. Although stem teleosts excel at discovering new body shapes, early crown-group taxa commonly display higher rates of shape evolution. However, the latter reflects low rates of shape evolution in stem teleosts relative to all other neopterygian taxa, rather than an exceptional feature of early crown teleosts. These results complement those emerging from studies of both extant teleosts as a whole and their sublineages, which generally fail to detect an association between genome duplication and significant shifts in rates of lineage diversification.Numbering ∼29,000 species, teleost fishes account for half of modern vertebrate richness. In contrast, their holostean sister group, consisting of gars and the bowfin, represents a mere eight species restricted to the freshwaters of eastern North America (1). This stark contrast between teleosts and Darwin''s original “living fossils” (2) provides the basis for assertions of teleost evolutionary superiority that are central to textbook scenarios (3, 4). Classic explanations for teleost success include key innovations in feeding (3, 5) (e.g., protrusible jaws and pharyngeal jaws) and reproduction (6, 7). More recent work implicates the duplicate genomes of teleosts (8–10) as the driver of their prolific phenotypic diversification (8, 11–13), concordant with the more general hypothesis that increased morphological complexity and innovation is an expected consequence of genome duplication (14, 15).Most arguments for enhanced phenotypic evolution in teleosts have been asserted rather than demonstrated (8, 11, 12, 15, 16; but see ref. 17), and draw heavily on the snapshot of taxonomic and phenotypic imbalance apparent between living holosteans and teleosts. The fossil record challenges this neontological narrative by revealing the remarkable taxonomic richness and morphological diversity of extinct holosteans (Fig. 1) (18, 19) and highlights geological intervals when holostean taxonomic richness exceeded that of teleosts (20). This paleontological view has an extensive pedigree. Darwin (2) invoked a long interval of cryptic teleost evolution preceding the late Mesozoic diversification of the modern radiation, a view subsequently supported by the implicit (18) or explicit (19) association of Triassic–Jurassic species previously recognized as “holostean ganoids” with the base of teleost phylogeny. This perspective became enshrined in mid-20th century treatments of actinopterygian evolution, which recognized an early-mid Mesozoic phase dominated by holosteans sensu lato and a later interval, extending to the modern day, dominated by teleosts (4, 20, 21). Contemporary paleontological accounts echo the classic interpretation of modest teleost origins (22–24), despite a systematic framework that substantially revises the classifications upon which older scenarios were based (22–25). Identification of explosive lineage diversification in nested teleost subclades like otophysans and percomorphs, rather than across the group as a whole, provides some circumstantial neontological support for this narrative (26).Open in a separate windowFig. 1.Phenotypic variation in early crown neopterygians. (A) Total-group holosteans. (B) Stem-group teleosts. (C) Crown-group teleosts. Taxa illustrated to scale.In contrast to quantified taxonomic patterns (20, 23, 24, 27), phenotypic evolution in early neopterygians has only been discussed in qualitative terms. The implicit paleontological model of morphological conservatism among early teleosts contrasts with the observation that clades aligned with the teleost stem lineage include some of the most divergent early neopterygians in terms of both size and shape (Fig. 1) (see, for example, refs. 28 and 29). These discrepancies point to considerable ambiguity in initial patterns of phenotypic diversification that lead to a striking contrast in the vertebrate tree of life, and underpins one of the most successful radiations of backboned animals.Here we tackle this uncertainty by quantifying rates of phenotypic evolution and capacity for evolutionary innovation for the first 160 million y of the crown neopterygian radiation. This late Permian (Wuchiapingian, ca. 260 Ma) to Cretaceous (Albian, ca. 100 Ma) sampling interval permits incorporation of diverse fossil holosteans and stem teleosts alongside early diverging crown teleost taxa (Figs. 1 and ​and2A2A and Figs. S1 and ​andS2),S2), resulting in a dataset of 483 nominal species-level lineages roughly divided between the holostean and teleost total groups (Fig. 2B and Fig. S2). Although genera are widely used as the currency in paleobiological studies of fossil fishes (30; but see ref. 31), we sampled at the species level to circumvent problems associated with representing geological age and morphology for multiple congeneric lineages. We gathered size [both log-transformed standard length (SL) and centroid size (CS); results from both are highly comparable (Figs. S3 and ​andS4);S4); SL results are reported in the main text] and shape data (the first three morphospace axes arising from a geometric morphometric analysis) (Fig. 2A and Figs. S1) from species where possible. To place these data within a phylogenetic context, we assembled a supertree based on published hypotheses of relationships. We assigned branch durations to a collection of trees under two scenarios for the timescale of neopterygian diversification based on molecular clock and paleontological estimates. Together, these scenarios bracket a range of plausible evolutionary timelines for this radiation (Fig. 2B). We used the samples of trees in conjunction with our morphological datasets to test for contrasts in rates of, and capacity for, phenotypic change between different partitions of the neopterygian Tree of Life (crown-, total-, and stem-group teleosts, total-group holosteans, and neopterygians minus crown-group teleosts), and the sensitivity of these conclusions to uncertainty in both relationships and evolutionary timescale. Critically, these include comparisons of phenotypic evolution in early crown-group teleosts—those species that are known with certainty to possess duplicate genomes—with rates in taxa characterized largely (neopterygians minus crown teleosts) or exclusively (holosteans) by unduplicated genomes. By restricting our scope to early diverging crown teleost lineages, we avoid potentially confounding signals from highly nested radiations that substantially postdate both genome duplication and the origin of crown teleosts (26, 32). This approach provides a test of widely held assumptions about the nature of morphological evolution in teleosts and their holostean sister lineage.Open in a separate windowFig. 2.(A) Morphospace of Permian–Early Cretaceous crown Neopterygii. (B) One supertree subjected to our paleontological (Upper) and molecular (Lower) timescaling procedures to illustrate contrasts in the range of evolutionary timescales considered. Colors of points (A) and branches (B) indicate membership in major partitions of neopterygian phylogeny. Topologies are given in Datasets S4 and S5. See Dataset S6 for source trees.Open in a separate windowFig. S1.Morphospace of 398 Permian–Early Cretaceous Neopterygii. Three major axes of shape variation are presented. Silhouettes and accompanying arrows illustrate the main anatomical correlates of these principal axes, as described in Open in a separate windowFig. S2.Morphospace of 398 Permian–Early Cretaceous Neopterygii, illustrating the major clades of (A) teleosts and (B) holosteans.Open in a separate windowFig. S3.Comparisons of size rates between (A) holosteans and teleosts, (B) crown teleosts and all other neopterygians, (C) crown teleosts and stem teleosts, (D) crown teleosts and holosteans, and (E) stem teleosts and holosteans. Comparisons were made using the full-size SL dataset, a CS dataset, and a smaller SL dataset pruned to exactly match the taxon sampling of the CS dataset. Identical taxon sampling leads the CS and pruned SL datasets to yield near identical results. Although the larger SL dataset results often differ slightly, the overall conclusion from each pairwise comparison (i.e., which outcome is the most likely in an overall majority of trees) is identical in all but one comparison (E, under molecular timescales).Open in a separate windowFig. S4.Comparisons of size innovation between (A) holosteans and teleosts, (B) crown teleosts and all other neopterygians, (C) crown teleosts and stem teleosts, (D) crown teleosts and holosteans, and (E) stem teleosts and holosteans. Comparisons were made using the full-size SL dataset, a CS dataset, and a smaller SL dataset pruned to exactly match the taxon sampling of the CS dataset. Comparisons of size innovation are presented for K value distributions of the three datasets resemble each other closely.     

Emerging single-phase state in small manganite nanodisks

In complex oxides systems such as manganites, electronic phase separation (EPS), a consequence of strong electronic correlations, dictates the exotic electrical and magnetic properties of these materials. A fundamental yet unresolved issue is how EPS responds to spatial confinement; will EPS just scale with size of an object, or will the one of the phases be pinned? Understanding this behavior is critical for future oxides electronics and spintronics because scaling down of the system is unavoidable for these applications. In this work, we use La_0.325Pr_0.3Ca_0.375MnO₃ (LPCMO) single crystalline disks to study the effect of spatial confinement on EPS. The EPS state featuring coexistence of ferromagnetic metallic and charge order insulating phases appears to be the low-temperature ground state in bulk, thin films, and large disks, a previously unidentified ground state (i.e., a single ferromagnetic phase state emerges in smaller disks). The critical size is between 500 nm and 800 nm, which is similar to the characteristic length scale of EPS in the LPCMO system. The ability to create a pure ferromagnetic phase in manganite nanodisks is highly desirable for spintronic applications.Owing to strong coupling between spin, charge, orbital, and lattice (1, 2), different electronic phases often coexist spatially in strongly correlated materials known as electronic phase separation (EPS) (3, 4). For colossal magnetoresistance (CMR) manganites, EPS has been observed to have strong influence on the global magnetic and transport properties (5, 6). Regarding the physical origin of EPS, it has been shown theoretically that quenched disorder can lead to inhomogeneous states in manganites (1, 3, 7). Once long-range effects such as coulombic forces (8), cooperative oxygen octahedral distortions (9), or strain effects (10) are included, calculations show infinitesimal disorder (8, 11) or even no explicit disorder (10) may lead to EPS. Within a phenomenological Ginzburg–Landau theory, it has been shown that EPS is intrinsic in complex systems as a thermodynamic equilibrium state (12).Although the details of the origin of the EPS remain as a matter of dispute, its very existence as a new form of electronic state has been well accepted. The length scale of the EPS has been observed to vary widely from nanometers to micrometers depending on many parameters that can affect the competition between different electronic phases (13–20). It is thus of great interest to examine whether the EPS state still exists as the system is scaled down, especially when the spatial dimension of the system is smaller than the length scale of the EPS domains.In this work, we use La_0.325Pr_0.3Ca_0.375MnO₃ (LPCMO) as a prototype system to show a spatial confinement-induced transition from the EPS state to a single ferromagnetic phase state. The LPCMO system is chosen because of its well-known large length scale of EPS domains (approximately a micrometer) (21), which allows us to conveniently fabricate LPCMO epitaxial thin films into disks with diameters that are smaller than the EPS domain size. In LPCMO bulk (21) and thin films (6, 22), the EPS state was observed to be the low-temperature ground state. Using magnetic force microscope, we observe that the EPS state remains to be the ground state in disks with the size of 800 nm in diameter or larger but vanishes in the 500-nm-diameter disks whose size is distinctly smaller than the characteristic length scale of the EPS domains. In the 500-nm disks, only the ferromagnetic phase can be observed at all temperatures below Curie temperature T_c, indicating that the system is in a single-phase state rather than a EPS state. Our results further indicate that the large length scale EPS in the LPCMO system does not cost extra Coulomb energy, which otherwise should lead to a scaling down of EPS with decreasing size of the LPCMO disks (23, 24).LPCMO films with 60-nm thickness were epitaxially grown on SrTiO₃(001) substrates by pulsed-laser deposition. The substrates were kept at 780 °C in oxygen atmosphere of 5 × 10⁻³ millibars during growth. Unit cell by unit cell growth was achieved as indicated by oscillations of intensity of reflection high-energy electron diffraction (RHEED). The films were postannealed to 950 °C for 3 h in flowing oxygen to reduce oxygen vacancy and make sure that the films have the same magnetic properties as the bulk. The LPCMO disks with diameters from 500 nm to 20 μm were fabricated from the epitaxial thin films by electron beam lithography with a negative tone resist (for details, see the sample fabrication method and Fig. S1 in the Supporting Information). Magnetic properties of the LPCMO disk arrays were carried out using superconducting quantum interference device (SQUID) and magnetic force microscope (MFM) measurements.Open in a separate windowFig. S1.(A and B) Schematics of LPCMO disks samples for magnetic property measurement (A) and MFM mapping (B). (C) Optical microscopic image of LPCMO disk array with specific diameter. (D) SEM image of LPCMO disk sample for MFM imaging.A distinct signature of the EPS state in the LPCMO system is the thermal hysteresis for temperature-dependent magnetic and transport properties. Fig. 1 A–D shows temperature dependent magnetic properties of LPCMO disks with different diameters. To enhance the measuring signal for SQUID, we fabricate disk arrays for each selected diameter (the optical microscopic image shown in Fig. 1B, Inset for the 1-μm disk array). Fig. 1 A–D shows temperature-dependent magnetization measured under 1,000 Oe in-plane field for 7-μm, 1-μm, 800-nm, and 500-nm disk arrays, respectively. Thermal hysteresis can be observed for disk arrays with size down to 800 nm, reflecting the fact that ferromagnetic metallic (FMM) and charge order insulating (COI) phases coexist during the first-order phase transition (7, 25). For the 500-nm disk array, however, no thermal hysteresis can be observed. This observation implies that the EPS state may no longer exist in the system (7, 26).Open in a separate windowFig. 1.Temperature dependence of magnetization (black lines for cooling and red lines for warming) under 1,000 Oe (A–D) and initial magnetization (red lines) and hysteresis loop (black lines) (E–H) at 5 K of arrays of LPCMO disks with sizes of 7 μm (A and E), 1 μm (B and F), 800 nm (C and G), and 500 nm (D and H) in diameter and an area of 3 mm × 3 mm. (B, Inset) The optical microscopic image of d = 1 μm array. (C and D, Insets) Zoomed-in M vs. T loop around the thermal hysteresis region.The lack of EPS state in the 500-nm disk array is supported by the field-dependent magnetization measurements. Fig. 1 E–H shows in-plane initial magnetization curves and magnetic hysteresis loops (M-H loops) for the disk arrays measured at 5 K after zero-field cooling. For 800-nm or larger disk arrays, there is a clear difference between the initial magnetization curves and the corresponding M-H loops due to the coexistence of FMM and COI phases. When the magnetic field is applied from the initial state, the magnetization of the FMM phase first quickly aligns along the field direction, leading to the low field fast rise of the initial magnetization curve. With increasing field, the COI phase is melted and transits into the FMM phase. Once transited, the FMM phase will mostly stay even if the field is reduced, giving rise to the difference between initial magnetization curve and the M-H loop. The difference, however, becomes smaller with decreasing size. For the 500-nm disk array, the initial magnetization curve and the M-H loop virtually superimpose each other, indicating no melting of COI phase occurs. Both the temperature- and field-dependent magnetization measurements show a transition from the EPS state to a single FMM state with decreasing size of the disk, and the critical size should be between 500 nm and 800 nm.The transition from the EPS state to a single FMM state can be seen in MFM images shown in Fig. 2 (for MFM imaging details, see micromagnetic mapping method in Supporting Information). Fig. 2A shows morphological appearance of LPCMO disks with different sizes acquired by atomic force microscope (AFM). Fig. 2 B–D shows the corresponding MFM images of the LPCMO disks acquired at different temperatures under a perpendicular magnetic field of 1T. Here, the perpendicular magnetic field is applied to yield some perpendicular magnetization components for MFM imaging because the easy magnetization axis is in the plane. In the present color scale, the contrast below zero (red or black) represents FMM phase, whereas the contrast above zero (green or blue) represents nonferromagnetic phase [i.e., COI phase based on previous knowledge of the LPCMO system (21, 22)]. Apparently, except the 500-nm disk, all other disks show distinct features of the EPS state (i.e., the coexistence of the FMM and COI phases). Although the portion of FMM phase increases noticeably with decreasing temperature, the system stays in the EPS state even at 10 K. The typical length scale of the EPS domains is around a micrometer, which is consistent with previous reports (21, 27).Open in a separate windowFig. 2.(A) AFM images of LPCMO disks with sizes of 500 nm, 1 μm, 2 μm, 3.8 μm, 5 μm, and 7 μm in diameter. (B–D) The MFM images of LPCMO disks under 1T field (external magnetic field direction is pointing perpendicularly to the sample surface plane) taken at 10 K (B), 100 K (C), and 180 K (D). The sizes of disks in MFM images are adjusted and corrected to have same scales for each size with the help of scanning electron microscope (SEM) images (shown in Fig. S2) and dash lines show the approximate physical boundary of disks. The negative value in MFM image indicates attractive force and positive value indicates repulsive force.In stark contrast to the larger disks, the 500-nm disk does not exhibit any features of EPS in Fig. 2. Instead, the whole disk is in a ferromagnetic phase with a magnetization profile peaking in the center. To ensure that the EPS state is not diminished by the magnetic field applied during MFM imaging, we took MFM images of the 500-nm disk at 10 K under different perpendicular magnetic fields from 0T to 1T, as shown in Fig. 3A. At 0T, signals with opposite sign can only be seen on two sides of the disk along the marked line (MFM images of 4 disks shown in Fig. S3). This pattern is a typical MFM image for an in-plane ferromagnetic single domain, because only the two ends of an in-plane magnetic dipole yield perpendicular field gradient (with opposite signs) for the MFM tip to detect. Once a perpendicular field of 0.15T is applied, the in-plane magnetization is driven out of plane, leading to a center peaked MFM contour. The MFM signal increases with increasing field, as shown in Fig. 3B by the marked line profiles extracted from Fig. 3A.Open in a separate windowFig. 3.(A) MFM images of 500-nm disks at 10 K under different magnetic field. The lines show the path the line profile extracted. (B) Line profiles extracted from MFM images partly shown in Fig. 3A. (C) Simulated results corresponding to line-profiles in B. (D and E) The simulated Z-component of stray field 100 nm above sample disks and magnetic structure of 500-nm disks under a different magnetic field. In E, the direction and size of arrows show the direction and relative value of in-plane component of magnetization, and the color of disk indicates the value of Z-component of magnetization presented by ratio of Z-component to total magnetization, as shown in the color bar.Open in a separate windowFig. S3.MFM images of 500-nm disks taken at 10 K under zero field after zero-field cooling.The field-dependent behavior of the MFM contrast of the 500-nm disk is in qualitative agreement with micromagnetic simulations. Based on the MFM observation, the 500-nm disk is in an in-plane, single-domain state. Using this model as input, we performed micromagnetic simulation and obtained the Z-component of magnetic stray field distribution at 100 nm above sample surface (Fig. 3D; for details, see the micromagnetic simulation method in the Supporting Information), which is virtually the signal detected by MFM tip. The corresponding magnetic structures under different magnetic fields are shown in Fig. 3E. The marked line profiles extracted from simulation (Fig. 3D) are shown in Fig. 3C alongside with the experimental MFM line profiles (Fig. 3B). The subtle differences between Fig. 3B and Fig. 3C are likely caused by the fact that experimental MFM images are convoluted from signals of both the LPCMO disks and the MFM tips (∼100 nm in size). The consistency of MFM images and simulation confirms that the 500-nm disk is in a ferromagnetic single-domain state with an in-plane easy magnetization axis.Finally, we show that the 500-nm disk is in a single FMM state at all temperatures. Fig. 4 shows MFM images of the 500-nm disk acquired every 20 K, from 20 K to 200 K under 1,000 Oe. Other than the center-peaked FMM phase, no traces of COI phase can be observed. The MFM signal decreases with increasing temperature, which is consistent with the behavior of the temperature-dependent magnetization shown in Fig. 1. Considering the fact that we have never observed pure COI phase in the 500-nm disks, we believe this phenomenon may be caused by the existence of the ferromagnetic metallic edge state in the LPCMO system (22), which assists the 500-nm disk to be in pure ferromagnetic state when a single state is energetically preferred in the 500-nm disk.Open in a separate windowFig. 4.MFM images of 500-nm disks taken every 20 K, from 20 K to 200 K.In summary, we discovered a spatial confinement-induced transition from a EPS state featuring coexistence of FMM and COI phases to a single FMM state in the LPCMO system. The critical size for the transition is between 500 nm and 800 nm, which is similar to the characteristic length scale of the EPS state in the LPCMO system. Combining the MFM data and the micromagnetic simulation, we conclude that the 500-nm LPCMO disk is in a single-domain ferromagnetic state at all temperatures below T_c. A similar conclusion can be reached for 300-nm LPCMO disks (shown in Fig. S4), although it needs to be studied further whether a new state would emerge if the disk size becomes a few tens of nanometers or smaller. Our work opens a way to control EPS without external field or introducing strain and disorder, which is potentially useful to design electronic and spintronic devices in complex oxides systems.Open in a separate windowFig. S4.MFM images of 300-nm disks taken at different temperatures and magnetic fields.     

A 2,300-year-old architectural and astronomical complex in the Chincha Valley,Peru

Recent archaeological research on the south coast of Peru discovered a Late Paracas (ca. 400–100 BCE) mound and geoglyph complex in the middle Chincha Valley. This complex consists of linear geoglyphs, circular rock features, ceremonial mounds, and settlements spread over a 40-km² area. A striking feature of this culturally modified landscape is that the geoglyph lines converge on mounds and habitation sites to form discrete clusters. Likewise, these clusters contain a number of paired line segments and at least two U-shaped structures that marked the setting sun of the June solstice in antiquity. Excavations in three mounds confirm that they were built in Late Paracas times. The Chincha complex therefore predates the better-known Nasca lines to the south by several centuries and provides insight into the development and use of geoglyphs and platform mounds in Paracas society. The data presented here indicate that Paracas peoples engineered a carefully structured, ritualized landscape to demarcate areas and times for key ritual and social activities.The Chincha Valley, located 200 km south of Lima, was one of the largest and most productive regions of southern coastal Peru (Fig. 1). Previous research identified a rich prehispanic history in the valley, beginning at least in the early first millennium BCE and continuing through the Inca period in the 16th century CE (1–3). The earliest settled villages were part of the Paracas culture, a widespread political and social entity that began around 800 BCE and continued up to around 100 BCE. Previous field surveys identified at least 30 major Paracas period sites in the valley (1, 3, 4), making Chincha one of the main centers of development for this early Andean civilization (5). As such, it is an ideal area to test models of social evolution in general and to define the strategies that early peoples used to construct complex social organizations within the opportunities and constraints provided by their environments.Open in a separate windowFig. 1.Map showing location of the Chincha Valley, southern coastal Peru.Previous research demonstrated a dense Paracas settlement in the lower valley that focused on large platform mound complexes (Fig. 2) (4, 6). Three seasons of systematic, intensive survey and excavations by our team confirm the existence of a rich and complex Paracas occupation in the midvalley area as well, including both mound clusters and associated geoglyph features. In short, our data indicate that (i): the Chincha geoglyphs predate the better-known Nasca drainage ones by at least three centuries; (ii) Paracas period peoples created a complex landscape by constructing linear geoglyphs that converge on key settlements; and (iii) solstice marking was one component underlying the logic of geoglyph and platform mound construction and use in the Chincha Valley during the Paracas period.Open in a separate windowFig. 2.Distribution of archaeological sites linked to Paracas period settlement in Chincha, Peru. Redrawn from Canziani (4).     

Late Miocene episodic lakes in the arid Tarim Basin,western China

The Tibetan Plateau uplift and Cenozoic global cooling are thought to induce enhanced aridification in the Asian interior. Although the onset of Asian desertification is proposed to have started in the earliest Miocene, prevailing desert environment in the Tarim Basin, currently providing much of the Asian eolian dust sources, is only a geologically recent phenomenon. Here we report episodic occurrences of lacustrine environments during the Late Miocene and investigate how the episodic lakes vanished in the basin. Our oxygen isotopic (δ¹⁸O) record demonstrates that before the prevailing desert environment, episodic changes frequently alternating between lacustrine and fluvial-eolian environments can be linked to orbital variations. Wetter lacustrine phases generally corresponded to periods of high eccentricity and possibly high obliquity, and vice versa, suggesting a temperature control on the regional moisture level on orbital timescales. Boron isotopic (δ¹¹B) and δ¹⁸O records, together with other geochemical indicators, consistently show that the episodic lakes finally dried up at ∼4.9 million years ago (Ma), permanently and irreversibly. Although the episodic occurrences of lakes appear to be linked to orbitally induced global climatic changes, the plateau (Tibetan, Pamir, and Tianshan) uplift was primarily responsible for the final vanishing of the episodic lakes in the Tarim Basin, occurring at a relatively warm, stable climate period.Once part of the Neo-Tethys Sea, indicated by the Paleogene littoral-neritic deposits with intercalated marine strata (1), the Tarim Basin in western China (Fig. 1) has experienced dramatic environmental and depositional changes during the Cenozoic. The eventual separation of the basin from the remnant sea probably occurred during the middle to late Eocene (1, 2), and since then, terrestrial sedimentation and environment have prevailed. Today, the basin has relatively flat topography with elevation ranging between 800 and 1,300 m above sea level (asl), surrounded by high mountain ranges with average elevation exceeding 4,000 m. Hyperarid climate prevails with mean annual precipitation <50 mm and evaporation ∼3,000 mm. Active sand dunes occupy 80% of the basin, forming the Taklimakan Desert, the largest desert in China and the second largest in the world (3). Thick desert deposits in the basin also provide a major source for dust storms occurring in East Asia (4).Open in a separate windowFig. 1.Digital elevation model of the Tarim Basin and surrounding mountain ranges. The studied 1,050-m sediment core is retrieved from Lop Nor (1), at a relatively low elevation (∼800 m asl) in the eastern basin. Also indicated are locations of previous studied sections, near Sanju (2), Kuqa (3), and Korla (4) in the basin.The Tibetan Plateau uplift, long-term global cooling, and the associated retreat of the remnant sea during the Cenozoic, through their complex interplay, may have all contributed to the enhanced aridification in the Asian interior and eventually the desert formation (2, 5–11). Although the onset of Asian desertification is proposed to have started in the earliest Miocene (12, 13), prevailing desert environment in the Tarim Basin, currently providing much of the Asian eolian dust sources (4), is only a geologically recent phenomenon (9). Previous lithological and pollen studies (9, 14, 15) suggest that the currently prevailing desert environment started probably at ∼5 Ma. However, what kind of environmental conditions prevailed before that and how the dramatic changes occurred largely remain elusive.To better decipher the aridification history, we used a 1,050-m-long, continuous sediment core retrieved from Lop Nor (39°46''0''''N, 88°23''19''''E) in the eastern basin (Fig. S1). The core site has a relatively low elevation (∼800 m asl) (Fig. 1) and thus effectively records basin environment. Previous study sites came from elevated basin margins (9, 14). The core mainly consists of lacustrine sediments with associated fluvial-eolian sands (Fig. 2). The core chronology was established previously based on paleomagnetic polarity (15). The total 706 remanence measurements on the continuous sediment profile allow straightforward correlation with the CK95 geomagnetic polarity timescale (16) and identification of 14 normal (N1–N14) and 13 reversal (R1–R13) polarity zones over the last 7.1 Ma (Fig. S1). The CK95 timescale is largely consistent (within a few ky) with that inferred from marine archives (17) over the last 5.23 Ma, and before 5.23 Ma a practical measure of chronological uncertainty was estimated to be within ∼100 ky (18), assuring that the terrestrial records can be directly compared with marine records and orbital changes at the timescale of >100 ky. The derived chronology yields an average sedimentation rate of ∼200 m/Ma at lower sections (>1.77 Ma), allowing high-resolution studies of the desertification history at the critical interval.Open in a separate windowFig. 2.Records of δ¹¹B, δ¹⁸O, TOC, and CaCO₃ changes in the Lop Nor profile. Representative photos show lithological changes mainly from lacustrine bluish gray argillaceous limestone to fluvial-eolian brown/red clayey silt, as shown in the lithological column with visual colors indicated (15). High levels and large fluctuations of proxy values occurred only before ∼4.9 Ma, and since then, those remain low and within a small range, largely similar to modern conditions. The transition from episodic occurrences of lacustrine phases to prevailing desert environments thus appears to be permanent and irreversible.Boron and oxygen isotopes from carbonates are important environmental indicators (19–25) and used here to infer the environmental evolution in the Tarim Basin. We also present total organic carbon (TOC), calcium carbonate (CaCO₃), ostracod, and grain size records to substantiate the isotopic evidence (SI Materials and Methods). Our δ¹¹B profile shows substantial, stepwise changes over the last ∼7.1 Ma (Fig. 2). Carbonate δ¹¹B values were −5.0 ± 1.6‰ (n = 3) before 6.0 Ma, increased rapidly to ∼11‰ at 4.9–6.0 Ma, and then stayed at roughly the same level (10.7 ± 2.2‰, n = 25) for the remaining 4.9 Ma. Higher-resolution δ¹⁸O, TOC, and CaCO₃ profiles generally confirm the pattern observed in the low-resolution δ¹¹B one (Fig. 2). The δ¹⁸O values remained low, ranging from −10‰ to −4‰ over the last 4.9 Ma. However, δ¹⁸O values frequently oscillated between −10‰ and 5‰ before that. Similarly, the TOC profile shows consistently low organic carbon content (0–0.2%) after 4.9 Ma and large fluctuations (0–1.0%) before then. The CaCO₃ profile also indicates consistently low values (0–25%) after 4.9 Ma and large fluctuations (0–50%) earlier (Fig. 2).The multiple proxy records strongly suggest that critical environmental changes must have occurred at ∼4.9 Ma. δ¹¹B values of carbonates from marine sources differ substantially from those of nonmarine carbonates (19–21). δ¹¹B values after 4.9 Ma are close to those from marine carbonates, but values before 6 Ma fall into the range of lacustrine carbonates (22). Positive δ¹⁸O values before 4.9 Ma also indicate lacustrine environments at that time. Carbonates from modern lakes in arid and semiarid regions of northwestern China show similar positive δ¹⁸O values (23), due to strong evaporation processes. High TOC and CaCO₃ contents (Fig. 2) further support that lacustrine environments existed in the basin before ∼4.9 Ma. δ¹⁸O values after 4.9 Ma are comparable to those in Cenozoic soil carbonates (24) and ancient marine carbonates in the Tarim Basin (25). However, the accompanying carbonate δ¹³C values throughout the record, ranging from −4‰ to 1‰ (Dataset S1), are significantly higher than those from Cenozoic soil carbonates reported (26), essentially ruling out the possibility of soil carbonate source. Using modern prevailing desert environment in the basin as an analog, the combined δ¹¹B and δ¹⁸O evidence thus suggests that the sediment deposits in the basin after 4.9 Ma must be eolian-fluvial in origin and their sources, at least carbonate grains, came from weathered ancient marine carbonates in nearby regions.Sedimentological and stratigraphic patterns in other exposed sections from different parts of the basin (9, 14) share great similarity with the Lop Nor core profile (Fig. S1). Episodic lacustrine mudstones and/or siltstones during the Late Miocene were present in all sections and were replaced by fluvial-eolian deposits later. Studies of ostracod assemblages (27) also suggest a shallow paleolake with brackish water environments in the northern basin during the Late Miocene. Changes in the depositional environment from our Lop Nor profile alone could be plausibly explained by a shift in basin center due to tectonic compressions, as evidenced from the slightly uplifted central basin (Fig. 1). However, similar temporal changes occurring basin-wide at ∼4.9 Ma argue against it. Instead, our results, together with previous studies (5, 14, 15, 27), suggest that paleolakes were widely present in the low lands of the basin during the Late Miocene, much different from currently prevailing desert environments with a few scattered small lakes. The existing evidence, although still limited (Fig. 1), would point to the occurrence of a possible megalake in the Tarim Basin during the Late Miocene.Three high-resolution records, δ¹⁸O, TOC, and CaCO₃, further suggest that lacustrine environments before ∼4.9 Ma were not permanent (Fig. 2). These large fluctuations indicate frequent switches between lacustrine and fluvial-eolian environments in the basin. High proxy values, δ¹⁸O in particular, appear to indicate lacustrine environments, whereas low values, similar to ones after 4.9 Ma, correspond to fluvial-eolian deposits. This is consistent with lithological features at this interval, showing argillaceous limestone intercalated with clayey layers (15), the occurrence of ostracod assemblages (Fig. 3) from lacustrine sediments, grain size changes (Fig. S2), and detrital carbonate grains identified in photomicrographs of fluvial-eolian deposits (Fig. S3).Open in a separate windowFig. 3.δ¹⁸O fluctuations linked to eccentricity and obliquity orbital variations at 4.5–7.1 Ma. Lacustrine phases (high δ¹⁸O) generally correspond to periods of high eccentricity and obliquity. Fluvial-eolian environment (low δ¹⁸O, highlighted with gray bars), developed more around 6.5, 6.1, 5.7, and 5.2 Ma, at a ∼400-ky eccentricity beat. The red bar indicates the last occurrence of lacustrine environment at ∼4.9 Ma. Occurrences of ostracod assemblages with total number >100 are also indicated.To further investigate such episodic changes, we performed spectral analysis on the δ¹⁸O record over the interval 4.5–7.1 Ma. Strong spectral power at orbital frequencies were identified, with periods of ∼400 ky throughout the interval, ∼41 ky particularly at 6–6.5 Ma, and ∼100 ky at 5–6 Ma and 6.5–7.1 Ma (Fig. S4). Precessional ∼20-ky power might also have existed but was relatively weak and discontinuous. Lacustrine phase as indicated by high δ¹⁸O values and ostracod assemblages generally occurred at periods of high eccentricity and obliquity (28) (Fig. 3). At 6–6.5 Ma, δ¹⁸O shows clear correspondence to orbital obliquity variation (Fig. 3A). Additionally, the number of low δ¹⁸O values (fluvial-eolian environment) occurred more around 6.5, 6.1, 5.6, and 5.2 Ma, at a ∼400-ky beat following orbital eccentricity variation (Fig. 3B). The cluster of high δ¹⁸O values (>0‰) at 4.9–5.0 Ma signals the last occurrence of lacustrine environments in the basin. Our orbital association thus allows us to precisely determine the timing of desert formation at ∼4.9 Ma, ∼400 ky (an eccentricity cycle) later than the age inferred from basin margins (9, 14) and yet all occurring at eccentricity minima (Fig. 3B). As high eccentricity and obliquity generally correspond to warm conditions at orbital timescales, lacustrine (wet) phase could be associated with warm periods, consistent with the notion that cooler conditions would reduce moisture in the atmosphere and enhance continental drying (10, 11). We recognize that the chronological uncertainty from the geomagnetic polarity timescale before 5.23 Ma, within ∼100 ky (18), could confound our association of wet phase with high obliquity, although it is unlikely affected at the 400-ky eccentricity beat. However, the opposite association, wet phase with low obliquity, and the combination with high eccentricity, would require a different, yet unknown mechanism that is inconsistent with the orbital theory of Pleistocene ice ages.Superimposed on the orbitally episodic changes, the δ¹⁸O record also shows a long-term trend of deteriorating lacustrine conditions at 4.9–7.1 Ma. As δ¹⁸O values indicate two depositional environments, lacustrine and fluvial-eolian, the range of δ¹⁸O changes (between −10‰ and 5‰) does not vary much over this period (Fig. 3). Rather, the duration of high δ¹⁸O vs. low values would reflect the long-term trend. The mean δ¹⁸O values over 40-ky and 400-ky intervals both show a decreasing trend, with dominant lacustrine phase before 6.1 Ma, more developed fluvial-eolian environment at 5.7–6.1 Ma, a return to slightly better lacustrine environment at 5.3–5.7 Ma, and lacustrine phase permanently vanished around 4.9 Ma (Fig. 4).Open in a separate windowFig. 4.Long-term δ¹⁸O changes compared with global climatic conditions and regional tectonic activities. Global benthic δ¹⁸O (29) and northwestern Pacific SST (30) records show minimal long-term changes at 4–7 Ma, whereas occurrences of detrital apatite fission track ages (35) in northern western Kunlun peaked, and the Tarim episodic lakes gradually vanished. The dark thick lines are their 400-ky running means, and the light green line on δ¹⁸O represents the 40-ky running mean.The gradual disappearance of the Tarim episodic lakes could be potentially explained by the two driving forces, plateau uplift (the Tibetan Plateau, Pamir Plateau, and Tianshan) and long-term global cooling. However, the long-term global climate was relatively warm and stable during this period (Fig. 4). The global benthic δ¹⁸O record (29) shows that much of the Miocene cooling occurred between 15 and 11 Ma, and the cooling between 8 and 5 Ma was minimal. Supporting this view, sea surface temperature records from the northwestern Pacific (30) show that almost no cooling occurred between 6 and 4 Ma (perhaps further to 3 Ma) (Fig. 4). Particularly, the mean global climate (29, 30) was even warmer at 4.1–4.5 Ma than at 5.7–6.1 Ma, whereas lacustrine phase permanently disappeared after ∼4.9 Ma (Fig. 4), indicating decoupling of the lake evolution from global climate. Therefore, long-term global cooling might have played a subordinate role in the lake disappearance.Instead, the growth of surrounding mountain ranges (Tibetan, Pamir, and Tianshan) may have blocked moisture from the west and south, changed air circulations, and eventually led to the permanent lake disappearance within the basin. Today, the Tarim Basin receives limited moisture from westerlies (31) through the Pamir and Tianshan Ranges (and perhaps from the Indian Ocean in summertime as well). Although the Indo–Eurasian convergence since the Late Eocene resulted in high elevations of the Tibetan Plateau and, to a lesser degree, surrounding mountains including Pamir and Tianshan by the mid-Miocene time (8), tectonic activities in broad areas around the Tarim Basin appear to be rejuvenated since the Late Miocene. Tectonic deformations during the Late Miocene–Early Pliocene inferred from growth strata, sedimentary facies changes, and low-temperature thermochronologic studies occurred in Tianshan to the north of the Tarim Basin, in the Kunlun Mountains to its south and the Pamir to its west (32–35). Syntectonic growth strata from the foreland basins of the Kunlun and Tianshan Ranges (32) show that strong crust shortening and potential mountain uplift initiated ∼6.5–5 Ma and lasted to the Early Pleistocene (Fig. S5), similarly reported in northern Pamir (33). Cenozoic sequences in the Pamir–Tianshan convergence zone, changing from an arid continental plain to an intermountain basin by ∼5 Ma, support surface uplift of the west margin of the Tarim Basin (34). Occurrence of detrital apatite fission track ages from West Kunlun Ranges (35) also peaked at ∼4.5 Ma (Fig. 4). Thus, the reactivated uplift of Pamir and West Kunlun Ranges and northward movement of Pamir during the Late Miocene–Pliocene would progressively block moisture into the basin and enhance regional aridity to a certain threshold to terminate lacustrine environments even during warm periods with favorable orbital configuration. Although inconsistencies indeed exist in linking regional climatic and environmental changes to tectonic events (8), the final vanishing of the Tarim episodic lakes is better explained by tectonic factors.Therefore, our multiple-proxy results consistently show that the Taklimakan sand sea began to form at ∼4.9 Ma and that the transition into prevailing desert environment was permanent and irreversible. The episodic occurrences of lacustrine environments at favorable climatic conditions (warm periods) during the Late Miocene suggest that uplifted mountain ranges then were not high enough to effectively block moisture from being transported into the basin. The transition from episodic lacustrine environments to prevailing desert deposits was gradual, from 7.1 Ma (or earlier beyond our record) to 4.9 Ma, for which we suggest that although long-term global cooling enhanced the overall aridification in the Asian interior during the Late Cenozoic (10, 11), plateau uplift played a more important role in finally drying up the episodic lakes within a relatively warm, stable climate period, thus decoupling regional climate temporally from a global trend. Our high-resolution records thus demonstrate that regional climate in the Tarim Basin reached a critical state in the Late Miocene, with the dual effects from global climate conditions and regional tectonic settings then. With global climate remaining relatively stable and warm entering the Pliocene, by ∼3–4 Ma (Fig. 4), drying up of episodic lakes at ∼4.9 Ma could thus be largely attributed to rejuvenated tectonic activities.