Arctic riparian shrub expansion indicates a shift from streams gaining water to those that lose flow

Shrub expansion has been observed across the Arctic in recent decades along with warming air temperatures, but tundra shrub expansion has been most pronounced in protected landscape positions such as floodplains, streambanks, water tracks, and gullies. Here we show through field measurements and laboratory analyses how stream hydrology, permafrost, and soil microbial communities differed between streams in late summer with and without tall shrubs. Our goal was to assess the causes and consequences of tall shrub expansion in Arctic riparian ecosystems. Our results from Toolik Alaska, show greater canopy height and density, and distinctive plant and soil microbial communities along stream sections that lose water into unfrozen ground (talik) compared to gaining sections underlain by shallow permafrost. Leaf Area Index is linearly related to the change in streamflow per unit stream length, with the densest canopies coinciding with increasingly losing stream sections. Considering climate change and the circumpolar scale of riparian shrub expansion, we suggest that permafrost thaw and the resulting talik formation and shift in streamflow regime are occurring across the Low Arctic. Shrub canopy height and density is greater along Arctic streams that lose water into unfrozen ground than along gaining streams, suggest analyses of field measurements from Toolik, Alaska.

O bservations of increased tundra shrub cover provide some of the most striking examples of ecological change in the Arctic region since the 1940s [1][2][3][4] . In North America, tundra shrub expansion has been most pronounced in protected landscape positions, such as floodplains, streambanks, water tracks, and gullies [5][6][7][8][9] , while shrub patches on exposed uplands have remained relatively stable 8,10 . Retrospective studies document at least two-dozen Arctic riparian sites with increased tall shrub cover 2,4,6 ( Fig. 1). Several highly disjunct, tall shrublands have also been identified along streams in colder tundra subzones of Arctic Canada, although trends at these sites are unknown [11][12][13][14] . Nearly a century ago, Polunin 15,16 showed that birches (Betula) in Lapland and Greenland were able to conduct water through at least a half meter of frozen soil and attained "tree" size at microsites where roots could access unfrozen materials beneath the frozen-surface layer. Accordingly, the growing observational record of tall shrub expansion across the Arctic may therefore symbolize soil thermal and hydrological change in permafrost landscapes where streams are shifting from gaining (i.e., where discharge increases downstream) to losing streams (i.e., where discharge decreases downstream due to infiltration of streamflow into the ground). Permafrost temperatures are warming 17 and cold, ice-rich near-surface permafrost is degrading 18,19 . Large rivers draining to the Arctic Ocean, which are underlain by taliks in the continuous permafrost zone 20 , have shown an increased winter discharge 21,22 , indicating an increased groundwater contribution. Nutrient export from Arctic waterways of all sizes has also increased 23,24 . Riparian tall shrub expansion, from low-to high-order streams, may therefore serve as sentinels of change across geographical scales.
Increasing summer air temperature is well documented as an overarching factor promoting shrub cover increase across the Arctic [25][26][27] . Local disturbances, such as fire and thermokarst, have also been linked to increased shrub abundance 8,28-31 against a backdrop of earlier snowmelt and a longer growing season 30,32,33 . Field studies suggest a coupling between shrub expansion and permafrost thaw 34,35 , where shrub-seedling recruitment occurs on mineral soils 30,32,36 after degradation of ice-rich permafrost and subsequent surficial disturbance. Further, shrubs are known to modify a wide range of ecosystem properties and processes, such as increasing snow depth 37 , nutrient cycling 38 , reducing albedo 39 , and altering carbon 40 and surface energy balance 41 , many of which ameliorate local conditions for shrub recruitment and thus promote further infilling of existing shrub stands (e.g., through snow trapping). Evapotranspiration losses on the other hand are the same (29 W m −2 ) between tall and low-shrub tundra 42 . Rising mean annual soil temperatures are documented following shrub expansion 43,44 , especially due to Fig. 1 Locations of documented pan-Arctic riparian tall shrub expansion and schematic representation of the relationship between hydrology, permafrost, vegetation, and soils. Riparian corridors (a) have primarily seen an expansion of tall (>1.5 m) shrubs across the Arctic region (b) in recent decades. Ice-rich permafrost thaw results in ground subsidence and exposure of mineral soil (c) that promote seed-based shrub recruitment. A schematic representation of the hydrology, permafrost, vegetation, and soil system where losing streams (discharge decrease downstream) have unlike gaining streams (discharge increase downstream), tall shrub cover, no near-surface permafrost, and an underlying talik (d). Field measurements were made near Toolik, Northern Alaska (b). Permafrost map (b) from ref. 86 .
increased snow accumulation 7,45 that, in turn, is simulated with numerical models to promote the development of taliks (thawed soil surrounded by permafrost and seasonally frozen ground) 46 .
Furthermore, warming soils and thawing permafrost increase nutrient availability and lead to shifts in soil microbial community structures [47][48][49] as well as changes in the abundances of functional genes, especially genes involved in nutrient cycling [50][51][52][53] . Soil microbes control major biogeochemical cycles and are crucial symbionts of plants and, as such, influence the establishment, productivity, and composition of plant communities 54,55 . Across the Arctic, different vegetation types have shown distinct microbial communities [56][57][58][59] , but microbial characteristics associated with riparian plant communities are relatively underexplored.
We assessed low-ordered streams that are bordered by tall shrubs (and those that are not) and identified ecohydrological linkages to subsurface biogeophysical properties and stream hydrology. We suggest that changes in riparian tall shrub extent could serve as a sentinel of climate change impacts in the Arctic, i.e., as a bioindicator of shifts in stream hydrological and biological functioning where the presence of tall shrubs indicates a losing stream underlain by a talik.

Results
Vegetation structure linked to stream-drainage type. Vegetation-species composition and structure showed distinct and consistent patterns of covariation with stream-drainage type. On losing-stream sections (i.e., where discharge decreases downstream) shrublands were characterized by tall (>1.5 m), canopy-forming feltleaf willows (Salix alaxensis), while on gaining sections, low-growing (<1-m) species such as diamondleaf willow (S. pulchra) and dwarf birch (Betula nana) were dominant. Across all 23 ecological field plots near Toolik, Alaska, the mean maximum shrub height averaged much taller on losing (307 cm, SD ± 147 cm) versus gaining-stream sections (53 cm, SD ± 28 cm). Similarly, Leaf Area index (LAI) averaged higher on losing (2.5, SD ± 0.6) versus gaining-stream sections (1.3, SD ± 0.4), and displayed a linear relationship with the percent change in discharge, %ΔQ (r 2 = 0.82, P < 0.01) (Fig. 2). Maximum canopy height had a logarithmic relationship with LAI with a positive linear relationship with LAI up to heights of <0.6 m (LAI 0.8-1.5), while the taller shrub canopies (>1.5 m) included a wider range of high LAI values (1.6-3.6) ( Supplementary Fig. 1). The highest LAI value (3.6) was collected at one of the tallest willow stands (~4.5 m).
Nonmetric multidimensional scaling (NMDS) ordination and hierarchical cluster analysis show that plant communities, including understory vegetation, are clustered by stream-section drainage type and not by a stream (Fig. 3 and Supplementary  Fig. 2). This finding is also supported by the Multiple Response Permutation Procedure (MRPP) analysis (Supplementary Table 1), which shows a significant effect by stream-section drainage type (A = 0.27, P < 0.001), but not by a stream (A = −0.01, P > 0.4), on plant functional types (Supplementary Table 1). Losing-stream sections were correlated with forb and tall shrub cover, and gaining-stream sections with Sphagnum and short shrub cover ( Fig. 3 and Supplementary Table 2).
The species composition of shrub communities encountered at the field sites resembled several riparian plant associations described by Schickhoff and others 60 on Alaska's North Slope. The tall willow communities, which were found along the losing-stream sections, corresponded to the association Epilobio-Salicetum alaxensis that is characterized by tall, canopy-forming shrublands dominated by feltleaf willow. We also encountered tall willow stands codominated by littletree willow (Salix arbusculoides), which is linked to circumneutral soils with deep active layers (1 m or more). Retrospective comparison of high-resolution imagery  showed that although tall willow communities were already present at our losing-stream sections in 1949, their abundance has increased over the 67-year record ( Supplementary Fig. 3). The low-statured communities associated with gaining-stream sections correspond to Schickhoff and others' Valeriano-Salicetum pulchrae, which is linked to acidic soils and shallow active layers.
Taliks found beneath losing-stream sections. The geophysical (very-low-frequency electromagnetic, VLF EM) survey results show taliks beneath losing-stream sections at both Aesthetic and Oksrukuyik Creeks ( Fig. 4 and Supplementary Fig. 4). The unfrozen zone above the permafrost extended up to 30 m below the losingstream sections. Three (Aesthetic Cr., ΔQ + 4.7% km −1 ) of the four cross-stream transects in gaining-stream sections showed no talik, while the fourth (Oksrukuyik Cr.) included a talik beneath the streambed (ΔQ + 1.2% km −1 ). Thaw depths exceeded 125 cm (i.e., the length of the metal probe) beneath all streams. The depth to permafrost was greater in the riparian zone compared to the surrounding upland tundra (103.1 compared to 50.7 cm).
Soil microbial community linked to stream-drainage type. Soil microbial analyses showed a clear separation of all studied communities (bacteria, fungi, and functional genes) across losing and gaining drainage types ( Fig. 3 and Supplementary Fig. 2). NMDS ordinations found two-dimensional solutions for which the two axes represented 93% (plants), 80% (fungi), 96% (bacteria), and 94% (functional gene) of the variance (Fig. 3). Axis 1 represented the transition of losing to gaining sites, and was mainly correlated with C:N, pH, gravimetric soil moisture thaw depth, Fe, and the cover of tall shrubs, low shrubs, forbs, Sphagnum, and litter. Axis 2 represented the transition between the two watersheds, Oksrukuyik and Aesthetic Creek, and was correlated mainly with pH and tall shrub cover (Supplementary Data 1). In the hierachical cluster analysis, soil fungal and bacterial communities were clearly separated between losing-and gaining-stream sections. Such clustering was less clear for the functional gene communities, in which half of the gaining sites (three out of six sites) clustered with the losing sites (Supplementary Fig. 2). Furthermore, the MRPP analysis also showed Included are the linear fit (P < 0.01, R 2 0.81) and the 95% confidence band (red shading) for seven (black dots) of the eight stream sections. The downstream-measurement location of the eighth stream section (gray square), here seen as an outlier, was within 100 m (upstream) of an aufeis field on Oksrukuyik Creek and was therefore likely to be disproportionately influenced by recent rain events as a temporary upwelling zone (referred to as "outlier" and also OLD, Supplementary Figs. 3 and 5).  Table 1). However, stream identity per se (i.e., Oksrukuyik vs. Aesthetic Creek) had no significant effect on the communities. Finally, along with losing-stream sections, we observed significantly higher soil thaw depth, bulk density, pH, and lower gravimetric soil moisture and Fe than in gaining sites (Supplementary Data 2).

Discussion
Our results show that Arctic riparian plant and soil microbial communities are linked to stream-drainage type (gaining or losing) with tall shrubs appearing exclusively along losing-stream corridors underlain by taliks. Loss of surface discharge is the result of infiltration of streamflow into large taliks (<30-m deep), adding subsurface flow and storage, which in turn can provide riparian access to water in winter and baseflow downstream beyond the rain-and snowmelt-fed seasons.
The congruent patterns between stream-drainage types for plant, fungal, bacterial, and functional gene communities coincide with a set of core soil environmental factors (thaw depth, gravimetric soil moisture, pH, Fe, and C:N). Previous studies across the Arctic support the correlation between environmental factors and plant, fungal, bacterial, and functional gene communities 56,61-64 . Our findings suggest that this set of core environmental factors might act to a certain degree in a universal fashion across the different organismal groups and furthermore shape the plant-microbial interaction and their functional potential. It is a rather striking result that all communities (plant, fungal, bacterial, and functional genes) correlated with these five environmental factors considering that the organisms are from three different kingdoms and ecologies. For the losing streams, riparian soil water mixing with talik groundwater would lead to higher pH and improved drainage would lower soil moisture and result in less anoxic conditions and warmer soil temperatures-all which would promote higher decomposition rates and lead to increased nutrient availability and lower C:N. Our study lends support to an increased thaw depth and pH, reduced gravimetric soil moisture and C:N, while soil nutrients did not show a clear distinction between the drainage types. Other studies of permafrost degradation have shown that plant communities are a stronger predictor for soil microbial communities (and their function) than soil physical and chemical factors 53 . Conversely, soil microbial communities (fungi) have been shown to affect plant community composition and biomass after permafrost thaw 65 . This has led to the suggestion that plant growth and composition could be directly linked with microbial composition through the process of coevolution 53 .
The mechanisms and processes behind the observed biogeophysical relationships build upon documented knowledge while suggesting new linkages that may challenge established interpretations of cause and effect in regard to permafrost thaw and tall shrub establishment (Fig. 5). Riparian tall shrub-seedling recruitment requires exposure of mineral soils 30,32,36 , which is facilitated by ground subsidence that is caused by ice-rich permafrost thaw (or mineral soil exposure due to fluvial erosion). Warming soils and thawing permafrost increase nutrient availability 66 and lead to shifts in soil microbial communities 67 that may occur both before and during riparian tall shrub establishment. The value of access to liquid water throughout most of the winter 16 suggests that a talik may need to form prior (or in conjunction with) tall shrub recruitment. Accordingly, rising soil temperatures documented following tall shrub expansion 43,44 , especially due to increased snow accumulation 7,45 , would therefore represent an amplification of an already-initiated soil warming. The appearance of riparian tall shrubs in early-order streams may therefore represent an easily detectable "symptom" (and a cause) of ice-rich permafrost thaw, including a shift from gaining to losing streams.
Tall shrubs create abrupt, conspicuous transitions in vegetation canopy structure on the Arctic landscape that make them ideal targets for retrospective studies using multitemporal remote sensing. Our findings, albeit from a limited number of streams in Arctic Alaska, suggest that remote sensing assessments of riparian tall shrub encroachment can be used to identify permafrost thaw,  Supplementary Fig. 4, which repeat these findings while also including a gaining stream with a talik beneath.

Fig. 5 Processes controlling riparian tall shrub recruitment and success.
A warming climate, through long-term gradual warming and seasonal extreme events (summer warmth and/or abundant snow) thaw ice-rich permafrost that, in turn, initiates a cascade of interlinked responses in the biogeophysical riparian corridor system. Accordingly, the appearance of riparian tall shrubs along early-order streams may represent an easily detectable "symptom" (and cause) of ice-rich permafrost thaw, including a shift from gaining to losing streams. COMMUNICATIONS EARTH & ENVIRONMENT | https://doi.org/10.1038/s43247-020-00050-1 ARTICLE COMMUNICATIONS EARTH & ENVIRONMENT | (2020) 1:50 | https://doi.org/10.1038/s43247-020-00050-1 | www.nature.com/commsenv talik development, and aquifer expansion. The growing literature concerning Arctic shrub encroachment suggests that some streams in the cold continuous permafrost region may have transitioned from gaining to losing discharge regimes in recent decades. The biophysical controls of increasing alkalinity 68 and baseflow 21,69 in large Arctic rivers remain unclear, but could partly be explained by the establishment and expansion of taliks beneath low-order streams. Changes in riparian tall shrub extent could therefore serve as a sentinel of climate change impacts in the Arctic, i.e., as a bioindicator of shifts in stream hydrological and biological functioning with cascading effects into downstream aquatic environments.

Methods
Study site. We undertook field studies during August 21-27, 2016, just prior to leaf-off in five subbasins (~2.5-115 km 2 ) within the Kuparuk and Sagavanirktok Rivers in the Arctic Foothills of the North Slope, Alaska, north of and near Toolik field station along the Dalton Hwy ( Supplementary Figs. 5 and 6). The drainages were located on two post-glacial surfaces: the Itkillik I (50-120 K years old), and the Sagavanirktok (125-780 K years old) 70,71 . Permafrost thickness ranges from 250 to 600 m 72 . At the Imnavait Creek basin, which is within 25 km of all study sites, mean summer (June-August) and winter (December-February) air temperatures are 8.6 and −18.3°C, respectively, with a mean annual air temperature of −7.2°C (2006-2011) 73 . Snow-water equivalent as measured from snow surveys at the end of the winter averages about 120 mm (range 69-185 mm, the year 1985-2012), while total summer precipitation ranges from 100 to 350 mm in a dry and wet year, respectively (average about 200 mm) 73 . In 2016, June through August precipitation was 197 mm with a total of 13 mm of rain, August 20 through 26 ( Supplementary Fig. 7). Snowmelt is the major hydrologic event of the year, although peak annual flows can occur after late summer storm events in smaller watersheds 73 . Aufeis (river icing) is common across the region but has been decreasing in extent 74 . Tall shrub expansion has been documented at numerous floodplains in the region with most dramatic increases in proximity to waterways, but has generally remained stable on uplands away from riparian corridors 2,8 .
Tundra vegetation in the region corresponds to bioclimate subzone E, the warmest, southernmost subzone of the Arctic tundra biome 75 . Deciduous shrubs constitute the dominant canopy-forming plant functional type across much of the landscape and are especially common in riparian zones. Widespread taxa include dwarf birch (Betula nana) and a variety of willows (Salix spp.) and ericaceous species (e.g., Vaccinium spp.). Shrubs seldom exceed 0.5-m height on upland tundra, but can reach heights of 5 m or more on larger floodplains 60 . Riparian shrubs capable of forming tall canopies (>1.5-m height) include feltleaf willow (Salix alaxensis), littletree willow (S. arbusculoides), diamondleaf willow (S. pulchra), and Siberian alder (Alnus viridis ssp. fruticosa). The dominant tall shrub present on well-developed floodplains-feltleaf willow-is an early successional species that require mineral soils for initial recruitment and is intolerant of anaerobic soils and acidic pH 36 .
Field measurements. The field effort included measurements and descriptions of stream discharge, permafrost, vegetation, soil characteristics, and soil microbial communities. A total of five different streams, including a total of eight stream sections with or without extensive thickets of tall shrubs (>1.5-m height), were measured for hydrology, permafrost, and vegetation, while only two streams (each with gaining-and losing-stream section) for soil microbial communities. The sites were located in second-to fourth-order stream networks.
Each stream section was bound by an upstream and downstream dischargemonitoring point located 400-3690 m (stream-distance) apart. The change in discharge between upstream and downstream points defined the stream section as either gaining (increasing discharge) or losing (decreasing discharge). The change in discharge was normalized to the length of the respective stream section (to % change in discharge per kilometer length of stream) by manually tracing the stream-section length using high-resolution satellite imagery. Discharge measurements for each stream section were made within 1 h of each other using a handheld flowmeter (HACH FH950, HACH Company, Loveland, CO). Two or more repeated stream-velocity and stream-depth transects, each including 15 or more point measurements, were sampled at each discharge-measurement site and averaged. Stream velocity was measured over a 40-s period and was obtained at two depths when water depth exceeded 40 cm. Two streams were selected for more detailed measurements of permafrost, vegetation, and soil microbial communities: Oksrukuyik Creek (~115 km 2 ) is on Itkillik I and an unnamed creek hereafter referred to as "Aesthetic Creek" (25 km 2 ) on the Sagavanirktok glacial surface. These two streams were selected because they both included sections with riparian low-shrub cover (upstream) and tall shrub cover (downstream) and were accessible by foot from the Dalton Hwy.
Thaw depths were measured on transects perpendicular to and across each stream section. Thaw depth (or depth to rock) was measured by inserting a 1.5-mlong metal rod to the point of resistance. The following spring (May 2017), a portable Very Low Frequency (VLF, 15-30 kHz) Electromagnetic (EM) system (GEM systems, Markham, Ontario, Canada) was used prior to snowmelt to assess frozen versus thawed ground across and along streams to a depth of about 50 m. We used an inversion code VLF2Dmf (EMTOMO) to translate the measured apparent resistivity values (ohm) to a resistivity profile of the ground. The device can be used to detect differences in electrical resistivity of the ground in 3D from 5to 10-m resolution with finer details (especially relatively conductive material such as taliks) near the ground surface (~2 m). The device records a "bubble," in which the diameter is equal to the data resolution, i.e., small conduits are not resolved but rather an average value of a volume of ground. Further, smaller conductors or taliks (surrounded by relative resistive material such as permafrost) are easier to detect than smaller permafrost or resistive areas (surrounded by relative conductive material). The literature shows cutoff values on average of about 250 Ohm m and lower for thawed soils. These values depend however strongly on the soil grain-size distribution and likely vary between 300 and 100 Ohm m 76,77 . We know that the areas outside the influence of the riparian zone have stable permafrost in almost all the observed sections. The values of the resistivity for those frozen sections are our indication of resistivity values of frozen soils (>100 Ohm m). Uncertainty in EMderived permafrost absence/presence estimates can be found in the work by Minsley et al. 77 .
Three ecological field plots were established along each gaining-and losingstream section near the discharge-monitoring sites. The plots were 20 m in length and oriented in riparian vegetation parallel to the adjacent stream section, with plot widths ranging 2-5 m depending on the width of floodplains and the size and configuration of riparian shrub patches. Ecological measurements focused on vegetation-species composition, canopy structure, and soil characteristics, including (1) visual species-cover estimation of all shrub species and dominant herbaceous and nonvascular taxa, (2) LAI, (3) maximum shrub canopy height, (4) thaw depth, and (5) soil chemistry. At each plot, we estimated the live areal cover for each species (0-100%), and independently for plant functional types (e.g., tall/ low/dwarf deciduous shrub, graminoid, and moss). For LAI, canopy height, and thaw depth, we made measurements at ten points spaced at 2-m intervals along the plot centerline. The LAI measurements were taken using an LAI-2000 Plant Canopy Analyzer (LI-COR Biosciences, Lincoln, NE). LAI is a unitless metric that represents the density of leaves, wood, and standing litter projecting above an optical sensor placed on the ground. The maximum canopy height represented the uppermost branch or leaf within 15 cm of each sampling point and was measured using a metal tape. All LAI measurements obtained within a stream section (~30) were averaged when compared to the percent change in discharge. When compared to maximum canopy height, averages for each ecological plot (five values) were used. Categorical descriptions of surficial geology and geomorphic landforms were recorded using standard classification systems developed for Arctic Alaska [78][79][80] . At each ecological field plot, we also collected a fixed-volume soil sample to determine bulk density and soil moisture.
Within each ecological field plot at the Oksrukuyik and Aesthetic creeks, we randomly chose five sites to sample soil for microbial communities and soil chemistry. At each site, we excised a 5 × 5-cm soil plug from the upper 10 cm of the soil using a sterile soil knife. The five individual plugs were combined in a clean zip-lock bag and stored on ice until we returned to the laboratory (within 8 h). In the laboratory, we homogenized and then subsampled the soils and stored them at −80°C until DNA extraction and soil analysis. Prior to DNA extraction, soils were lyophilized (Labconco FreeZone 6 Plus lyophilizer) and ground at 4°C, using 0.8cm steel beads on a Genie Vortex-2 (Scientific Industries). Visible roots were removed with forceps prior to the ball milling of the soils. Samples for soil chemistry analysis were dried at 65°C until dry and were sieved using a 2-mm sieve. Samples for bulk density and soil moisture analysis were weighed before and after drying at 105°C.
Soil laboratory analyses. The dried and sieved soils were sent to Brookside Laboratories, Inc. (New Bremen, OH) for nutrient and soil texture analysis following standard protocols (Supplementary Methods). Total soil DNA was extracted from 5 g of soil using the PowerMax Soil DNA Isolation Kit (MO BIO Laboratories, Inc., Carlsbad, CA, USA). DNA qualities were checked on a Thermo Scientific NanoDrop OneC (Fisher Scientific, Madison, WI), while DNA was quantified with Quibit 2.0 fluorometer (Invitrogen, Carlsbad, CA). Extracted DNAs were sent to the Institute for Environmental Genomics (University of Oklahoma) for sequencing of fungi and bacteria on an Illumina MiSeq Sequencer, and to assess the presence and relative abundance of functional genes with the GeoChip 5.0 functional gene microarray. The GeoChip 5.0 functional microarray contains 180,000 probes, which cover~1500 functional genes involved in geochemical cycling (N, C, S, and P), metal homeostasis, stress response, organic contaminant degradation, antibiotic resistance, secondary metabolism, virulence factors, and genes specific to viruses, fungi, and protists.
Statistical analyses. Differences in means of the various soil factors among the stream-drainages types and streams were tested with a two-way ANOVA followed by the Tukey HSD test in RStudio (R version 3.3.2). Furthermore, we ran a nonmetric multidimensional scaling (NMDS) ordination in PCORD5 to examine the relationships between plant, fungal, bacterial, and functional gene communities and soil properties (e.g., pH, gravimetric soil moisture, C:N ratio, surface organic thickness, bulk density, and nutrient concentrations) across the different streams and their drainage types. Axes in the NMDS are described in terms of the environmental factors that most influence the distribution of communities (plants, fungi, bacteria, and genes) within the ordination space. Sequence numbers were rarified by the site (stream*drainage type*stream section) and we used the abundance-based Sorensen dissimilarity index (Bray-Curtis) as the dissimilarity measure. To investigate whether plant, fungal, bacterial, and functional gene communities varied statistically from random assembly among the streams and drainage types, we applied multiple-response permutation procedures (MRPP) using the Bray-Curtis index. Furthermore, we performed hierarchical cluster analysis using the flexible beta-linkage method (beta = −0.25) to show the clustering among the different sites. Finally, we used Indicator species analysis to determine plants that were indicators for the different drainage types. MRPP, cluster analysis, and Indicator Species Analysis were performed in PCORD5.0. For more details, see Supplementary Methods. Imagery time-series analyses. We compared very-high-resolution (VHR) imagery from 1949 to 2016 to determine whether the extent of tall shrubs on the two streams with losing-stream sections has been increasing in recent decades (Supplementary Fig. 6). We also attempted this for low shrubs, but low shrubs did not produce sufficient contrasting shadow in 1949 imagery to reliably perform a retrospective comparison.