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).     

Groundwater sapping as the cause of irreversible desertification of Hunshandake Sandy Lands,Inner Mongolia,northern China

In the middle-to-late Holocene, Earth’s monsoonal regions experienced catastrophic precipitation decreases that produced green to desert state shifts. Resulting hydrologic regime change negatively impacted water availability and Neolithic cultures. Whereas mid-Holocene drying is commonly attributed to slow insolation reduction and subsequent nonlinear vegetation–atmosphere feedbacks that produce threshold conditions, evidence of trigger events initiating state switching has remained elusive. Here we document a threshold event ca. 4,200 years ago in the Hunshandake Sandy Lands of Inner Mongolia, northern China, associated with groundwater capture by the Xilamulun River. This process initiated a sudden and irreversible region-wide hydrologic event that exacerbated the desertification of the Hunshandake, resulting in post-Humid Period mass migration of northern China’s Neolithic cultures. The Hunshandake remains arid and is unlikely, even with massive rehabilitation efforts, to revert back to green conditions.Earth’s climate is subject to abrupt, severe, and widespread change, with nonlinear vegetation–atmosphere feedbacks that produced extensive and catastrophic ecosystem shifts and subsequent cultural disruption and dispersion during the Holocene (1–7). In the early and middle Holocene, northern China’s eastern deserts, including much of the currently sparsely vegetated and semistabilized Hunshandake (Figs. 1 and ​and2),2), were covered by forests (8), reflecting significantly wetter climate associated with intensification of monsoon precipitation by up to 50% (6).Open in a separate windowFig. 1.Geographical location of the Hunshandake Sandy Lands (A) and its area (encircled by red line in B). The black rectangle in B marks the location of the enlarged maps C and D on the Right, and the green rectangle shows the location of Fig. 2. Map C shows the localities of water samples, and map D shows the localities of sections with stratigraphy presented in Fig. 3. The sand–paleosol section P (Fig. 3) is on the southern margin, and the site Bayanchagan marks the coring site of ref. 8. Rivers with headwaters in the Hunshandake likely formed by groundwater sapping are marked in blue. Drainages to the southwest and west are currently undergoing groundwater sapping, with substantial spring-driven flow found at the current river base level.Open in a separate windowFig. 2.(Left) Holocene lakes and channels in the Hunshandake and lake extent at selected epochs. Upper, middle, and lower lakes are indicated by points A, B, and C, respectively. Xilamulun River (point D) drains to the east. Groundwater-sapping headcuts at the upper reaches of incised canyons (point E) suggest a mid-Holocene interval of easterly surface flow, followed by groundwater drainage beginning at the ca. 4.2 ka event. Northern and central channels at point E are currently abandoned, and groundwater sapping has migrated to the southerly of the three channels shown. (Right) Cross-sections of the predrainage shift, northerly drainage into Dali Lake (Localities shown on the Left), showing the increase in widths of channels downstream (Vertical exaggeration ∼30:1).Monsoonal weakening, in response to middle-to-late Holocene insolation decrease, reduced precipitation, leading to a green/sandy shift and desertification across Inner Mongolia between ca. 5,000 and 3,000 y (years) ago (6). However, variations in the timing of this transition (9, 10) suggest local/regional thresholds or possibly environmental tipping by stochastic fluctuations. The impacts of this wet-to-dry shift in the Hunshandake, expressed as variations in surface and subsurface hydrology coincident with the termination of the formation of thick and spatially extensive paleosols, and the impacts of a ca. 4.2 ka (1 ka = 1,000 years) mid-Holocene desiccation of the Hunshandake on the development of early Chinese culture remain poorly understood and controversial (6, 11). Here we report for the first time to our knowledge on variability in a large early-to-middle Holocene freshwater lake system in China’s Hunshandake Sandy Lands and associated vegetation change, which demonstrates a model of abrupt green/desert switching. We document a possible hydrologic trigger event for this switching and discuss associated vegetation and hydrologic disruptions that significantly impacted human activities in the region.     

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.     

Amphitheater-headed canyons formed by megaflooding at Malad Gorge,Idaho

Many bedrock canyons on Earth and Mars were eroded by upstream propagating headwalls, and a prominent goal in geomorphology and planetary science is to determine formation processes from canyon morphology. A diagnostic link between process and form remains highly controversial, however, and field investigations that isolate controls on canyon morphology are needed. Here we investigate the origin of Malad Gorge, Idaho, a canyon system cut into basalt with three remarkably distinct heads: two with amphitheater headwalls and the third housing the active Wood River and ending in a 7% grade knickzone. Scoured rims of the headwalls, relict plunge pools, sediment-transport constraints, and cosmogenic (³He) exposure ages indicate formation of the amphitheater-headed canyons by large-scale flooding ∼46 ka, coeval with formation of Box Canyon 18 km to the south as well as the eruption of McKinney Butte Basalt, suggesting widespread canyon formation following lava-flow diversion of the paleo-Wood River. Exposure ages within the knickzone-headed canyon indicate progressive upstream younging of strath terraces and a knickzone propagation rate of 2.5 cm/y over at least the past 33 ka. Results point to a potential diagnostic link between vertical amphitheater headwalls in basalt and rapid erosion during megaflooding due to the onset of block toppling, rather than previous interpretations of seepage erosion, with implications for quantifying the early hydrosphere of Mars.Landscapes adjust to perturbations in tectonics and base level through upstream propagation of steepened river reaches, or knickzones, thereby communicating environmental signals throughout a drainage basin (e.g., ref. 1). Nowhere are knickzones more important and apparent than in landscapes where canyon heads actively cut into plateaus, such as tributaries of the Grand Canyon, United States, and the basaltic plains of Mars (e.g., refs. 2–4). Here the stark topographic contrast between low-relief uplands and deeply incised canyons sharply delineates canyon rims and planform morphology. Canyon heads can have varied shapes from amphitheaters with vertical headwalls to more pointed planform shapes with lower gradients, and a prominent goal in geomorphology and planetary science is to link canyon morphology to formation processes (e.g., refs. 4–8), with implications for understanding the history of water on Mars.Amphitheater-headed canyons on Mars are most likely cut into layered basalt (9, 10), and canyon-formation interpretations have ranged widely from slow seepage erosion to catastrophic megafloods (4–6, 11, 12). Few studies have been conducted on the formation of amphitheater-headed canyons in basalt on Earth, however, and instead, terrestrial canyons in other substrates are often used as Martian analogs. For example, groundwater sapping is a key process in forming amphitheater-headed canyons in unconsolidated sand (e.g., refs. 8, 13, 14), but its importance is controversial in rock (5, 12, 15). Amphitheater-headed canyons are also common to plateaus with strong-over-weak sedimentary rocks (3, 16); however, here the tendency for undercutting is so strong that canyon-head morphology may bear little information about erosional processes, whether driven by groundwater or overland flow (e.g., refs. 3, 5, 17). Canyons in some basaltic landscapes lack strong-over-weak stratigraphy, contain large boulders that require transport, and show potential for headwall retreat by block toppling (18–21), all of which make extension of process–form relationships in sand and sedimentary rocks to basalt and Mars uncertain.To test the hypothesis of a link between canyon formation and canyon morphology in basalt, we need field measurements that can constrain formation processes for canyons with distinct morphologies, but carved into the same rock type. Here we report on the origin of Malad Gorge, a canyon complex eroded into columnar basalt with markedly different shaped canyon heads. Results point to a potential diagnostic link between canyon-head morphology and formative process by megaflood erosion in basalt.Malad Gorge is a tributary to the Snake River Canyon, Idaho, within the Snake River Plain, a broad depression filled by volcanic flows that erupted between ∼15 Ma and ∼2 ka (22, 23). The gorge sits at the northern extent of Hagerman Valley, a particularly wide (∼7 km) part of the Snake River Canyon (Fig. 1). Malad Gorge is eroded into the Gooding Butte Basalt [⁴⁰Ar/³⁹Ar eruption age: 373 ± 12 ka (25)] which is composed of stacked lava beds, each several meters thick with similar well-defined columns bounded by cooling joints and no apparent differences in strength between beds. The Wood (or Malad) River, a major drainage system from the Sawtooth Range to the north, drains through Malad Gorge before joining the Snake River. The Wood River is thought to have been diverted from an ancestral, now pillow lava-filled canyon into Malad Gorge by McKinney Butte basalt flows (24) [⁴⁰Ar/³⁹Ar eruption age: 52 ± 24 ka (25) (Fig. 1).Open in a separate windowFig. 1.Shaded relief map of the study region (50-m contour interval) showing basalt flows (23), their exposure age sample locations, and the path of the ancestral Wood River following Malde (24) (US Geological Survey).Malad Gorge contains three distinct canyon heads herein referred to as Woody’s Cove, Stubby Canyon, and Pointed Canyon (Fig. 2A). Woody’s Cove and Stubby have amphitheater heads with ∼50-m-high vertical headwalls (Fig. 2C), and talus accumulation at headwall bases indicates long-lived inactive fluvial transport (Fig. 3 A and B). Woody’s Cove, the shortest of the three canyons, lacks major spring flows and has minor, intermittent overland flow partially fed by irrigation runoff that spills over the canyon rim. Stubby has no modern-day overland flow entering the canyon, and springs emanate from a pool near its headwall (Fig. 3B). In contrast, Pointed Canyon is distinctly more acute in planform morphology, contains a 7% grade knickzone composed of multiple steps rather than a vertical headwall (Figs. 2C and ​and3C),3C), and extends the farthest upstream.Open in a separate windowFig. 2.Malad Gorge topography (10-m contour interval) and aerial orthophotography (US Geological Survey). (A) Overview map and (B) close-up for Stubby and Pointed canyons showing mapped bedrock scours (white arrows), exposure age sample locations (red circles) with age results, location of the uppermost active knickpoint (black circle), abandoned bedrock channels (blue dashed lines), and grain-size analysis sites (blue squares). The blue star shows the reconstructed location of the headwall of Pointed Canyon at 46 ka (see Discussion and Fig. 5). (C) Longitudinal profile along Stubby and Pointed canyons from their confluence (shown as white lines in B) with local slope, S, averaged over regions demarked by dashed lines (Fig. 4A shows close-up of profile in Stubby Canyon).Open in a separate windowFig. 3.Photographs of (A) headwall of Woody’s Cove (person for scale, circled), (B) ∼50-m-high headwall of Stubby Canyon, (C) downstream-most waterfall at Pointed Canyon knickzone (12-m-high waterfall with overcrossing highway for scale), (D) fluted and polished notch at the rim of Stubby Canyon (notch relief is 10 m), (E) upstream-most waterfall at Pointed Canyon knickzone (within the southern anabranch of Fig. S2), and (F) upstream-most abandoned channel in Fig. 2B and Fig. S2 (channel relief is ∼10 m). White coloring on the headwalls in A and B is likely residue from irrigation runoff.Early work attributed the amphitheater-headed canyons in this region—Malad Gorge, Box Canyon, located 18 km south of Malad Gorge (Fig. 1), and Blue Lakes Canyon located 42 km to the SE—to formation by seepage erosion because of no modern overland flow and the occurrence of some of the largest springs in the United States in this region (7). Because spring flows (e.g., ∼10 m³/s in Box Canyon; US Geological Survey gauge 13095500) are far deficient to move the boulders that line the canyon floors, Stearns (7) reasoned that the boulders must chemically erode in place. This explanation is improbable, however, given the young age of the Quaternary basalt (25), spring water saturated in dissolved solids (19), and no evidence of rapid chemical weathering (e.g., talus blocks are angular and have little to no weathering rinds). Instead of groundwater sapping, Box Canyon was likely carved by a large-scale flood event that occurred ∼45 ka based on ³He cosmogenic exposure age dating of the scoured rim of the canyon headwall (19, 26). In addition, Blue Lakes Canyon was formed during the Bonneville Flood [∼18–22 ka (27, 28)], one of the world’s largest outburst floods that occurred as a result of catastrophic draining of glacial lake Bonneville (21). In both cases, canyon formation was inferred to have occurred through upstream headwall propagation by waterfall erosion.Herein we aim to test whether the amphitheater-headed canyons at Malad Gorge also owe their origin to catastrophic flooding, whether Pointed Canyon has a different origin, and whether canyon morphology is diagnostic of formation process. To this end we present field observations, sediment-size measurements, hydraulic modeling, and cosmogenic exposure ages of water-scoured rock surfaces and basalt-flow surfaces (Methods and Tables S1 and S2).     

Rapid seafloor changes associated with the degradation of Arctic submarine permafrost

Repeated high-resolution bathymetric surveys of the shelf edge of the Canadian Beaufort Sea during 2- to 9-y-long survey intervals reveal rapid morphological changes. New steep-sided depressions up to 28 m in depth developed, and lateral retreat along scarp faces occurred at multiple sites. These morphological changes appeared between 120-m and 150-m water depth, near the maximum limit of the submerged glacial-age permafrost, and are attributed to permafrost thawing where ascending groundwater is concentrated along the relict permafrost boundary. The groundwater is produced by the regional thawing of the permafrost base due to the shift in the geothermal gradient as a result of the interglacial transgression of the shelf. In contrast, where groundwater discharge is reduced, sediments freeze at the ambient sea bottom temperature of ∼−1.4 °C. The consequent expansion of freezing sediment creates ice-cored topographic highs or pingos, which are particularly abundant adjacent to the discharge area.

The effects of on-going terrestrial permafrost degradation (1–3) have been appraised by comparison of sequential images of Arctic landscapes that show geomorphic changes attributed primarily to thermokarst activity induced by recent atmospheric warming and ongoing natural periglacial processes (4–8). While the existence of extensive relict submarine permafrost on the continental shelves in the Arctic has been known for years (9, 10), the dynamics of submarine permafrost growth and decay and consequent modifications of seafloor morphology are largely unexplored.Throughout the Pleistocene, much of the vast continental shelf areas of the Arctic Ocean experienced marine transgressions and regressions associated with ∼125-m global sea level changes (11). Extensive terrestrial permafrost formed during sea-level low stands when the mean annual air temperatures of the exposed shelves were less than −15 °C (11, 12). Exploration wells drilled on the continental shelf in the Canadian Beaufort Sea show that relict terrestrial permafrost occurs in places to depths >600 m below seafloor (mbsf) and forms a seaward-thinning wedge beneath the outer shelf (10, 13, 14) (Fig. 1 A and B). The hydrography of the Canadian Beaufort Sea slows the degradation of the relict permafrost because a cold-water layer with temperatures usually near -1.4 °C blankets the seafloor from midshelf depths down to ∼200-m water depth (mwd) (15, 16) (Fig. 1B). As the freezing point temperature of interstitial waters is also controlled by salinity and sediment grain size, partially frozen sediments occur in a zone delimited by the ∼−2 °C and 0 °C isotherms (17, 18) (Fig. 1B).Open in a separate windowFig. 1.Map and cross-section showing the relationship between shelf edge morphology and the subsurface thermal structure along the shelf edge in the Canadian Beaufort Sea. (A) Shows the location of the study area with respect to estimates of submarine permafrost density (see key) and thickness, modified after 14. Thin contours indicate permafrost thickness in meters. Thicker contour is 120-m isobath marking the shelf edge. Area of repeat mapping coverage shown in Fig. 2 is indicated with a red box. (B) Shows a schematic cross-section with contours of selected subsurface isotherms modified after 15 along line x-x’ in A. The dotted blue line illustrates a thermal minimum (T-min) running through the relict permafrost isotherm and Beaufort Sea waters (16). Green shading indicates relict permafrost. Turquoise arrows show inferred flow of water from permafrost thawing along the base of the relict permafrost to the seafloor. The brown area indicates the zone where relict Pleistocene permafrost is predicted to have thawed with consequent movement of liberated groundwater, associated latent heat transfer and thaw consolidation causing surface settlement. Dashed brown lines define the subbottom limits for methane hydrate stability zone (MHSZ) which starts at ∼240 m below the sea surface and extends into the subsurface depending on the pressure and temperature gradient. The red box indicates the area shown in more detail in C with the same color scheme. The area of denuded seafloor in C is flanked by PLFs (dark-blue fill). Red arrows indicate the direction of heat transfer along the seaward edge of relict permafrost wedge.Distinctive surface morphologies characterize terrestrial permafrost areas. Conical hills (3 to 100 m in diameter) called pingos are common in the Arctic (19, 20). Pingos are formed due to freezing of groundwater. They characteristically contain lenses of nearly pure ground ice that cause heaving of the ground surface. Positive relief features with similar dimensions, referred to as pingo-like features (PLFs), are scattered across the Canadian Beaufort shelf (21, 22). On land, permafrost thawing, where there is ground ice in excess of the sediment pore space, can induce sediment consolidation (23), and surface subsidence results in widespread thermokarst landforms. Among the more dramatic occurrences are retrogressive thaw slumps (4–8, 24). These form where ice-rich permafrost experiences surface thaw causing thaw settlement and release of liquified sediment flows. Because of the loss of volume associated with thawing of massive ground ice, thaw slumps can quickly denude permafrost landscapes.During the first systematic multibeam mapping surveys in 2010 covering part of the shelf edge and slope in the Canadian Beaufort Sea, a band of unusually rough seafloor morphology between ∼120 and ∼200 mwds (25) was discovered along a ∼95-km-long stretch of the shelf. Subsequently, three additional multibeam surveys covering small characteristic areas (Fig. 2A) were conducted to understand the processes responsible for the observed morphologies. Here, we document the unique morphologies and seafloor change in this area and explore how the seafloor features may be related to subsea permafrost degradation and formation.Open in a separate windowFig. 2.(A) Shows bathymetry of a small section of the shelf edge indicated in Fig. 1A, with a color scale going from white (128 m) to blue (200 m) and contours at 120, 140, 170, and 200 mbsf. Outlines of areas resurveyed in 2013 (blue), 2017 (turquoise), and 2019 (green) are superimposed on the 2019 and regional 2010 survey. Colored symbols indicate locations of cores with porewater data using same key as Fig. 5B. The red star indicates the location of a temperature tripod deployed in the period 2015 to 2016. The location of ROV dive tracks (blue paths) are indicated. (B) Covers the same area as A with polygons identifying sites where changes were noted between surveys as follows: 2010 to 2013 (green), 2013 to 2017 (purple), 2017 to 2019 (black), and 2010 and 2019 (red). (C) Shows the same area, colored according to the difference in bathymetry between the 2019 survey and an idealized smooth surface extending between the top of the shelf edge scarp and the layered sediments occurring between the numerous PLFs. This is used to estimate the volume of material that eroded assuming the earlier Holocene seafloor corresponded with this idealized surface. Three zones of topography are labeled. Red boxes are locations of Figs. 3 and ​and5.5. (D) Shows Chirp profiles with the position of profiles shown in C and Fig. 5A. Light-green backdrop in X-X’ indicates possible void produced by retrogressive slide retreat used to calculated volume loss. Also indicated are TL, tilted layers; P, pingo-like-feature; and DR, diffuse reflector.     

Greenland ice sheet motion insensitive to exceptional meltwater forcing

Changes to the dynamics of the Greenland ice sheet can be forced by various mechanisms including surface-melt–induced ice acceleration and oceanic forcing of marine-terminating glaciers. We use observations of ice motion to examine the surface melt–induced dynamic response of a land-terminating outlet glacier in southwest Greenland to the exceptional melting observed in 2012. During summer, meltwater generated on the Greenland ice sheet surface accesses the ice sheet bed, lubricating basal motion and resulting in periods of faster ice flow. However, the net impact of varying meltwater volumes upon seasonal and annual ice flow, and thus sea level rise, remains unclear. We show that two extreme melt events (98.6% of the Greenland ice sheet surface experienced melting on July 12, the most significant melt event since 1889, and 79.2% on July 29) and summer ice sheet runoff ∼3.9σ above the 1958–2011 mean resulted in enhanced summer ice motion relative to the average melt year of 2009. However, despite record summer melting, subsequent reduced winter ice motion resulted in 6% less net annual ice motion in 2012 than in 2009. Our findings suggest that surface melt–induced acceleration of land-terminating regions of the ice sheet will remain insignificant even under extreme melting scenarios.Surface melting and runoff from the Greenland ice sheet (GrIS) has increased during the last 30 y (1–3) coincident with Northern Hemisphere warming (4, 5) resulting in unprecedented melt extents (6) and widespread dynamic thinning, which has penetrated up to 120 km into the ice sheet interior (7). One potential dynamic thinning mechanism is surface melt–induced acceleration of ice sheet motion (termed hydrodynamic coupling) during summer (8–11). Observations of GrIS ice motion during the summer show considerable variability over a range of timescales (12). Rapid variations in meltwater input from the ice sheet surface to the glacier bed result in periods when the subglacial drainage system is more highly pressurized, leading to an increase in basal sliding (13, 14). This mechanism explains both multiday increases in ice motion at the beginning of the melt season, analogous to the “spring events” observed at alpine glaciers (15), and increases in velocity at other times when meltwater is delivered to the bed at a rate faster than the subglacial drainage system can expand to accommodate.However, as the drainage-system capacity gradually expands in response to increased melting, the subglacial water pressure falls and higher velocities can therefore only be caused by much larger meltwater pulses than earlier in the melt season (16). This feedback mechanism has been invoked previously to suggest that the ice sheet could flow more slowly in a warmer year, but observations have either been limited to close to the ice sheet margin (17) or have been unable to resolve the seasonal behavior responsible for the velocity variations (18). A recent study, which incorporated seasonal ice flow and melt observations extending beyond the equilibrium line, showed that summer velocity enhancement is negated by subsequent reductions in winter flow rates (19), but the bounding conditions of the study have since been exceeded by the exceptional melting observed in 2012 (20). Moreover, the current paucity of field observations is a significant impediment to modeling the impact of coupled hydrodynamics on net ice-mass loss (21).The recent trend of warmer summers in Greenland is related to an increase in the frequency of anticyclonic conditions (22). Persistent anticyclonic conditions during summer 2012 resulted in extreme runoff volumes from the Greenland ice sheet (23), compounded by unprecedented melt extent in July 2012 associated with low-level liquid clouds (24), which led to flood damage such as the destruction of the Watson River bridge in Kangerlussuaq, west Greenland. These conditions resulted in a year during which ice sheet–wide runoff set a new record at ∼3.9σ above the 1958–2011 mean (23). National Centers for Environmental Prediction/National Center for Atmospheric Research reanalysis 1,000-mb temperature anomalies above Kangerlussuaq (25 km west of our site 1) relative to the 1981–2010 mean were +2.2 °C during May–August 2012, compared with ±0.3 °C during May–August 2009. The 2012 melt season is therefore a surrogate for potential future melting and forms a natural test for quantifying the effect of extreme meltwater supply on ice motion compared with the “average” melt year of 2009.We used global positioning system (GPS) records to observe ice motion during 2009 and 2012 at seven sites along a transect on a land-terminating margin of the GrIS, at ∼67°N (Fig. 1). Air temperature and annual ablation were also measured at each site. The lowest three sites on the transect are located within Leverett Glacier’s inferred hydrological catchment, from which we measured bulk runoff (25). Transect dynamics during 2012 (Fig. 2) had similar characteristics to previous years (10, 19): initiation of meltwater-induced acceleration during multiday “spring” speed-up events, followed by shorter duration spikes in velocity superimposed on gradually declining background seasonal velocities, which fell below premelt season velocities by the end of summer.Open in a separate windowFig. 1.Location of the transect on the western margin of the GrIS. Stars indicate sites where ice motion, temperature, and seasonal melting were measured. The triangle indicates where proglacial discharge was measured and the GPS base station located. Circles indicate locations of ‘KAN’ PROMICE/GAP weather stations along the K transect. Contours (in meters) are from a digital elevation model (DEM) of the ice sheet surface produced from interferometric synthetic aperture radar (InSAR) (29). The inferred hydrological catchment of Leverett Glacier, delineated in light gray, was calculated from the ice sheet surface DEM. Inset shows surface and bed elevation along our transect as measured by IceBridge ATM (ILATM2) and MCoRDS (IRMCR2) in 2010 and 2011, respectively (30).Open in a separate windowFig. 2.Transect observations during 2012. (A–E). Daily (24 h) along-track ice velocities (stepped black lines) and positive degree days (gray bars) for each transect site at which daily measurements were made. (F) Discharge hydrograph for Leverett Glacier (in cubic meters per second), with cumulative discharge between May 7 and August 27 (marked by gray box). The associated catchment is shown on Fig. 1. (A–F) Gray shading defines peak velocity response to July 12 and July 29 melt events (see text).Here we concentrate on two specific aspects of hydrodynamic coupling during 2012 to give insight into the likely dynamic behavior of the ice sheet in a warming climate. We examine (i) the ice flow response to the extreme melt events of July 12 and July 29 (20), and (ii) the impact of unprecedented melt volumes (23) on total annual ice motion.Enhanced ice flow lasting ∼2 d was associated with both extreme melt events (shaded periods in Figs. 2–4), with several characteristics common to both events. First, peak velocities occurred in advance of satellite-observed peak ice sheet melt extent (20), while proglacial discharge was still rising—2 d in advance of July 12, and 1 d in advance of July 29. Second, velocities increased at every site along the transect during the enhanced ice flow period. Third, at the majority of sites, velocities were lower after the enhanced ice flow period than before it (Fig. 2).Open in a separate windowFig. 4.Observations around July 29 melt event. See A–F in Fig. 3 for details.Before the July 12 melt event, sites up to 1,482 m above sea level (site 6) experienced positive air temperatures every day from June 10 (Fig. 2). Peak velocities during July 9–10 were coincident with a 2.3 °C increase in mean air temperature at our transect sites and a 73% increase in mean wind speed at Programme for Monitoring the Greenland Ice Sheet (PROMICE)/Greenland Analogue Project (GAP) K-transect sites (Fig. 1) compared with the previous 8 d. The mean daily transect velocity during July 9–10 was 61% greater than during the preceding 8 d, with sites 3 and 4 (794 and 1,061 m a.s.l, respectively) experiencing the highest peak velocities of 103 and 77% greater, respectively, than the previous 8 d (Fig. 3). By July 12, ice velocities were falling despite peaks in both ice sheet–wide melting and proglacial river discharge (∼800 m³ s⁻¹; in excess of double that observed both at the start of the melt event and in previous years; ref. 19). Sites 1 and 2 returned to daily velocities within 10 m y⁻¹ of July 1–8 mean velocities; and sites 3, 4, and 6 decreased to velocities at least 30 m y⁻¹ slower than July 1–8 mean velocities.Open in a separate windowFig. 3.Observations around July 12 melt event. (A–E) Near-surface air temperatures (dashed lines), daily (24 h) along-track ice velocities (stepped black lines) and short-term along-track ice velocities (gray lines) for each site at which daily measurements were made. Periods with inadequate quality observations removed. (F) Discharge hydrograph for Leverett Glacier (in cubic meters per second).(A–F) Gray shading defines the peak velocity response to the melt event (see text).In contrast to the July 12 melt event, a period of falling air temperatures in the previous 15 d leading up to the July 29 melt event (as low as −7 °C at site 6) resulted in falling discharge to a minimum of 240 m³ s⁻¹ on July 25, the lowest since June 18 (Fig. 2). During July 27–28, the mean transect air temperature rose by 4.4 °C compared with the previous 8 d, with associated, although lagged, increases in discharge and velocity (Fig. 4). Mean transect ice velocity on July 27–28 was 116% greater than during the preceding 8 d. At sites 2, 3, and 4, the velocity perturbation was short-lived, lasting ∼2 d before an abrupt drop in velocities that returned to within 20 m y⁻¹ of preevent velocities. Site 6 slowed down more gradually after the July 28 peak. Unlike the July 12 melt event, river discharge remained close to its event peak of ∼400 m³ s⁻¹ for 8 d following the July 29 melt event (Fig. 2F).Increased ice velocities in the lead up to the July 12 and 29 melt events were clearly caused by a rapid increase in the rate of meltwater supply to the ice sheet bed forced by changes in the rate of surface melting. Although antecedent melt conditions and the absolute volumes of meltwater associated with each event were different, the nature and style of meltwater forcing, overwhelming the capacity of the hydrological system and leading to ice acceleration, were very similar, replicating responses observed previously (12, 16).We estimated the potential contribution of each melt event to summer ice displacement by comparison with estimates of the projected ice displacement that would have occurred in the absence of the melt events. We used mean ice velocities at each site during the 8 d preceding each 2-d period of enhanced ice flow to estimate what the total displacement would likely have been through the 2-d enhanced ice flow period and the following 8 d in the absence of the enhanced ice flow period (see Materials and Methods for more information). Observations during the corresponding time periods are shown in Figs. 3 and ​and4.4. On average, the July 12 melt event forced only 7% more ice displacement over the 10-d period, and the July 29 melt event, which was preceded by lower melt rates than the July 12 event, forced 34% more ice displacement over the equivalent 10-d period. These findings reinforce the importance of antecedent melt rates (as opposed to simply meltwater volume) and thus drainage system efficiency in controlling the short-term dynamic response to variations in meltwater supply (26).The second exceptional characteristic of 2012 was ice sheet–wide runoff of ∼3.9σ above the 1958–2011 mean (23). For comparison, Fig. 5 shows observations collected along our transect in the average melt year of 2009. Exceptional melting during 2012 resulted in a mean of 117% more ablation relative to 2009 along our transect (Fig. 5A) with bulk runoff from the local ice sheet margin (2.20 × 10⁹ m³) 113% greater than 2009 (19). Summer velocities (Fig. 5B) at all but the lowest two sites were also higher in 2012 than in 2009. However, winter velocities at all sites were on average 11% lower in 2012 than in 2009, resulting in 6% less net annual ice motion along the transect in 2012–2013 than in 2009–2010 (Fig. 5). These observations support previous findings that stronger melting results in faster summers, but that faster summers are then offset by subsequent slower winter ice flow due to the evolution of a larger, more extensive subglacial drainage system that drains high basal water pressure regions (19). Our findings also support ice sheet modeling results (21), which suggest that enhanced basal lubrication will not cause substantial net mass loss from the ice sheet, and provide the observations which Shannon et al. (21) had stated were currently “insufficient to determine whether changes in subglacial hydraulics will limit the potential for the speedup of flow.”Open in a separate windowFig. 5.(A) Annual (May 1–April 30) ablation in water equivalent meters for sites 2–7 in 2009 and 2012. (B) Summer (Sum, May 1–August 31), winter (Win, September 1–April 30), and annual (Ann, May 1–April 30) velocities for each site in 2009 and 2012.Our findings demonstrate that despite the exceptional melting observed in 2012, annual ice motion along our transect was not enhanced relative to an average melt year (2009). These findings suggest that although hydrologically forced ice motion influences short-term and seasonal ice dynamics, land-terminating margins of the Greenland ice sheet are insensitive dynamically over annual timescales to melt volumes that are commensurate with temperature projections for 2100 (27). Furthermore, our data demonstrate that the importance of hydrologically forced ice motion over annual timescales can only be understood with reference to both summer and winter seasonal velocities due to their significant interannual variability. We also note that the effects of surface melt and oceanic forcing mechanisms on the dynamics of marine terminating glaciers in a warming climate remain unclear and should be a priority for future research.     

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.     

North American tree migration paced by climate in the West,lagging in the East

Tree fecundity and recruitment have not yet been quantified at scales needed to anticipate biogeographic shifts in response to climate change. By separating their responses, this study shows coherence across species and communities, offering the strongest support to date that migration is in progress with regional limitations on rates. The southeastern continent emerges as a fecundity hotspot, but it is situated south of population centers where high seed production could contribute to poleward population spread. By contrast, seedling success is highest in the West and North, serving to partially offset limited seed production near poleward frontiers. The evidence of fecundity and recruitment control on tree migration can inform conservation planning for the expected long-term disequilibrium between climate and forest distribution.

Effective planning for the redistribution of habitats from climate change will depend on understanding demographic rates that control population spread at continental scales. Mobile species are moving, some migrating poleward (1, 2) and/or upward in elevation (3, 4). Species redistribution is also predicted for sessile, long-lived trees that provide the resource and structural foundation for global forest biodiversity (5–7), but their movement is harder to study. Contemporary range shifts are recognized primarily where contractions have followed extensive die-backs (8) or where local changes occur along compact climate gradients in steep terrain (9, 10). Whether migration capacity can pace habitat shifts of hundreds of kilometers on decade time scales depends on seed production and juvenile recruitment (Fig. 1A), which have not been fitted to data in ways that can be incorporated in models to anticipate biogeographic change (11–13). For example, do the regions of rapid warming coincide with locations where species can produce abundant seed (Fig. 1B)? If so, does seed production translate to juvenile recruitment? Here, we combine continent-wide fecundity estimates from the Masting Inference and Forecasting (MASTIF) network (13) with tree inventories to identify North American hotspots for recruitment and find that species are well-positioned to track warming in the West and North, but not in parts of the East.Open in a separate windowFig. 1.Transitions, hypothesized effects on spread, and sites. (A) Population spread from trees (BA) to new recruits is controlled by fecundity (seed mass per BA) followed by recruitment (recruits per seed mass). (B) The CTH that warming has stimulated fecundity ahead of the center of adult distributions, which reflect climate changes of recent decades. Arrows indicate how centroids from trees to fecundity to recruitment could be displaced poleward with warming climate. (C) The RSH that cold-sensitive fecundity is optimal where minimum temperatures are warmer than for adult trees and, thus, may slow northward migration. The two hypotheses are not mutually exclusive. B and C refer to the probability densities of the different life stages. (D) MASTIF sites are summarized in SI Appendix, Table S2.2 by eco-regions: mixed forest (greens), montane (blues), grass/shrub/desert (browns), and taiga (blue-green).Suitable habitats for many species are projected to shift hundreds of kilometers in a matter of decades (14, 15). While climate effects on tree mortality are increasingly apparent (16–19), advances into new habitats are not (20–23). For example, natural populations of Pinus taeda may be sustained only if the Northeast can be occupied as habitats are lost in the South (Fig. 2). Current estimates of tree migration inferred from geographic comparisons of juvenile and adult trees have been inconclusive (2, 7, 21, 24, 25). Ambiguous results are to be expected if fecundity and juvenile success do not respond to change in the same ways (20, 26–29). Moreover, seedling abundances (7, 30) do not provide estimates of recruitment rates because seedlings may reside in seedling banks for decades, or they may turn over annually (31–33). Another method based on geographic shifts in population centers calculated from tree inventories (3, 34) does not separate the effects of mortality from recruitment, i.e., the balance of losses in some regions against gains in others. The example in Fig. 2 is consistent with an emerging consensus that suitable habitats are moving fast (2, 14, 15), even if population frontiers are not, highlighting the need for methods that can identify recruitment limitation on population spread. Management for forest products and conservation goals under transient conditions can benefit from an understanding of recruitment limitation that comes from seed supply, as opposed to seedling survival (35).Open in a separate windowFig. 2.Suitable habitats redistribute with decade-scale climate change for P. taeda (BA units m² /ha). (Suitability is not a prediction of abundance, but rather, it is defined for climate and habitat variables included in a model, to be modified by management and disturbance [e.g., fire]. By providing habitat suitability in units of BA, it can be related it to the observation scale for the data.) Predictive distributions for suitability under current (A) and change expected from mid-21^st-century climate scenario Representative Concentration Pathway 4.5 (B) showing habitat declines in the Southwest and East. Specific climate changes important for this example include net increases in aridity in the southeast (especially summer) and western frontier and warming to the North. Occupation of improving habitats depends on fecundity in northern parts of the range and how it is responding. Obtained with Generalized Joint Attribute Modeling (see Materials and Methods for more information).We hypothesized two ways in which fecundity and recruitment could slow or accelerate population spread. Contemporary forests were established under climates that prevailed decades to centuries ago. These climate changes combine with habitat variables to affect seeds, seedlings, and adults in different ways (36, 37). The “climate-tracking hypothesis” (CTH) proposes that, after decades of warming and changing moisture availability (Fig. 3 A and B), seed production for many species has shifted toward the northern frontiers of the range, thus primed for poleward spread. “Fecundity,” the transition from tree basal area (BA) to seed density on the landscape (Fig. 1A), is taken on a mass basis (kg/m² BA) as a more accurate index of reproductive effort than seed number (38, 39). “Recruitment,” the transition from seed density to recruit density (recruits per kg seed), may have also shifted poleward, amplifying the impact of poleward shifts in fecundity on the capacity for poleward spread (Fig. 1B). Under CTH, the centers for adult abundance, fecundity, and recruitment are ordered from south to north in Fig. 1B as might be expected if each life-history stage leads the previous stage in a poleward migration.Open in a separate windowFig. 3.Climate change and tracking. (A) Mean annual temperatures since 1990 have increased rapidly in the Southwest and much of the North. (Zero-change contour line is in red.) (B) Moisture deficit index (monthly potential evapotranspiration minus P summed over 12 mo) has increased in much of the West. (Climate sources are listed in SI Appendix.) (C) Fecundity (kg seed per BA summed over species) is high in the Southeast. (D) Recruits per kg seed (square-root transformed) is highest in the Northeast. (E and F) Geographic displacement of 81 species show transitions in Fig. 1A, as arrows from centroids for adult BA to fecundity (E) and from fecundity to recruitment (F). Blue arrows point north; red arrows point south. Consistent with the RSH (Fig. 1B), most species centered in the East and Northwest have fecundity centroids south of adult distributions (red arrows in E). Consistent with the CTH, species of the interior West have fecundity centroids northwest of adults (blue arrows). Recruitment is shifted north of fecundity for most species (blue arrows in F). SI Appendix, Fig. S2 shows that uncertainty in vectors is low.The “reproductive-sensitivity hypothesis” (RSH) proposes that recruitment may limit population growth in cold parts of the range (Fig. 1C), where fecundity and/or seedling survival is already low. Cold-sensitive reproduction in plants includes late frost that can disrupt flowering, pollination, and/or seed development, suggesting that poleward population frontiers tend to be seed-limited (40–44). While climate warming could reduce the negative impacts of low temperatures, especially at northern frontiers, these regions still experience the lowest temperatures. The view of cold-sensitive fecundity as a continuing rate-limiting step, i.e., that has not responded to warming in Fig. 1C, is intended to contrast with the case where warming has alleviated temperature limitation in Fig. 1B. Lags can result if cold-sensitive recruitment naturally limits growth at high-latitude/high-altitude population frontiers (Fig. 1C). In this case, reproductive sensitivity may delay the pace of migration to an extent that depends on fecundity, recruitment, or both at poleward frontiers. The arrows in Fig. 1C depict a case where optimal fecundity is equator-ward of optimal growth and recruitment. The precise location of recruitment relative to fecundity in Fig. 1C will depend on all of the direct and indirect effects of climate, including through seed and seedling predators and disturbances like fire. Fig. 1C depicts one of many hypothetical examples to show that climate variables might have opposing effects on fecundity and recruitment.Both CTH and RSH can apply to both temperature and moisture; the latter is here quantified as cumulative moisture deficit between potential evapotranspiration and precipitation, D=∑m=112(PETm−Pm) for month m, derived from the widely used Standardized Precipitation Evapotranspiration Index (45). Whereas latitude dominates temperature gradients and longitude is important for moisture in the East, gradients are complicated by steep terrain in the West, with temperature tending to decline and moisture increase with elevation.We quantified the transitions that control population spread, from adult trees (BA) to fecundity (seeds per BA) to recruitment (recruits per kg seed) (Fig. 1A–C). Fecundity observations are needed to establish the link between trees and recruits in the migration process. They must be available at the tree scale across the continent because seed production depends on tree species and size, local habitat, and climate for all of the dominant species and size classes (13, 46). These estimates are not sufficient in themselves, because migration depends on seed production per area, not per tree. The per-area estimates come from individual seed production and dispersal from trees on inventory plots that monitor all trees that occupy a fixed sample area. Fecundity estimates were obtained in the MASTIF project (13) from 211,000 (211K) individual trees and 2.5 million (2.5M) tree-years from 81 species. We used a model that accommodates individual tree size, species, and environment and the codependence between trees and over time (Fig. 1C). In other words, it allows valid inference on fecundity, the quasisynchronous, quasiperiodic seed production typical of many species (47). The fitted model was then used to generate a predictive distribution of fecundity for each of 7.6M trees on 170K forest inventory plots across the United States and Canada. Because trees are modeled together, we obtain fecundity estimates per plot and, thus, per area. BA (m² /ha) of adult trees and new recruits into the smallest diameter class allowed us to determine fecundity as kg seed per m² BA and recruitment per kg seed, i.e., each of the transitions in Fig. 1A.Recruitment rates, rather than juvenile abundances, come from the transitions from seedlings to sapling stages. The lag between seed production and recruitment does not allow for comparisons on an annual basis; again, residence times in a seedling bank can span decades. Instead, we focus on geographic variation in mean rates of fecundity and recruitment.We summarized the geographic distributions for each transition as 1) the mean transition rates across all species and 2) the geographic centroids (central tendency) for each species as weighted-average locations, where weights are the demographic transitions (BA to fecundity, fecundity to recruitment, and BA to recruitment). We analyzed central tendency, or centroids (e.g., refs. 3 and 34) because range limits cannot be accurately identified on the basis of small inventory plots (21). If fecundity is not limiting poleward spread (CTH of Fig. 1B), then fecundity centroids are expected to be displaced poleward from the adult population. If reproductive sensitivity dominates population spread (RSH of Fig. 1C), then fecundity and/or recruitment centroids will be displaced equator-ward from adult BA. The same comparisons between fecundity and recruitment determine the contribution of recruitment to spread.     

Human-induced ecological cascades: Extinction,restoration, and rewilding in the Galápagos highlands

Oceanic islands support unique biotas but often lack ecological redundancy, so that the removal of a species can have a large effect on the ecosystem. The larger islands of the Galápagos Archipelago once had one or two species of giant tortoise that were the dominant herbivore. Using paleoecological techniques, we investigate the ecological cascade on highland ecosystems that resulted from whalers removing many thousands of tortoises from the lowlands. We hypothesize that the seasonal migration of a now-extinct tortoise species to the highlands was curtailed by decreased intraspecific competition. We find the trajectory of plant community dynamics changed within a decade of the first whaling vessels visiting the islands. Novel communities established, with a previously uncommon shrub, Miconia, replacing other shrubs of the genera Alternanthera and Acalypha. It was, however, the introduction of cattle and horses that caused the local extirpation of plant species, with the most extreme impacts being evident after c. 1930. This modified ecology is considered the natural state of the islands and has shaped subsequent conservation policy and practice. Restoration of El Junco Crater should emphasize exclusion of livestock, rewilding with tortoises, and expanding the ongoing plantings of Miconia to also include Acalypha and Alternanthera.

Owing to their isolation, island systems often have simplified food webs, with less ecological redundancy than mainland systems (1). Islands also account for a disproportionately high number of critically endangered and recently extinct species (2), with the majority of those extinctions occurring soon after the initial colonization by humans (3). Here, we describe the cascading ecological effects of a whaling-precipitated vegetation change in the highlands of the Galápagos Islands, mediated through the depletion of tortoises and the introduction of livestock. Tortoises have repeatedly proven to be successful colonists of remote oceanic islands (e.g., Aldabra, the Mascarene Islands, the Seychelles, and the Galápagos Islands), where they evolve to become the dominant, and usually the only, megaherbivore. On the semidesert islands of the Galápagos, a classic coevolutionary adaptive arms race of leg and neck elongation between saddle-backed Chelonoidis spp. giant tortoises and trunk height in Opuntia spp. cactus highlights the struggle between predator and prey (4). But in the limited food webs on the islands, the Opuntia cactus are reliant on the tortoises for the dispersal of their seeds (5, 6). The other main morphotype of Chelonoidis tortoises, the dome-shelled tortoise, does not show such tight coevolution as they have a broader diet and are not confined to the lowlands. During the cool months of October to February, tagged tortoise telemetry shows that dome-shelled tortoises migrate up to about 420-m elevation (7). From evolutionary data, it is evident that this elevation has not been a physiological boundary as modern tortoise populations within the archipelago are resident on high calderas at elevations >1,000 m when cut off from the lowlands by lava fields (8). Thus, the height of the observed migration may reflect energetic tradeoffs that could be influenced by forage quality and competition (9). Prior paleoecological analyses from the island of Santa Cruz used the occurrence of the coprophilous fungus Sporormiella spp. to suggest that tortoises visited bogs at c. 700-m elevation until c. 700 to 500 y ago (10).The Galápagos Islands were only visited by pirates and wayward sailors prior to the 1780s (Fig. 1) (12). However, rising demand for whale oil and overexploitation of Atlantic whale stocks forced whalers to explore the Pacific Ocean (13). Whalers began operating around the Galápagos in the 1790s, and by the early 1800s they were regularly visiting the islands (13). Each whaling voyage lasted about a year, and tortoises were an essential source of fresh meat, with each whaler commonly taking aboard 200 to 300 tortoises, each weighing 150 to 200 kg (14). The tortoises could be stored in the hold of the vessel, living without food or water for up to a year (14). It is estimated that more than 200,000 tortoises were removed by an estimated 700 whaling ships that visited the islands between 1800 and 1870 (13). Although comprehensive documentation of tortoise harvesting does not exist, excerpts of 79 logbooks of whaling ships that captured ∼13,000 tortoises provide insights into the harvest between 1831 and 1870 (8). These logs record the overharvesting of tortoises, manifested in the time taken to capture 200 animals more than doubling between 1830 and 1860 (13). The weight of the tortoises and the harsh terrain meant that sailors generally hunted close to the coast. Hence, depletions were strongest close to natural harbors and water sources and on “tortoise roads”—well-traveled tracks created by the animals. The replacement of whale oil with kerosene in the 1860s led to the collapse of whaling around the Galápagos.Open in a separate windowFig. 1.Map of the Galápagos Islands showing annual precipitation (11) and elevation (Inset). The area in the Inset of western San Cristóbal is shown by the red rectangle.The slow reproduction and long lifespans of tortoises condemned some tortoise species to extinction (15). By the 1880s, scientists mounted a series of expeditions (1888, 1891, 1897, 1899, and 1905 to 1996) to capture tortoises but only found them on the islands of Isabela and Pinzón (Fig. 1) (16). Thus, by the late 1800s, giant tortoises were close to extinct on all but two islands in the archipelago. Some recovery of giant tortoises has taken place, but modern populations are thought to exist at about 10% of their former density (17).El Junco Crater Lake (hereafter El Junco) lies at 660-m elevation on the island of San Cristóbal (Fig. 2), and is the only permanent freshwater lake in the Galápagos Archipelago. A sedimentary core was recovered from the center of El Junco in 2004 and analyzed for its fossil pollen and spore content (18, 19). We use spores of the dung fungus Sporormiella (expressed as a percentage of the terrestrial pollen sum) as a proxy for megaherbivore presence (20, 21), thereby documenting megaherbivore history around the lake. The lake has also provided a detailed record of climate change (18, 22), but significant changes in temperature do not align closely either with Sporormiella abundance or vegetation change (SI Appendix). While we do not discount the possibility that recent warming may influence restoration efforts, it probably does not account for the species declines documented here.Open in a separate windowFig. 2.Data on tortoise capture from 79 whaling ship logbooks for the period from 1831 to 1870 CE (13). Tortoise harvesting effort is based on our interpretation of the ships’ logbooks to tabulate the number of tortoises captured per day of hunting. The logbook entries are not detailed enough to quantify this as tortoises per person per day, as the number of sailors in the hunting party was seldom defined. Also shown is the number of ships represented by the logbooks each year.The El Junco paleoecological record shows a reduction in Sporormiella abundance in the 1790s that is coincident with the onset of whaling activity (Fig. 3). For >1,000 y prior to the 1790s, when tortoise populations are unaffected by humans, the fossil pollen spectra of El Junco are rich in two shrubby genera, Alternanthera and Acalypha (SI Appendix), both of which are eaten by tortoises (19, 23). Between 1790 and 1927, however, Alternanthera gradually declines (Fig. 3), while Miconia increases by a factor of four between 1790 and 1900. Detrended correspondence analysis (DCA) (Fig. 3), which depicts the overall trajectory of the vegetation and had been largely stable for 1,500 y (SI Appendix), produces an inflection point at c. 1790 and a new trajectory that continued into modern times.Open in a separate windowFig. 3.Percentage occurrence of selected pollen and spores from the El Junco paleoecological record. Sporormiella percentages of the pollen sum and concentrations (spores per cubic centimeter) are shown. The axis 1 scores of the DCA of the fossil pollen flora for the period from c. 1560 to 2000 CE are plotted against time. The periods of whaling and the settlement of the San Cristóbal highlands are indicated. The background is an inferred grazing pressure around El Junco based on our data and the literature, with green representing tortoise grazing and gray that of introduced livestock. The dashed lines mark inflection in vegetation history.By the 1840s, tortoises were rare on San Cristóbal (16) and were unlikely to account for increasing amounts of Sporormiella found between 1840 and 1890. Rather, these spores probably represent the spread of introduced livestock. Donkeys introduced by whalers to help transport tortoises became feral on some islands, though no specifics are known of their introduction to San Cristóbal. The introduction of pigs and goats to San Cristóbal is thought to have been in the 1830s (24). Cattle were probably first introduced to San Cristóbal in 1842 (25), when the first settlement was made at Wreck Bay (now Puerto Baquerizo Moreno), though by 1854 there were enough feral cattle on Isabel and San Cristóbal for a small company to be formed to exploit them for tallow, but that enterprise had failed by 1860 (26). In 1869, the village of El Progreso (Fig. 1) was founded at c. 200-m elevation, about 8 km from El Junco. A zooarchaeological survey of the associated midden did not find giant tortoise (27), consistent with observations that the animals were already rare in this part of the island. A governmental visitor to the island in 1890 estimated there to be about 10,000 wild cattle on San Cristóbal. We therefore infer that the peaks of Sporormiella seen in the late 1800s come from these exotic livestock. An example of an endemic plant species that appears to have thrived under this regime of disturbance is the tree fern Cyathea weatherbyana, the spores of which reach four times the terrestrial pollen sum.A sharp increase in Sporormiella, evident at c. 1927, may reflect the arrival in 1926 and settlement of Norwegian colonists at Campo Noruega, just 2 km from El Junco (Fig. 1) (28). Only one family of Norwegians continued to inhabit the area after 1930, but the vegetation around El Junco clearly transitioned from mildly to extremely impacted at this time (Fig. 3). From the 1930s until the 1960s, San Cristóbal was the most populous island in the archipelago, but very little information on land use or herd densities exists for this period. After 1933, both Alternanthera and Acalypha gradually decline to background (<5%) levels. Cyathea and Poaceae show marked increases in abundance that peak in 1953 and 1973, respectively (Fig. 3). Eckhardt (29) cites a resident stating that hunting reduced cattle densities in the 1940s, and that open grassy areas were recolonized by shrub. In the last 30 y of the record, Psidium begins to increase from <2% in the 1970s to >12% in 1986 (Fig. 3). The Psidium pollen type includes both the endemic Psidium galapageium, a small, fairly uncommon tree of midelevations, and the exotic Psidium guajava, whose fruit (guava) is spread by cattle, horses, and tortoises (30). By outcompeting all other shrubs, P. guajava forms near-monodominant stands (31).     

From the Cover: Evidence for superconductivity in Li-decorated monolayer graphene

Monolayer graphene exhibits many spectacular electronic properties, with superconductivity being arguably the most notable exception. It was theoretically proposed that superconductivity might be induced by enhancing the electron–phonon coupling through the decoration of graphene with an alkali adatom superlattice [Profeta G, Calandra M, Mauri F (2012) Nat Phys 8(2):131–134]. Although experiments have shown an adatom-induced enhancement of the electron–phonon coupling, superconductivity has never been observed. Using angle-resolved photoemission spectroscopy (ARPES), we show that lithium deposited on graphene at low temperature strongly modifies the phonon density of states, leading to an enhancement of the electron–phonon coupling of up to λ ? 0.58. On part of the graphene-derived π^?-band Fermi surface, we then observe the opening of a Δ ? 0.9-meV temperature-dependent pairing gap. This result suggests for the first time, to our knowledge, that Li-decorated monolayer graphene is indeed superconducting, with T_c ? 5.9 K.Although not observed in pure bulk graphite, superconductivity occurs in certain graphite intercalated compounds (GICs), with T_c values of up to 11.5 K in the case of CaC₆ (1, 2). The origin of superconductivity in these materials has been identified in the enhancement of electron–phonon coupling induced by the intercalant layers (3, 4). The observation of a superconducting gap on the graphitic π^?-bands in bulk CaC₆ (5) suggests that realizing superconductivity in monolayer graphene might be a real possibility. This prospect has, indeed, attracted intense theoretical and experimental efforts (6–12). In particular, recent density functional theory calculations have suggested that, analogous to the case of intercalated bulk graphite, superconductivity can be induced in monolayer graphene through the adsorption of certain alkali metals (8).Although the Li-based GIC—bulk LiC₆—is not known to be superconducting, Li-decorated graphene emerges as a particularly interesting case with a predicted superconducting T_c of up to 8.1 K (8). The proposed mechanism for this enhancement of T_c is the removal of the confining potential of the graphite C₆ layers, which changes both the occupancy of the Li 2s band (or the ionization of the Li) and its position with respect to the graphene layer. These modifications lead to an increase of the electron–phonon coupling constant from λ = 0.33 to λ = 0.61, in going from bulk to monolayer LiC₆. It has been argued that the LiC₆ monolayer should exhibit the largest values of both λ and T_c among all alkali–metal–C₆ superlattices (8). Nevertheless, although there is thorough experimental evidence for adatom-enhanced electron–phonon coupling in graphene (7, 11, 13), superconductivity has not yet been observed in decorated monolayer graphene.Angle-resolved photoemission spectroscopy (ARPES) measurements of the electronic dispersion of pristine and Li-decorated graphene at 8 K, characterized by the distinctive Dirac cones at the corners of the hexagonal Brillouin zone (Fig. 1E), are shown in Fig. 1 A and B. Li adatoms electron-dope the graphene sheet by charge transfer doping, leading to a shift of the Dirac point to higher binding energies. As evidenced by the evolution of the graphene sheet carrier density in Fig. 1F, this trend begins to saturate after several minutes of Li deposition. Concomitantly, we observe the emergence of a new spectral weight (Fig. 1E) at the Brillouin zone center (the comparison of the Γ-point ARPES dispersion for pristine and 10-min Li-decorated graphene is shown in Fig. 1 C and D). The origin of this spectral weight is probably the Li-2s band expected for this system (8) superimposed with the folded graphene bands caused by a Li superstructure, which were observed in Li and Ca bulk GIC systems (5, 14). This spectral weight, which disappears above  ～ 50 K and is not recovered on subsequent cooling, is associated with the strong enhancement of electron–phonon coupling (discussed later, see Fig. 3 and SI Appendix).Open in a separate windowFig. 1.Charge transfer doping of graphene by lithium adatoms. Dirac-cone dispersion measured by ARPES at 8 K (A) on pristine graphene and (B) after 3 min of Li evaporation along the K-point momentum cut indicated by the white line in the Fermi surface plot in E. The Dirac cone–Fermi surface was measured at this specific K point and then replicated at the other K points by symmetry (note that high-symmetry points are here defined for the Brillouin zone of pristine graphene and not of 3×3R30° reconstructed Li-graphene, which is, instead, the notation in ref. 8). The point at which the spectroscopic gap is studied is indicated by the shaded white circle. The Dirac point, (A) already located below E_F on pristine graphene because of the charge transfer from the SiC substrate, further shifts to higher energies with (B) Li evaporation. The presence of a single well-defined Dirac cone indicates a macroscopically uniform Li-induced doping. (C) Although no bands are present at the Γ-point on pristine graphene, spectral weight is detected on 10-min Li-decorated graphene in D and E. As illustrated in the 8 K sheet carrier density plot vs. Li deposition time in F, which accounts for the filling of the π??? Fermi surface, the spectral weight at Γ is observed for charge densities n_2D ? 9 × 10¹³ cm^-2 (but completely disappears if the sample temperature is raised above  ～ 50 K and is not recovered on subsequent cooling) (SI Appendix).Open in a separate windowFig. 3.Analysis of electron–phonon coupling in Li-decorated graphene. (A) Dirac dispersion from 3-min Li-decorated graphene along the k-space cut indicated in the Fermi surface plot in E that exhibits kink anomalies caused by electron–phonon coupling (white line indicates MDC dispersion). (B–D) MDC dispersion and bare bands obtained from the self-consistent Kramers–Kronig bare-band fitting (KKBF) routine (20, 21) for several Li coverages (Methods and SI Appendix); the real part of the self-energy Σ′ is shown in Right (orange indicates Σ′ from the KKBF routine analysis, and black indicates Σ′ corresponding to the Eliashberg function presented below). (F–H) Eliashberg function α²F(ω) from the integral inversion of Σ′(ω) (22) and electron–phonon coupling constant λ = 2∫dω?α²F(ω)/ω (Methods and SI Appendix); in H, the theoretical results from ref. 8 for a LiC₆ monolayer are also shown (gray shading). (I) Experimentally determined contribution to the total electron–phonon coupling (black circles) from phonon modes in the energy ranges 100–250 meV (blue shading and white circles) and 0–100 meV (orange shading); the coupling of low-energy modes strongly increases with Li coverage.Next, we use high-resolution, low-temperature ARPES to search for the opening of a temperature-dependent pairing gap along the π^?-band Fermi surface as a direct spectroscopic signature of the realization of a superconducting state in monolayer LiC₆. To increase our experimental sensitivity, as illustrated in Fig. 2A, using the approach introduced for FeAs (17) and cuprate (18) superconductors, we perform an analysis of ARPES energy distribution curves (EDCs) integrated in dk along a 1D momentum–space cut perpendicular to the Fermi surface. This method also provides the added benefit that the integrated EDCs can be modeled in terms of a simple Dynes gap function (19) multiplied by a linear background and the Fermi–Dirac distribution, all convolved with a Gaussian resolution function (Methods and Eq. 4). As shown in Fig. 2A and especially, Fig. 2B, for data from the k-space location indicated by the white circles in Figs. 1E and ​and3E,3E, a temperature dependence characteristic of the opening of a pairing gap can be observed near E_F. The leading-edge midpoints of the Li-graphene spectra move away from E_F (Fig. 2B) in cooling from 15 to 3.5 K, which is at variance with the case of Au spectra crossing precisely at E_F according to the Fermi–Dirac distribution (Fig. 2D). Fitting these data with Eq. 4 returns a gap value of Δ = 0.9 ± 0.2 meV at 3.5 K (with Γ ? 0.09 meV). [Note that the parameter Γ in the Dynes fitting function is not treated as a free-fitting parameter, because the broadening of the coherence peaks and filling in of the gap are dominated by the experimental energy resolution. However, setting this parameter to small realistic values (Γ ～ 0.1Δ) improves the fit at the center of the gap (i.e., at E = 0 in the symmetrized data) without affecting the value of the gap itself.] Given its small value compared with the experimental resolution, the gap opening is best visualized in the symmetrized data in Fig. 2C, which minimize the effects of the Fermi function. Finally, we note that the gap appears to be anisotropic and is either absent or below our detection limit along the K ? M direction (SI Appendix, Fig. S4).Open in a separate windowFig. 2.Spectroscopic observation of a pairing gap in Li-decorated graphene. (A) Dirac dispersion from 10-min Li-decorated graphene measured at 15 and 3.5 K at the k-space location indicated by the white circles in Figs. 1E and ​and3E;3E; the temperature dependence is here evaluated for EDCs integrated in the 0.1-Å^-1 momentum region about k_F shown by the white box in Lower, with (Upper) the only changes occurring near E_F. (D) Although Au spectra cross at E_f as described by the Fermi–Dirac distribution, (B) the crossing points of the Li-graphene spectra are shifted away from E_F (cyan dashed line) because of the pull back of the leading edge at 3.5 K. A fit to the Dynes gap equation (Methods) yields a gap of Δ ? 0.9 meV at 3.5 K (and 0 meV at 15 K). The superconducting gap opening is best visualized in the symmetrized data in C [i.e., by taking I(ω) + I(?ω), which minimizes the effects of the Fermi function, even in the case of finite energy and momentum resolutions (15, 16); blue and red symbols in C represent the smoothed data, whereas the light shading gives the rmsds of the raw data]. The qualitatively similar behavior observed on polycrystalline niobium—and returning a superconducting gap Δ ? 1.4 meV—is shown in E and F.The detection of a temperature-dependent anisotropic gap at the Fermi level with a leading-edge profile described by the Dynes function—with its asymmetry about E_F and associated transfer of spectral weight to just below the gap edge—is suggestive of a superconducting pairing gap. The phenomenology would, in fact, be very different in the case of a Coulomb gap, which is typically observed in disordered semiconductors (23–25) because of the combination of disorder with long-range Coulomb interactions. A Coulomb gap would lead to a rigid shift of the spectra leading edge (isotropic in momentum) and result in a vanishing of the momentum-integrated density of states at E_F. Similarly, the observed gap is unlikely to have a charge density wave origin, because the observed gap is tied to the Fermi energy as opposed to a particular high-symmetry wavevector (the latter might occur at the M points, when graphene is doped all of the way to the Van Hove singularity, resulting in a highly nested hexagonal Fermi surface; or the K points, in the case of a 3 ×3 R 30° reconstruction, leading to a Dirac point gap). Finally, we note that these measurements do not allow us to speculate on the precise symmetry of the gap along a single Dirac cone–Fermi surface or the relative phases of the gap on the six disconnected Fermi pockets. As such, our results do not rule out any of the recent proposals for a possible unconventional superconducting order parameter in graphene (9, 26, 27).To further explore the nature of the gap observed on Li-decorated graphene (and also show our ability to resolve a gap of the order of 1 meV), in Fig. 2 E and F, we show as a benchmark comparison the analogous results from a bulk polycrystalline niobium sample—a known conventional superconductor with T_c ? 9.2?K. The Dynes fit of the integrated EDCs Fermi edge in Fig. 2E determines the gap to be Δ = 1.4 ± 0.2 meV (with Γ ? 0.14 meV), in excellent agreement with reported values (28). Although the leading-edge shift (Fig. 2E) and the dip in the symmetrized spectra (Fig. 2F) are more pronounced than for Li-graphene owing to the larger gap, the behavior is qualitatively very similar. This similarity provides additional support to the superconducting origin of the temperature-dependent gap observed in Li-decorated graphene.If the spectroscopic gap observed in Li-graphene is, indeed, a superconducting gap, the responsible mechanism may likely be electron–phonon coupling, which was predicted by the theory for monolayer Li-graphene (8) and also, seen experimentally for the bulk GIC CaC₆ (5). In direct support of this scenario, we present a detailed analysis of the graphene π^?-bands in Fig. 3, showing that the Li-induced enhancement of the electron–phonon coupling is, indeed, sufficient to stabilize a low-temperature superconducting state. Graphene doped with alkali adatoms always shows a strong kink in the π^?-band dispersion at a binding energy of about 160 meV (11). For the Li-graphene studied here, the same effect is seen in the momentum distribution curve (MDC) dispersions and the corresponding real part of the self-energy Σ′ in Fig. 3 B–D. This structure stems from the coupling to carbon in-plane (C_xy) phonons (4, 8). Despite the apparent strength of this kink, the interaction with these phonon modes contributes little to the overall coupling parameter because of their high energy (note that ω is a weighting factor in the integral calculation of λ) (Methods). As illustrated by the white circles in Fig. 3I, the contribution to λ from these high-energy (100–200 meV) modes is determined to be 0.14 ± 0.05, and it remains approximately constant for all Li coverages studied here. This value is, however, too small to stabilize a superconducting state in this system (8, 11).With increasing Li coverage and the appearance of the spectral weight at Γ, significant modifications to the low-energy part of the dispersion ( ? 100 meV) become apparent (Fig. 3 B–D). With 10 min of Li deposition (Fig. 3D), an additional kink is visible at a binding energy of ∼30 meV along with the associated peak in the real part of the self-energy Σ′. The extracted (Methods) Eliashberg functions and energy-resolved λ(ω) in Fig. 3 F–H show that, at high Li coverage, phonon modes at energies below 60 meV are coupling strongly to the graphene electronic excitations. The phonon modes in this energy range are of Li in-plane (Li_xy) and C out-of-plane (C_z) character (4, 8). This assignment is in agreement with predictions (8) as shown by the direct comparison between theory and experiment in Fig. 3H. [As for the theoretical and experimental Eliashberg functions α²F(k, ω) in Fig. 3H, the agreement may, at first glance, appear not as good as the one for λ(ω). We note, however, that, in this regard, the relevant information is in the macroscopic energy distribution of the α²F(k, ω) weight rather than in its detailed structure.] As for the total electron–phonon coupling λ for each coverage (black circles in Fig. 3I), our values measured on the π^?-band Fermi surface at an intermediate location between Γ ? K and K ? M directions (Fig. 3E) provide an effective estimate for the momentum-averaged coupling strength. [The electron–phonon coupling parameter increases monotonically along the π^?-band Fermi surface in going from the Γ ? K to the K ? M direction as observed in both decorated graphene (11) and intercalated graphite (29). Empirically, the value measured at the intermediate Fermi crossing corresponds to the momentum-averaged coupling strength along the π^?-band Fermi surface.] Remarkably, the value λ = 0.58 ± 0.05 observed at the highest Li coverage (Fig. 3I) is comparable with λ = 0.61 predicted for monolayer LiC₆ (8) as well as λ ? 0.58 observed for bulk CaC₆ (29)—it is, thus, large enough for inducing superconductivity in Li-decorated graphene. It is also significantly larger than the momentum-averaged results previously reported for both Li and Ca depositions on monolayer graphene [λ ? 0.22 and λ ? 0.28, respectively (11)]. We note that achieving such a large λ-value is critically dependent on the presence of the spectral weight observed at Γ when Li is deposited on graphene at low temperatures, presumably forming an ordered structure on the surface and not intercalating. As shown in SI Appendix, we find λ = 0.13 ± 0.05 after the same sample is annealed at 60 K for several minutes, destroying the Li order and associated Γ-spectral weight.Taken together, our ARPES study of Li-decorated monolayer graphene provides evidence for the presence of a temperature-dependent pairing gap on part of the graphene-derived π^? Fermi surface. The detailed evolution of the density of states at the gap edge as well as the phenomenology analogous to the one of known superconductors, such as Nb—as well as CaC₆ and NbSe₂, which also show a similarly anisotropic gap around the K point (30–34)—indicate that the pairing gap observed at 3.5 K in graphene is most likely associated with superconductivity. Based on the Bardeen–Cooper–Schrieffer gap equation, Δ = 3.5?k_b?T_c, we estimate a superconducting transition temperature T_c ? 5.9 K, remarkably close to the value of 8.1 K found in density functional theory calculations (8). This work constitutes the first, to our knowledge, experimental realization of superconductivity in graphene—the most prominent electronic phenomenon still missing among the remarkable properties of this single layer of carbon atoms.     

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.     

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.     

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.     

From the Cover: Montsechia,an ancient aquatic angiosperm

The early diversification of angiosperms in diverse ecological niches is poorly understood. Some have proposed an origin in a darkened forest habitat and others an open aquatic or near aquatic habitat. The research presented here centers on Montsechia vidalii, first recovered from lithographic limestone deposits in the Pyrenees of Spain more than 100 y ago. This fossil material has been poorly understood and misinterpreted in the past. Now, based upon the study of more than 1,000 carefully prepared specimens, a detailed analysis of Montsechia is presented. The morphology and anatomy of the plant, including aspects of its reproduction, suggest that Montsechia is sister to Ceratophyllum (whenever cladistic analyses are made with or without a backbone). Montsechia was an aquatic angiosperm living and reproducing below the surface of the water, similar to Ceratophyllum. Montsechia is Barremian in age, raising questions about the very early divergence of the Ceratophyllum clade compared with its position as sister to eudicots in many cladistic analyses. Lower Cretaceous aquatic angiosperms, such as Archaefructus and Montsechia, open the possibility that aquatic plants were locally common at a very early stage of angiosperm evolution and that aquatic habitats may have played a major role in the diversification of some early angiosperm lineages.When did early angiosperms begin to diversify ecologically? This question is currently unanswered. Age estimates of the divergence of crown-group angiosperms using molecular clock data vary considerably, although it is in the range of (max. 210–) often accepted, 150–140 (min. 130) million years (1–7). Parsimony reconstruction of early angiosperm habit suggests that they may have been shrubs living in “damp, dark, and disturbed” habitats (8). In contrast, many living aquatic angiosperms are basal in angiosperm phylogenies [e.g., Nymphaeales in Amborella, Nymphaeales and Illiciales, Trimeniaceae-Austrobaileya (ANITA) or Ceratophyllales with the eudicots as commonly understood]. In the fossil record, we have found an aquatic angiosperm, Montsechia vidalii (Zeiller) Teixeira, which is an atypical plant fossil found in the Barremian (130–125 million years ago) freshwater limestone in the Pyrenees and Iberian Range in Spain. Montsechia (Fig. 1) lacks roots (no proximal or adventitious roots were found in more than 1,000 shoots examined) and shows flexible axes and two types of phyllotaxy and leaf morphology. The cuticle is very thin with rare stomata. The fruit is closed with a pore near the distal tip, indehiscent, and contains one unitegmic seed developed from an orthotropous and pendent ovule (Figs. 2 and ​and3).3). Cladistic analysis of these characters places Montsechia on the stem lineage basal to extant Ceratophyllum or a clade formed by Ceratophyllum and Chloranthaceae (Fig. 4) suggesting that mesangiosperms (non-ANITA angiosperms) existed 125 million years ago, as indicated by the tricolpate pollen record. Montsechia is well-adapted to a submerged aquatic habit. Montsechia is contemporaneous with another aquatic plant fossil, Archaefructus, indicating that some of the earliest angiosperms were fully aquatic very early in their ecological diversification.Open in a separate windowFig. 1.Long- and short-leaved forms of Montsechia vidalii. (A) The long-leaved specimen shows very flexuous branches and opposite, long leaves. LH02556. (Scale bar, 10 mm.) (B) The short-leaved specimen shows regularly developed lateral branches and tiny leaf rosettes. LH07198. (Scale bar, 10 mm.)Open in a separate windowFig. 2.Fruit and seed of Montsechia vidalii. The fruit shows a small apical pore (po). The funicle (f) of the single, upside-down seed (orthotropous pendent) is attached from the hilum (h) to the placenta (pl). (Scale bar, 500 µm.)Open in a separate windowFig. 3.Reconstructions of Montsechia vidalii. (A) The long-leaved form shows the opposite leaves and branches. (B) The short-leaved form shows the alternate phyllotaxy of leaves and branches bearing pairs of ascidiate, nonornamented fruits. (C and D) The fruit shows a small apical pore and a single seed developed from an orthotropous pendent ovule. The funicle arises from the placenta (near the micropyle) to the hilum (near the pollination pore). (C) Lateral view. (D) Front view. Diagram by O. Sanisidro, B.G., and V.D.-G.Open in a separate windowFig. 4.Most parsimonious position of Montsechia in a simplified tree derived from the matrix by Endress and Doyle (26) using the J & M backbone. Taxa in blue are considered ancestrally water-related (27). Diagram by C.C. and B.G.     

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.