Philippa Kaur

Abulafia D (2011) The Great Sea: a human history of the Mediterranean. Oxford University Press, New York
Google Scholar 

Allen PS, Meyer SE (2002) Ecology and ecological genetics of seed dormancy in downy brome. Weed Sci 50:241–247
CAS  Google Scholar 

Balfourier F, Charmet G, Ravel C (1998) Genetic differentiation within and between natural populations of perennial and annual ryegrass (Lolium perenne and L. rigidum). Heredity 81:100–110
Google Scholar 

Blondel, J, Medail F (2009) Mediterranean biodiversity and conservation. In: Woodward JC (ed) The physical geography of the Mediterranean, Oxford University Press, New York, p 615–650

Boston ESM, Montgomery WI, Hynes R, Prodohl PA (2015) New insights on postglacial colonization in western Europe: the phylogeography of the Leisler’s bat (Nyctalus leisleri). Proc R Soc B 28:20142605. https://doi.org/10.1098/rspb.2014.2605
Article  Google Scholar 

Brown T, Jones M, Powell W, Allaby R (2009) The complex origins of domesticated crops in the Fertile Crescent. Trends Ecol Evol 24:103–109
PubMed  Google Scholar 

Comes HP (2004) The Mediterranean region—a hot spot of plant ragwort, Senecio gallicus (Asteraceae). Evolution 52:355–67
Google Scholar 

Comps B, Gomory D, Letouzey J, Thiebault B, Petit RJ (2001) Diverging trends between heterozygosity and allelic richness during postglacial colonizing in the European beech. Genetics 157:389–397
CAS  PubMed  PubMed Central  Google Scholar 

Dallman PR (1998) Plant life in the world’s Mediterranean climates: the Mediterranean Basin, South Africa, Australia, Chile, and California. Oxford University Press, New York
Google Scholar 

Dykhuizen DE, Hartl DL (1983) Functional effects of PGI allozymes in Escherichia coli. Genetics 105:1–18
CAS  PubMed  PubMed Central  Google Scholar 

Erskine W, Smartt J, Muehlbauer FJ (1994) Mimicry of lentil and the domestication of common vetch and grass pea. Econ Bot 48:326–332
Google Scholar 

Esler KJ, Jacobsen AL, Pratt RB (2018) Characteristics of Mediterranean-type ecosystems. In: Esler KJ, Jacobsen AL, Pratt RB (eds) The biology of Mediterranean-type ecosystems. Oxford University Press, New York, p 23–65
Google Scholar 

Garcia P, Vences FJ, Perez de la Vega M, Allard RW (1989) Allelic and genotypic composition of ancestral Spanish and colonial Californian gene pools of Avena barbata: evolutionary implications. Genetics 122:687–694
CAS  PubMed  PubMed Central  Google Scholar 

Greuter W (1991) Botanical diversity, endemism, rarity, and extinction in the Mediterranean area: an analysis based on the published volumes of the Med-checklist. Bot Chron 10:63–79
Google Scholar 

Hewitt GM (1999) Post-glacial re-colonization of European biota. Biol J Linn Soc 68:87–112
Google Scholar 

Hewitt GM (2000) The genetic legacy of Quarternary ice ages. Nature 405:907–913
CAS  PubMed  Google Scholar 

Hewitt GM (2011) Mediterranean peninsulas: the evolution of hotspots. In: Zachos FE, Habel JC (eds) Biodiversity hotspots, Springer, Berlin, p 123–147

Huttanus TD, Novak SJ, Mack RN (2011) Propagule pressure and introduction pathways of Bromus tectorum (Cheatgrass; Poaceae) in the central United States. Int J Plant Sci 172:783–794
Google Scholar 

Jaramillo-Correa JP, Beaulieu J, Bousquet J (2001) Contrasting evolutionary forces driving population structure at expressed sequence tag polymorphisms, allozymes and quantitative traits in white spruce. Mol Ecol 10:2729–2740
CAS  PubMed  Google Scholar 

Kahler A, Allard R (1981) Worldwide patterns of genetic variation among four esterase loci in barley (Hordeum vulgare L.). Theor Appl Genet 59:101–111
CAS  PubMed  Google Scholar 

Kahler A, Allard R, Krzakowa M, Wehrhahn C, Nevo E (1980) Associations between isozyme phenotypes and environment in the slender wild oat (Avena barbata) in Israel. Theor Appl Genet 56:31–47
CAS  PubMed  Google Scholar 

Mack RN (2011) Fifty years of ‘Waging war on cheatgrass’: research advances, while meaningful control languishes. In: Richardson DM (ed) Fifty years of invasion ecology: the legacy of Charles Elton, Wiley-Blackwell, Oxford, p 253–265

McKone MJ (1985) Reproductive biology of several bromegrasses (Bromus): breeding system, pattern of fruit maturation, and seed set. Am J Bot 72:1334–1339
Google Scholar 

Medail F, Diadema K (2009) Glacial refugia influence plant diversity patterns in the Mediterranean Basin. J Biogeogr 36:1333–1345
Google Scholar 

Medail F, Quezel P (1997) Hot-spots analysis for the conservation of plant biodiversity in the Mediterranean Basin. Ann Mo Bot Gard 84:112–127
Google Scholar 

Morgenroth U (1999) Southern Iberia and the Mediterranean trade-routes. Oxf J Archaeol 18:395–401
Google Scholar 

Myers N, Mittermeier RA, Mittermeier CG, da Fonseca GAB, Kent J (2000) Biodiversity hotspots for conservation priorities. Nature 403:853–858
CAS  PubMed  Google Scholar 

Nei M (1973) Analysis of gene diversity in subdivided populations. Proc Natl Acad Sci USA 70:3321–3323
CAS  PubMed  Google Scholar 

Nei M (1977) F-statistics and analysis of gene diversity in subdivided populations. Ann Hum Genet 41:225–233
CAS  PubMed  Google Scholar 

Nei M (1978) Estimation of average heterozygosity and genetic distance from a small number of individuals. Genetics 89:583–590
CAS  PubMed  PubMed Central  Google Scholar 

Nei M, Maruyama T, Chakraborty R (1975) The bottleneck effect and genetic variability in populations. Evolution 29:1–10
Google Scholar 

Nevo E, Baum B, Beiles A, Johnson DA (1998) Ecological correlates of RAPD DNA diversity of wild barley, Hordeum spontaneum, in the Fertile Crescent. Genet Res Crop Evol 45:151–159
Google Scholar 

Nevo E, Beiles A, Zohary D (1986) Genetic resources of wild barley in the Near East: structure, evolution and application in breeding. Biol J Linn Soc 27:355–380
Google Scholar 

Nieto Feliner G (2014) Patterns and processes in plant phylogeography in the Mediterranean Basin. A review. Perspect Plant Ecol Evol Syst 16:265–278
Google Scholar 

Norwich JJ (2006) The Middle Sea: a history of the Mediterranean. Doubleday, New York
Google Scholar 

Novak SJ, Mack RN (1993) Genetic variation in Bromus tectorum (Poaceae): comparison between native and introduced populations. Heredity 71:167–176
Google Scholar 

Novak SJ, Mack RN (2016) Mating system, introduction and genetic diversity of Bromus tectorum in North America, the most notorious product of evolution within Bromus section Genea. In: Germino M, Brown C, Chambers J (eds) Exotic annual Bromus grasses in semiarid ecosystems of the Western US: assessing causes, consequences, and management alternatives, Springer: Zurich, p 99–132

Novak SJ, Mack RN, Soltis DS (1991) Genetic variation in Bromus tectorum (Poaceae): population differentiation in its North America range. Am J Bot 78:1150–1162
Google Scholar 

Oja T (1999) Allozyme diversity and interspecific differentiation of the two diploid bromegrass species, Bromus tectorum L. and B. sterilis L. (Poaceae). Plant Biol 1:679–686
CAS  Google Scholar 

Pawlak AR, Mack RN, Busch JW, Novak SJ (2015) Invasions of Bromus tectorum (L.) into California and the American Southwest: rapid, multi-directional and genetically diverse. Biol Invasions 17:287–306
Google Scholar 

Petit RJ, Aquinagalde I, De Beaulieu BJL, Bittkau C, Brewer S, Cheddadi R et al. (2003) Glacial refugia: hotspots but not melting pots of genetic diversity. Science 300:1563–1565
CAS  PubMed  Google Scholar 

Pierson EA, Mack RN (1990) The population biology of Bromus tectorum in forests: distinguishing the opportunity for dispersal from environmental restriction. Oecologia 84:519–525
PubMed  Google Scholar 

Pinhasi R, Fort J, Ammerman AJ (2005) Tracing the origin and spread of agriculture in Europe. PLoS Biol 3:e410. https://doi.org/10.1371/journal.pbio.0030410
CAS  Article  PubMed  PubMed Central  Google Scholar 

Post GE, Dinsmore JE (1933) Flora of Syria, Palestine and Sinai: a handbook of the flowering plants and ferns, native and naturalized from the Taurus to Ras Muhammad and from the Mediterranean Sea to the Syrian desert. American Press, Beirut
Google Scholar 

Provan J, Bennett KD (2008) Phylogeographic insights into cryptic glacial refugia. Trends Ecol Evol 23:564–571
PubMed  Google Scholar 

Rowe G, Sweet M, Beebee TJC (2017) An introduction to molecular ecology, 3rd edn. Oxford University Press, Oxford
Google Scholar 

Schachner LJ, Mack RN, Novak SJ (2008) Bromus tectorum (Poaceae) in midcontinental United States: population genetic analysis of an ongoing invasion. Am J Bot 95:1584–1595
PubMed  Google Scholar 

Schmidt T, Varga Z (2012) Extra-Mediterranean refugia: the rule and not the exception? Front Zoo 9:22. http://www.frontiersinzoology.com/content/9/1/22

Semple EC (1921) Geographic factors in the ancient Mediterranean grain trade. Ann Assoc Am Geographers 11:47–74
Google Scholar 

Sirkkomaa S (1983) Calculations on the decrease of genetic variation due to the founder effect. Hereditas 99:11–20
CAS  PubMed  Google Scholar 

Soltis DE, Haufler CH, Darrow DC, Gastony GJ (1983) Starch gel electrophoresis of ferns: a compilation of grinding buffers, gel and electrode buffers, and staining schedules. Am Fern J 73:9–27
Google Scholar 

Swofford DL, Selander RB (1989) BIOSYS-1. A computer program for the analysis of allelic variation in population genetics and biochemical systematics. University of Illinois, Urbana, Illinois, Release 1.7
Google Scholar 

Taberlet P, Fumagalli L, Wust-Saucy AG, Cosson JF (1998) Comparative phylogeography and postglacial colonization routes in Europe. Mol Ecol 7:453–464
CAS  PubMed  Google Scholar 

Thompson JD (2005) Plant evolution in the Mediterranean. Oxford University Press, New York
Google Scholar 

Valdiosera CE, Garcia N, Anderung C, Dalen L, Crégut‐Bonnoure E, Kahlke R-D (2007) Staying out in the cold: glacial refugia and mitochondrial DNA phylogeography in ancient European brown bears. Mol Ecol 16:5140–5148
CAS  PubMed  Google Scholar 

Valliant MT, Mack RN, Novak SJ (2007) Introduction history and population genetics of the invasive grass Bromus tectorum (Poaceae) in Canada. Am J Bot 94:1156–1169
CAS  PubMed  Google Scholar 

Volis S, Shulgina I, Ward D, Mendlinger S (2003) Regional subdivision in wild barley allozyme variation: adaptive or neutral? J Heredity 94:341–351
CAS  Google Scholar 

Widmer A, Lexer C (2001) Glacial refugia: sanctuaries for allelic richness, but not for gene diversity. Trends Ecol Evol 16:267–269
CAS  PubMed  Google Scholar 

Willcox G, Fornite S, Herveux L (2008) Early Holocene cultivation before domestication in northern Syria. Veg Hist Archaeobotany 17:313–25
Google Scholar 

Willkomm M, Lange J (1870) Prodromus florae Hispanicae seu synopsis methodica omnium plantarum in Hispania sponte nascentium vel frequentius cultarum quae innotuerunt. Stuttgartie: Germania

Wright S (1965) The interpretation of population structure by F-statistics with special regard to systems of mating. Evolution 19:395–420
Google Scholar 

Wright S (1978) Evolution and the genetics of population. Variability within and among natural populations. University of Chicago Press, Chicago
Google Scholar 

Zeder MA (2008) Domestication and early agriculture in the Mediterranean Basin: origins, diffusion, and impact. Proc Natl Acad Sci USA 105:11597–11604
CAS  PubMed  Google Scholar 

Zohary M, Feinbrun-Dothan N (1966) Flora Palaestina. Israel Academy of Science and Humanities, Jerusalem

Zohary D, Hopf M (2000) Domestication of plants in the old world: the origin and spread of cultivated plants in West Asia, Europe, and the Nile Valley. Oxford University Press, New York
Google Scholar  More

Age-length (AL) datasets
Alaska Department of Fish & Game (ADF&G) monitors the number, body size, sex, and age of Alaska salmon harvested in a variety of fisheries and on their return breeding migration from the ocean to freshwater. Age and body length (AL) data have been collected on mature adults from commercial, subsistence, and sport harvests, escapement (spawning population) projects, and test fisheries since the early 1900’s. ADF&G data has historically been archived in regional offices; however, for this project we were able to compile all available data from across the state (Supplementary Figs. S7–S10) into a single dataset, representing over 14 million raw AL samples.
The majority of Alaska salmon fisheries target mature adults during their breeding migration into freshwater. Data from commercial harvests represent the largest proportion (57%) of measurements and are generally collected from marine waters and near river mouths. Although many Alaska salmon fishing districts are designed to operate as terminal fisheries, targeting fish destined for their river of origin, even terminal fisheries can intercept salmon returning to other Alaskan populations, and many other districts are non-terminal. Because most commercial salmon fisheries in Alaska catch a combination of fish from the target stock and intercepted fish returning to other populations, commercial samples often include a mix of fish from different populations within a river drainage and outside the drainage (e.g., Southeast Alaska troll fishery may be >80% non-local fish at times). Commercial samples from some fisheries targeting wild salmon could include a relatively low but unknown proportion of hatchery-origin salmon, which could not be excluded from our analyses without individual-level information on origin (hatchery or wild). Samples from escapement enumeration projects (sampling projects that count the number of mature adults that ‘escape’ the fishery and return to freshwater) make up the next highest proportion of AL measurements (33%). Escapement projects collect AL data from fish sampled in the freshwater environment, close to or on the spawning grounds, generally at counting towers, weirs, or fences. A variety of other sampling project types (test fishing, subsistence catch, sport catch) make up the remaining portion of these data, with no single project type representing more than 5% of the samples. ADF&G recorded the name of the sampling project, generally as the name of a given river (e.g., Fish Creek) or district (e.g., Togiak District), which we refer to as sampling locations. To ensure as much as possible that methods of data collection were consistent across locations and species, we excluded data collected from projects other than commercial harvest and escapement monitoring from statistical analyses.
Age and length (AL) measurements were collected by ADF&G personnel using standard methods56. Briefly, fish length is collected to the nearest millimeter using a measuring tape or a manual or electronic measuring board, depending on project and year. Fish age was most commonly estimated by ADF&G scientists reading growth annuli on scales57. For many AL measurements, specimen sex was also recorded, predominantly using external characteristics for sex determination. Sex determination with external characteristics in ocean-phase fish is frequently unreliable58. Because most of our data come from commercial harvests that occur in ocean-phase fish prior to the development of obvious external secondary sexual characteristics, we did not analyze the sexes separately. However, other studies examining length at age with reliable sex determination have shown similar trends in size and age for males and females33,59. As in Lewis et al.19, we assume our results reflect similar trends in male and female salmon.
To ensure data were of high quality, a number of quality assurance checks were established, and data failing those checks were excluded from analysis. These checks include ensuring that ages and lengths were within reasonable bounds for each species, that sample dates were reasonable, that data were not duplicated, and that data were all of the same length measurement type (mid-eye to fork of tail). Because mid-eye to fork length was by far the most commonly used length measurement type (85% of samples) within the data, and the vast majority of sample protocols use mid-eye to fork measurements, we assumed that observations where no length measurement type was reported (0.08% of samples) were mid-eye to fork. No other unique length measurement type accounts for more than 2% of samples. We also excluded any samples that measured fewer than ten fish for a given year/location combination. After these extensive checks, we were left with measurements on over 12.5 million individual salmon.
A wide variety of gear types were used to collect samples. The three most common gear types included gillnet, seine, and weir. Sampling methods within projects did not change systematically over time; however, for at least some projects, changes did occur, such as changes in gillnet mesh materials and sizes (for commercial harvest60) or sampling location within a watershed (for escapement projects). Some of these methodology changes are sporadically reflected in the data (e.g., mesh size), whereas others are not included and difficult to capture (e.g., weir location changes). Given the inconsistency in data and metadata associated with these fine-scale methodology changes, and the spatial and temporal scale of this dataset, changes in mesh size, gear type, or fine scale location changes (movement of a project within the same river system) were not included in our analyses.
Consistency in salmon size declines
To quantify the spatial and temporal extent of body size change, we estimated the average length of fish for each species in each sampling location and return year (the year when the fish was caught or sampled on its return migration to freshwater), which we interpret as putative biological populations (henceforth referred to as populations). For each population, we averaged these annual means to find the mean body length during a baseline period before 1990 and recent period after 2010. The pre-1990 period included all data collected before 1990, though relatively little data was available before 1980. Comparing data from two discrete time periods avoids potential edge effects that would be introduced in dividing a consecutive time series. Only populations for which we had data in both periods were included (100 sockeye, 34 Chinook, 32 chum, and 13 coho salmon populations). We established a criterion of at least 3 years of data for each population during each time period for inclusion in this analysis. Although somewhat arbitrary, we chose 1990 as the end of the early period to ensure a large number of populations had sufficient data to be included, while still being early enough to provide a meaningful baseline for comparison with current data. Because our goal was to investigate trends experienced by resource users in Alaska, we included data from some stocks that are known to capture salmon that originated from areas other than Alaska. For example, estimates for Chinook salmon from Southeast Alaska are likely influenced by the inclusion of troll-caught Chinook salmon, which are largely composed of salmon originating from British Columbia (B.C.) and the U.S. West Coast. For visualization, the results of this analysis were then scaled up to the level of the fisheries management areas established by ADF&G (Fig. 1).
To quantify and visualize continuous changes in body size across time, we fit general additive models (GAMs) to annual mean population body length for each species. To avoid convergence problems due to small sample sizes, data collected before 1975 were excluded from this analysis. In contrast to previous studies that assumed monotonic linear changes in size18,19, year was included as a nonlinear smoothed term because preliminary analyses suggested that the rate of length change varied through time. We included data from all populations for which observations from five or more years were available (276 sockeye salmon populations, 202 Chinook salmon populations, 183 chum salmon populations, 142 coho salmon populations). We knew a priori that salmon populations differ in average body size, so to preserve original units (mm) while controlling for variation in absolute body length among populations, we included two fixed factors: population and region. We assigned regions based on terrestrial biomes and the drainage areas of major watershed (shown numbered on Fig. 1, colored by ADF&G management region). Repeating these GAMs on escapement data alone provided equivalent results (Supplementary Fig. S11), which confirms that our results are not due to an artifact of sampling procedures through time.
To visualize changes in age structure and size-at-age, we fit very similar GAMs to age and length-at-age data. As above we included fixed effects for population and region, as well as a nonlinear year effect. Using the same dataset as the previously described GAMs, we used either mean freshwater age, mean saltwater age, or mean length-at-age as the response variable. For length-at-age, we separately fit GAMs for the four most common age classes in each species, except coho salmon, for which sufficient data was available for only three age classes.
To determine the extent to which patterns of body size change are consistent across space within a species, we re-fit these GAMs by replacing the main year effect by either a region-by-year or population-by-year interaction and compared model fit using AIC. These nonlinear interactions allow regions or populations to differ in their patterns of length change through time. These models are more data intensive than the previous GAMs, so we included data from all populations for which our time series consisted of any 20 or more years of data (123 sockeye salmon populations, 37 Chinook salmon populations, 38 chum salmon populations, 14 coho salmon populations).
Contributions of declining age versus growth
To partition the contribution of changes in population age structure versus size-at-age to changes in mean population length, we used the chain rule61. We used the discrete time analog of the chain rule

$${Delta}left( {xy} right) = y{Delta}x + x{Delta}y,$$
(1)

and assume that change in mean length is a function of changes in population age structure, p(a), and mean length-at-age, x(a). For each species and population, age structure in year t was calculated as the proportion of individuals in each age a. Mean length in year t is given by

$$x_t = {Sigma}_ap_tleft( a right)x_tleft( a right),$$
(2)

and the year-to-year change in length is given by

$${Delta}x_t = x_{left( {t + 1} right)} – x_t = {Sigma}_ap_tleft( a right)x_tleft( a right) + {Delta}p_tleft( a right)x_tleft( a right),$$
(3)

where

$$p_t(a) = 1/2left[ {p_{t + 1}(a) + p_t(a)} right],$$
(4)

and

$${Delta}p_t(a) = left[ {p_{t + 1}(a) – p_t(a)} right].$$
(5)

Solving these formulas year-to-year for each species in each population, we estimated the proportion of change in mean length due to changes in age structure and size-at-age. We included all populations for which we had five or more years of data (though change can only be estimated for consecutive years of data) and averaged the results across populations in each region.
Causes of age and size changes
To identify potential causes of change in salmon body size, we quantified associations with a variety of indices describing physical and biological conditions in Alaska’s freshwater and marine salmon habitats. Each candidate explanatory variable was selected based on existing biological hypotheses or inclusion in previous analyses of salmon size or population dynamics.
We considered several ocean climate indicators as potential causes of change in salmon size over time. Pacific Ocean conditions are often quantified using large-scale climate indices such as the Pacific Decadal Oscillation (PDO), El Niño Southern Oscillation (ENSO), and NPGO. These large-scale indices of ocean conditions, as proxies for climate and marine environment, have been shown to affect the survival and productivity of Pacific salmon in the North Pacific Ocean62,63. PDO, NPGO64, and MEI65,66 indices were all accessed and downloaded online (PDO, http://research.jisao.washington.edu/pdo/; NPGO, http://www.o3d.org/npgo/npgo.php, accessed 2018-02-07; MEI, https://www.esrl.noaa.gov/psd/enso/mei/, accessed 2018-02-08; MEIw, https://www.beringclimate.noaa.gov/, accessed 2018-02-08). In this analysis, winter means of NPGO and MEI were used in addition to an annual mean of MEI. Two ice cover metrics were also used to capture ocean climate conditions. Bering Sea ice cover and retreat were downloaded from https://www.beringclimate.noaa.gov/, originally derived from the National Snow and Ice Data Center data. Bering Sea ice cover index represents the winter anomaly, relative to 1981–2000 mean. Bering Sea ice retreat is an index representing number of days with ice cover after March 15.
Sea surface temperature (SST) was also explored as a potential cause of the changes in salmon size and age. SST has proven to be closely linked to salmon productivity. Mueter et al.67 found that regional-scale SST predicted survival rates better than large-scale climate indices such as the PDO. They concluded that survival rates were largely driven by environmental conditions at regional spatial scales. SST was extracted from the Extended Reconstructed Sea Surface Temperature (ERSST) version 468. To approximate SST values close to the river mouths which juvenile salmonids are most likely to experience after ocean entry, a double layer of the grid cells tracing the coastline of Alaska were extracted and the mean summer SST was calculated for each region.
Because in situ fluvial temperature measurements are sparse, both spatially and temporally, compared to the coverage of the AL dataset, air temperature was used as a proxy for temperature during the freshwater life stages. Air temperature data were extracted and sorted from remote-sensed satellite observations into multi-monthly regional means by season69.
Finally, we considered the potential for competition with other salmon to influence salmon size by including the abundances of several highly abundant salmon species as explanatory covariates. Using data compiled by Ruggerone and Irvine39, we evaluated the abundance of adult pink, chum, and sockeye salmon returning to Asia and North America as a proxy for the abundance of adult salmon of each species in the North Pacific. In addition, we also considered the more localized abundance of pink, chum, and sockeye salmon returning to Alaska, because salmon body size has been shown to vary with salmon abundance in the year of return migration in some species70 at finer spatial scales. The abundances of coho and Chinook salmon were not included, because they occur at much lower abundance than sockeye, chum, and pink salmon.
We also explored marine mammal abundances as potential predictor variables, but found that the data available precluded rigorous statistical comparison with our time series of salmon size and age structure. For example, the only estimates of orca abundance available for our study area (that from Southeast Alaska and Prince William Sound) show steady, near monotonic increases through our study period71,72. Statistically, this leads to insufficient replication and high collinearity with year effects. Although caution is warranted in interpretations of any models for which the assumptions are so obviously violated, we note that preliminary analyses including marine mammal abundance were not dramatically superior in terms of variance explained or model fit. Because of these limitations, we determined that a reliable test of the effect of marine mammal predation was not possible for Alaska.
Ultimately, we only selected covariates with an absolute correlation among covariate time series of less than 0.61. By establishing this threshold for absolute pairwise covariate correlation we sought to include only covariates for which separate associations with salmon size could be identified. The final set of covariates included in our analyses were: (1) ocean climate indicators (PDO, NPGO, MEI, winter MEI (MEIw), and Bering Sea ice cover index); (2) sea surface temperature (SST); (3) air temperature as proxy for freshwater temperature; and (4) ocean salmon abundance (abundance of Alaska sockeye, pink, and chum salmon, and North Pacific wide abundance of sockeye, pink, and chum salmon).
To test hypothesized associations between temporal trends in the average body size (length) of salmon and environmental conditions, we fit a series of Bayesian hierarchical models to data describing size trends across sampling locations for each species. Because the chain rule analysis showed that changes in age structure explained greater interannual body size variation than did changes in size-at-age, we analyzed age-aggregated mean body length. Time series, starting in 1975, of annual mean length by species for each sampling location (l) and environmental covariates were mean-variance (Z) standardized prior to model fitting. Models of the form

$$L_{i,t} = mathop {{Sigma}}limits_c ( {beta _{l,c} * X_{t – delta _{c,}c}} ) + sleft( t right) + varepsilon _{l,t},$$
(6)

were fit to each salmon species separately using Bayesian methods, where Ll,t is the standardized length at each location (l) in each return or observation year (t), βl,c are coefficients describing the effect of each covariate (c) on average length at each location, and (X_{t – delta _{c,}c}) is the standardized value of each covariate in each year. The reference year for each covariate is specified relative to the return year, or year in which salmon length compositions are observed (t), by a species and covariate-specific offset δc that associates covariate effects with the hypothesized period of interaction in each species’ life history (Supplementary Table S2). Location-specific covariate effects are structured hierarchically such that parameters describing the effect of each covariate on observed changes in average length were subject to a normally-distributed prior whose hyperparameters (group-level means and standard deviations for each covariate) were estimated directly from the data:

$$beta _{l,c} sim {mathrm{Normal}}left( {mu _c,tau _c ^{2}} right),$$
(7)

This hierarchical structure permitted us to quantify both the average (group-level) association between length observations at each sampling location (l) and hypothesized covariates (i.e., the hyperparameter μc), and the level of among-location variation in these effects (i.e., (τ_c^{2})). Prior distributions for model parameters were generally uninformative, with the exception of the prior on the group-level mean covariate effects (μc) which included a mild penalty toward zero,

$$mu _c sim {mathrm{Normal}}left( {0,1} right).$$
(8)

The prior distribution of the group-level (hyper) standard deviation of covariate effects was broad and truncated at zero,

$$tau _c sim {mathrm{Normal}}left( {0,10} right)left[ {0,} right],$$
(9)

allowing the model to freely estimate the appropriate level of among-location variability in covariate effects.
Observation error was assumed to be normally distributed εl,t ~ Normal(0, σε2), with a common observation error variance (σε2) estimated as a free parameter and subject to a broad prior distribution

$$sigma _varepsilon sim {mathrm{Normal}}left( {0,10} right)left[ {0,} right].$$
(10)

Each species-specific model also included a smoothed nonlinear year effect s(t) describing residual trends in length across time that were shared among sampling (observation) locations but were not explained by the covariates. The degree of nonlinearity for the univariate smooth s(t) quantifying the common residual trend in length is controlled by the variance term (σs) for the coefficients forming the spline73, for which a broad zero-truncated prior distribution was defined:

$$sigma _s sim {mathrm{Normal}}left( {0,10} right)left[ {0,} right].$$
(11)

Hierarchical Bayesian models describing the temporal trend in location-specific salmon length were fit using the brms package73,74 in R (R Core Team 2018), which generates posterior samples using the No U-Turn Sampler implemented in the Stan software platform75. Three independent chains were run for 20,000 iterations with a 50% burn-in and saving every tenth posterior sample, resulting in 3000 posterior samples. Convergence of all chains was diagnosed by ensuring potential scale reduction factors (R̂) for each parameter were More

1.
Flombaum P, Gallegos JL, Gordillo RA, Rincón J, Zabala LL, Jiao N, et al. Present and future global distributions of the marine Cyanobacteria Prochlrococcus and Synechococcus. Proc Natl Acad Sci. 2013;110:9824–9.
CAS  PubMed  Google Scholar 
2.
Li WKW, Url S. Primary production of prochlorophytes, cyanobacteria, and eucaryotic ultraphytoplankton: Measurements from flow cytometric sorting. Limnol Oceanogr. 1994;39:169–75.
CAS  Google Scholar 

3.
Dvořák P, Casamatta DA, Poulíčková A, Hašler P, Ondřej V, Sanges R. Synechococcus: 3 billion years of global dominance. Mol Ecol. 2014;23:5538–51.
PubMed  Google Scholar 

4.
Morel A, Ahn YH, Partensky F, Vaulot D, Claustre H. Prochlorococcus and Synechococcus: a comparative study of their optical properties in relation to their size and pigmentation. J Mar Res. 1993;51:617–49.
CAS  Google Scholar 

5.
Partensky F, Blanchot J, Vaulot D. Differential distribution and ecology of Prochlorococcus and Synechococcus in oceanic waters: a review. Bull l’Institut océanographique. 1999;19:457–75.
Google Scholar 

6.
Biller SJ, Berube PM, Lindell D, Chisholm SW. Prochlorococcus: the structure and function of collective diversity. Nat Rev Microbiol. 2015;13:13–27.
CAS  PubMed  Google Scholar 

7.
Sun Z, Blanchard JL. Strong Genome-Wide Selection Early in the Evolution of Prochlorococcus Resulted in a Reduced Genome through the Loss of a Large Number of Small Effect Genes. PLoS ONE. 2014;9:e88837.
PubMed  PubMed Central  Google Scholar 

8.
Larsson J, Nylander JAA, Bergman B. Genome fluctuations in cyanobacteria reflect evolutionary, developmental and adaptive traits. BMC Evol Biol. 2011;11:187.
PubMed  PubMed Central  Google Scholar 

9.
Dufresne A, Garczarek L, Partensky F. Accelerated evolution associated with genome reduction in a free-living prokaryote. Genome Biol. 2005;6:R14.1–R14.10.
Google Scholar 

10.
Scanlan DJ, Ostrowski M, Mazard S, Dufresne A, Garczarek L, Hess WR, et al. Ecological genomics of marine picocyanobacteria. Microb Mol Biol Rev. 2009;73:249–99.
CAS  Google Scholar 

11.
Scanlan DJ. Marine Picocyanobacteria. In: Whitton B. (eds). Ecology of Cyanobacteria II: Their Diversity in Space and Time. Springer Netherlands: Dordrecht, Netherlands 2012, pp. 503–33.

12.
Sánchez-Baracaldo P, Hayes PK, Blank CE. Morphological and habitat evolution in the Cyanobacteria using a compartmentalization approach. Geobiology. 2005;3:145–65.
Google Scholar 

13.
Wang K, Wommack KE, Chen F. Abundance and distribution of Synechococcus spp. and cyanophages in the Chesapeake Bay. Appl Environ Microbiol. 2011;77:7459–68.
CAS  PubMed  PubMed Central  Google Scholar 

14.
Callieri C. Picophytoplankton in freshwater ecosystems: the importance of small-sized phototrophs. Freshw Rev. 2008;1:1–28.
Google Scholar 

15.
Stockner JG. Phototrophic picoplankton: an overview from marine and freshwater ecosystems. Limnol Oceanogr. 1988;33:765–75.
CAS  Google Scholar 

16.
Callieri C, Stockner JG. Freshwater autotrophic picoplankton: a review. J Limnol. 2002;61:1–14.
Google Scholar 

17.
Honda D, Yokota A, Sugiyama J. Detection of seven major evolutionary lineages in cyanobacteria based on the 16S rRNA gene sequence analysis with new sequences of five marine Synechococcus strains. J Mol Evol. 1999;48:723–39.
CAS  PubMed  Google Scholar 

18.
Rippka R, Deruelles J, Waterbury JB. Generic assignments, strain histories and properties of pure cultures of cyanobacteria. J Gen Microbiol. 1979;111:1–61.
Google Scholar 

19.
Wilmotte AMR, Stam WT. Genetic relationships among cyanobacterial strains originally designated as ‘Anacystis nidulans’ and some other Synechococcus strains. J Gen Microbiol. 1984;130:2737–40.
Google Scholar 

20.
Coutinho F, Tschoeke DA, Thompson F. Comparative genomics of Synechococcus and proposal of the new genus Parasynechococcus. PeerJ. 2016;4:e1522 1–18.
Google Scholar 

21.
Robertson BR, Tezuka N, Watanabe MM. Phylogenetic analyses of Synechococcus strains (cyanobacteria) using sequences of 16S rDNA and part of the phycocyanin operon reveal multiple evolutionary lines and reflect phycobilin content. Int J Syst Evol Microbiol. 2001;51:861–71.
CAS  PubMed  Google Scholar 

22.
Zwirglmaier K, Jardillier L, Ostrowski M, Mazard S, Garczarek L, Vaulot D, et al. Global phylogeography of marine Synechococcus and Prochlorococcus reveals a distinct partitioning of lineages among oceanic biomes. Environ Microbiol. 2008;10:147–61.
PubMed  Google Scholar 

23.
Rocap G, Distel DL, Waterbury JB, Chisholm SW. Resolution of Prochlorococcus and Synechococcus ecotypes by using 16S-23S ribosomal DNA internal transcribed spacer sequences. Appl Environ Microbiol. 2002;68:1180–91.
CAS  PubMed  PubMed Central  Google Scholar 

24.
Toledo G, Palenik B. Synechococcus diversity in the California Current as seen by RNA polymerase (rpoC1) gene sequences of isolated strains. Appl Environ Microbiol. 1997;63:4298–303.
CAS  PubMed  PubMed Central  Google Scholar 

25.
Fuller NJ, Marie D, Partensky F, Vaulot D, Post AF, Scanlan DJ. Clade-specific 16S ribosomal DNA oligonucleotides reveal the predominance of a single marine Synechococcus clade throughout a stratified water column in the Red Sea. Appl Environ Microbiol. 2003;69:2430–43.
CAS  PubMed  PubMed Central  Google Scholar 

26.
Dufresne A, Ostrowski M, Scanlan DJ, Garczarek L, Mazard S, Palenik BP, et al. Unraveling the genomic mosaic of a ubiquitous genus of marine cyanobacteria. Genome Biol. 2008;9:R90.
PubMed  PubMed Central  Google Scholar 

27.
Huang S, Wilhelm SW, Harvey HR, Taylor K, Jiao N, Chen F. Novel lineages of Prochlorococcus and Synechococcus in the global oceans. ISME J. 2012;6:285–97.
CAS  PubMed  Google Scholar 

28.
Callieri C, Coci M, Corno G, Macek M, Modenutti B, Balseiro E, et al. Phylogenetic diversity of nonmarine picocyanobacteria. FEMS Microbiol Ecol. 2013;85:293–301.
CAS  PubMed  Google Scholar 

29.
Crosbie ND, Pockl M, Weisse T. Dispersal and phylogenetic diversity of nonmarine picocyanobacteria inferred from 16S rRNA gene and cpcBA-intergenic spacer sequence analyses. Appl Environ Microbiol. 2003;69:5716–21.
CAS  PubMed  PubMed Central  Google Scholar 

30.
Jasser I, Królicka A, Karnkowska-Ishikawa A. A novel phylogenetic clade of picocyanobacteria from the Mazurian lakes (Poland) reflects the early ontogeny of glacial lakes. FEMS Microbiol Ecol. 2011;75:89–98.
CAS  PubMed  Google Scholar 

31.
Ernst A, Becker S, Wollenzien UIA, Postius C. Ecosystem-dependent adaptive radiations of picocyanobacteria inferred from 16S rRNA and ITS-1 sequence analysis. Microbiology. 2003;149:217–28.
CAS  PubMed  Google Scholar 

32.
Palenik B, Brahamsha B, Larimer FW, Land M, Hauser L, Chain P, et al. The genome of a motile marine Synechococcus. Nature. 2003;424:1037–42.
CAS  PubMed  Google Scholar 

33.
Palenik B, Ren Q, Dupont CL, Myers GS, Heidelberg JF, Badger JH, et al. Genome sequence of Synechococcus CC9311: Insights into adaptation to a coastal environment. Proc Natl Acad Sci. 2006;103:13555–9.
CAS  PubMed  Google Scholar 

34.
Stuart RK, Dupont CL, Johnson DA, Paulsen IT, Palenik B. Coastal strains of marine Synechococcus species exhibit increased tolerance to copper shock and a distinctive transcriptional response relative to those of open-ocean strains. Appl Environ Microbiol. 2009;75:5047–57.
CAS  PubMed  PubMed Central  Google Scholar 

35.
Stuart RK, Brahamsha B, Busby K, Palenik B. Genomic island genes in a coastal marine Synechococcus strain confer enhanced tolerance to copper and oxidative stress. ISME J. 2013;7:1139–49.
CAS  PubMed  PubMed Central  Google Scholar 

36.
Fucich D, Marsan D, Sosa A, Chen F. Complete genome sequence of Subcluster 5.2 Synechococcus sp. strain CB0101, isolated from the Chesapeake Bay. Microbiol Resour Announc. 2019;8:6–8.
Google Scholar 

37.
Marsan D, Place A, Fucich D, Chen F. Toxin-antitoxin systems in estuarine Synechococcus strain CB0101 and their transcriptomic responses to environmental stressors. Front Microbiol. 2017;8:1–11.
Google Scholar 

38.
Page R, Peti W. Toxin-antitoxin systems in bacterial growth arrest and persistence. Nat Chem Biol. 2016;12:208–14.
CAS  PubMed  Google Scholar 

39.
Unterholzner SJ, Poppenberger B, Rozhon W. Toxin-antitoxin systems. Mob Genet Elem. 2013;3:e26219 1–13.
Google Scholar 

40.
Makarova KS, Wolf YI, Koonin EV. Comprehensive comparative-genomic analysis of type 2 toxin-antitoxin systems and related mobile stress response systems in prokaryotes. Biol Direct. 2009;4:19.
PubMed  PubMed Central  Google Scholar 

41.
Kaneko T, Nakamura Y, Sasamoto S, Watanabe A, Kohara M, Matsumoto M, et al. Structural analysis of four large plasmids harboring in a unicellular cyanobacterium, Synechocystis sp. PCC 6803. DNA Res. 2003;10:221–8.
CAS  PubMed  Google Scholar 

42.
Chen Y, Holtman CK, Magnuseon RD, Youderian PA, Golden SS. The complete sequence and functional analysis of pANL, the large plasmid of the unicellular freshwater cyanobacterium Synechococcus elongatus PCC 7942. Plasmid. 2011;23:1–7.
Google Scholar 

43.
Chen F, Wang K, Kan J, Suzuki MT, Wommack KE. Diverse and unique picocyanobacteria in Chesapeake Bay, revealed by 16S-23S rRNA internal transcribed spacer sequences. Appl Environ Microbiol. 2006;72:2239–43.
CAS  PubMed  PubMed Central  Google Scholar 

44.
Xie Y, Wei Y, Shen Y, Li X, Zhou H, Tai C, et al. TADB 2.0: An updated database of bacterial type II toxin-antitoxin loci. Nucleic Acids Res. 2018;46:D749–D753.
CAS  PubMed  Google Scholar 

45.
O’Leary NA, Wright MW, Brister JR, Ciufo S, Haddad D, McVeigh R, et al. Reference sequence (RefSeq) database at NCBI: Current status, taxonomic expansion, and functional annotation. Nucleic Acids Res. 2016;44:D733–D745.
PubMed  Google Scholar 

46.
Agarwala R, Barrett T, Beck J, Benson DA, Bollin C, Bolton E, et al. Database resources of the national center for biotechnology information. Nucleic Acids Res. 2017;45:D12–D17.
CAS  Google Scholar 

47.
Nordberg H, Cantor M, Dusheyko S, Hua S, Poliakov A, Shabalov I, et al. The genome portal of the Department of Energy Joint Genome Institute: 2014 updates. Nucleic Acids Res. 2014;42:26–31.
Google Scholar 

48.
Shao Y, Harrison EM, Bi D, Tai C, He X, Ou HY, et al. TADB: a web-based resource for Type 2 toxin-antitoxin loci in bacteria and archaea. Nucleic Acids Res. 2011;39:606–11.
Google Scholar 

49.
Marchler-Bauer A, Bo Y, Han L, He J, Lanczycki CJ, Lu S, et al. CDD/SPARCLE: Functional classification of proteins via subfamily domain architectures. Nucleic Acids Res. 2017;45:D200–D203.
CAS  PubMed  Google Scholar 

50.
R Core Team. R: A Language and Environment for Statistical Computing. 2018. Vienna, Austria.

51.
Wickham H. ggplot2: elegant graphics for data analysis. New York: Springer-Verlag; 2009.
Google Scholar 

52.
Bertelli C, Laird MR, Williams KP, Lau BY, Hoad G, Winsor GL, et al. IslandViewer 4: Expanded prediction of genomic islands for larger-scale datasets. Nucleic Acids Res. 2017;45:W30–W35.
CAS  PubMed  PubMed Central  Google Scholar 

53.
Winther KS, Gerdes K. Enteric virulence associated protein VapC inhibits translation by cleavage of initiator tRNA. Prok Natl Acad Sci. 2011;108:7403–7.
CAS  Google Scholar 

54.
Harms A, Brodersen DE, Matarai N, Gerdes K. Toxins,targets,and triggers: an overview of toxin-antitoxin biology. Mol Cell. 2018;70:768–84.
CAS  PubMed  Google Scholar 

55.
Kopfmann S, Roesch S, Hess W. Type II toxin–antitoxin systems in the unicellular cyanobacterium Synechocystis sp. PCC 6803. Toxins. 2016;8:228.
PubMed Central  Google Scholar 

56.
Leplae R, Geeraerts D, Hallez R, Guglielmini J, Drze P, Van Melderen L. Diversity of bacterial type II toxin-antitoxin systems: a comprehensive search and functional analysis of novel families. Nucleic Acids Res. 2011;39:5513–25.
CAS  PubMed  PubMed Central  Google Scholar 

57.
Xia K, Bao H, Zhang F, Linhardt RJ, Liang X. Characterization and comparative analysis of toxin – antitoxin systems in Acetobacter pasteurianus. J Ind Microbiol Biotechnol. 2019;46:869–82.
CAS  PubMed  Google Scholar 

58.
Biller SJ, Berube PM, Berta-Thompson JW, Kelly L, Roggensack SE, Awad L, et al. Genomes of diverse isolates of the marine cyanobacterium Prochlorococcus. Sci Data. 2014;1:1–11.
Google Scholar 

59.
Palenik B, Barahamsha B, Larimer FW, Land M, Hauser L, Chain P, et al. The genome of a motile marine Synechococcus. Nature. 2003;424:1037–42.
CAS  PubMed  Google Scholar 

60.
Scanlan DJ, West NJ. Molecular ecology of the marine cyanobacterial genera Prochlorococcus and Synechococcus. FEMS Microbiol Ecol. 2002;40:1–12.
CAS  PubMed  Google Scholar 

61.
McWilliam H, Li W, Uludag M, Squizzato S, Park YM, Buso N, et al. Analysis tool web services from the EMBL-EBI. Nucleic Acids Res. 2013;41:W597–W600.
PubMed  PubMed Central  Google Scholar 

62.
Sevin EW, Barloy-Hubler F. RASTA-bacteria: a web-based tool for identifying toxin-antitoxin loci in prokaryotes. Genome Biol. 2007;8:R155.1–R155.14.
Google Scholar 

63.
Robson J, McKenzie JL, Cursons R, Cook GM, Arcus VL. The vapBC operon from mycobacterium smegmatis Is an autoregulated toxin-antitoxin module that controls growth via inhibition of translation. J Mol Biol. 2009;390:353–67.
CAS  PubMed  Google Scholar  More

Isolation and identification of four airborne microalgal strains belonging to three genera from Dongsha Island in the South China Sea
Four microalgal strains belonging to three genera were successfully isolated in the collection trip in December 2013 and all the four could not grow in 2f seawater medium10 (in both 33 and 22‰ salt levels), indicating they were freshwater species. As shown in Fig. 1, based on the morphology revealed by Scanning Electron Microscopy and the phylogenetic analysis, the four strains were named Scenedesmus sp. DSA1 (stands for Dongsha airborne #1), Coelastrella sp. DSA2, Coelastrella sp. DSA3, and Desmodesmus sp. DSA6 (hereafter referred to as DSA1, DSA2, DSA3 and DSA6, respectively). All four strains are members of the family Scenedesmaceae in Chlorophyta. Sequences of the 18S rDNA and ITS1-5.8S-ITS2 fragments of the four airborne microalgal strains were deposited into GenBank under the accession numbers KX818834–KX818841.
Figure 1

A Scanning electron microscopy (SEM) micrographs of the four airborne microalgae in the early stationary phase. Upper panel, left and right: Scenedesmus sp. DSA1 and Coelastrella sp. DSA2. Lower panel, left and right: Coelastrella sp. DSA3 and Desmodesmus sp. DSA6. B Phylogenetic analysis of the four airborne microalgal strains (underlined) and their related species based on the fusion sequences of ITS1 and ITS2 of each species. The numbers in the parentheses are accession numbers of each sequence in GenBank. The numbers at the nodes indicate bootstrap values (expressed as percentage) with 500 replicates.

Full size image

The airborne green microalgae had better UV tolerance than the waterborne green microalgae
When traveling in the air, airborne microalgae are inevitably exposed to much higher levels of UV radiation compared to microalgae in the aquatic environments. In order to survive aerial travel, these cells must have a mechanism(s) to protect themselves against the damaging effects of the radiation, a mechanism that is of less concern to waterborne microalgal populations. To verify whether the four airborne microalgae resisted UV radiation better than waterborne microalgae, these microorganisms were spread onto agar plates using the top agar method11, and exposed to UV-B radiation to determine their survival rates. The survival rates were compared to those of the two waterborne microalgae, Desmodesmus sp. F512 and Neodesmus sp. UTEX 2219-413 (hereafter referred to as F5 and 2219-4, respectively), both are members of Scenedesmaceae as well. As shown in Fig. 2A, the survival rates of the two waterborne microalgae were about 68 and 46%, respectively, after 1 min of the UV-B exposure, compared to more than 90% for all four airborne microalgae. After 3 min of the exposure, however, all of the waterborne cells died but the survival rates of the airborne cells were about 98, 73, 26, and 11% for DSA3, DSA2, DSA1 and DSA6, respectively. Significant differences could still be observed among the survival rates of DSA3, DSA2, and the other four species after 5 min of the UV-B exposure.
Figure 2

A UV-B stress tolerance of the four airborne microalgae and the two waterborne microalgae used as controls. About 500 cells in the early stationary phase were spread onto each agar plate and irradiated with 302 nm UV for the specified durations. The survival rate of each strain was defined as the colony numbers on the UV-treated plates compared to those on the non-treated plates (n = 6, mean ± SE). F5, Desmodesmus sp. F5; 2219-4, Neodesmus sp. UTEX 2219-4. B Different autofluorescence intensities from the cell wall of the six microalgal strains.

Full size image

For algal cells, cell wall is the first barrier to defend the cells against UV attack. It is well known that the cell walls of terrestrial plants are autofluorescent when excited by UV14. In this context, the four airborne microalgal strains would have better UV tolerance if their cell walls could absorb UV, and therefore are autofluorescent as well. To examine this possibility, the pigment-free cells were examined under an epifluorescence microscope. Indeed, as shown in Fig. 2B, the cell walls of the four airborne and the two waterborne green microalgae were able to emit autofluorescence when excited by UV under the microscope. Furthermore, the autofluorescence intensities emitted from the six microalgal strains varied, with DSA3 and DSA2 being the strongest, DSA1 and DSA6 in the middle, and F5 and 2219-4 the weakest as detected using the same parameters. The varied intensities suggested that these microalgae had different levels of UV tolerance because their cell walls absorbed different amounts of UV energy.
The UV tolerance of the microalgae was positively correlated with their cell wall thickness
There appeared to be a correlation between the autofluorescence intensities emitted by the cell walls and the survival rates of the six microalgal strains exposed to the UV-B radiation for 3 min. It was intriguing to investigate whether the cell wall thickness of these microalgae played a role in the survival rates. To measure the cell wall thickness of the six microalgal strains, these cells were fixed and sectioned, and then observed using Transmission Electron Microscopy (TEM). As shown in Fig. 3A, the cell wall thickness of the six strains varied, with DSA3 and DSA2 being the thickest, DSA1 and DSA6 intermediate, and F5 and 2219-4 the thinnest. This relation well-correlated with that of the autofluorescence intensities from their cell walls. When the survival rates under 3 min of the UV-B exposure were plotted against the cell wall thickness of the six strains, a good correlation (r = 0.99, p  More

Plant and fungal materials, maintenance, and preparation
The winter wheat cultivar Zheng 9023 and the spring wheat cultivar Yongliang 4 were purchased from the commercial seed market in Wuhan City and Minqin County in China, respectively. The barley cultivar Huadamai 14 was provided from Prof. Dongfa Sun in Huazhong Agricultural University, Wuhan, China. The oat cultivar Mengmai 2 was donated by Prof. Jun Zhao in Neimeng Agricultural University. The maize cultivar Zhengdan 958 was purchased from the commercial seed market in Wuhan City. The rice cultivar LTH was donated by Prof. Youliang Peng of China Agricultural University. All seeds were surface-sterilized with a 0.5% sodium hypochlorite solution (NaClO) before sowing or S. sclerotiorum treatment.
S. sclerotiorum strain DT-8, which was originally isolated from a sclerotium collected from a diseased rapeseed, is a hypovirulent strain infected by a DNA mycovirus. Strain DT-8VF, a virulent strain, is a virus-free derivative of DT-8 [21]. Strain DT-8VFRFP is a derivative of DT-8VF labeled with the mCherry fluorophore by DNA transformation; it shows normal virulence. The wheat Fusarium head blight (FHB) pathogen, F. graminearum strain PH-1, was used to inoculate wheat spikes. An S. sclerotiorum virulent strain Ep-1PNA367 and three hypovirulent strains, AH98, SCH941, and T1-1-20, were also used to investigate their potential endophytic growth on wheat. AH98 is infected by a negative-stranded RNA virus [22], while SCH941 and T1-1-20 are infected by various other mycoviruses. All S. sclerotiorum strains and F. graminearum strain PH-1 were grown on potato dextrose agar (PDA) or potato dextrose broth (PDB) at 20 °C; Magnaporthe oryzae strain 131 was grown on PDA at 28 °C, or grown on tomato–oat medium to produce conidia, and was stored on PDA slants at 4 °C.
Microscopic observation
To observe the growth of S. sclerotiorum in wheat by confocal microscopy, seeds were surface-sterilized with NaClO and sown on half-strength Murashige and Skoog (MS) agarose medium amended with 25 mM sucrose for 8 days, then root crowns were inoculated with mycelia of the strain DT-8VFRFP. Wheat seedlings were maintained at 20 °C for 4 days under a 12-h photoperiod. Seedlings were subsequently washed three times with PBS for 10 min each, and then root sections were separated and incubated in a 1:100 dilution of wheat germ agglutinin conjugated to FITC (Sigma) for 1 h at room temperature according to the manufacturer’s instructions. The roots were visualized with a LEICA confocal microscope (LEICA SP8) using the 488-nm line of a 25-mW Argon ion laser for FITC and the 561-nm line of a 20-mW solid-state laser for mCherry.
For further observation with a transmission electron microscopy (TEM), 5-mm root segments from wheat seedlings grown on MS medium for 15 days after inoculation with strain DT-8VFRFP were fixed in 0.4% (v/v) glutaraldehyde solution overnight at 4 °C. After washing in PBS buffer, roots were dehydrated with a graded ethanol series. Samples were then embedded in Epon-821 and polymerized at 60 °C. Thin sections (50 nm) were cut using a Leica ULTRACUT UCT ultramicrotome with a diamond knife.
For TEM immunodetection, wheat seedling roots were fixed in 4% (w/v) freshly depolymerized paraformaldehyde and 0.4% (v/v) glutaraldehyde in 1× PBS, pH 7.4, for 1 h at 4 °C. The samples were then embedded using an LR white embedding kit (Fluka) and polymerized at 50 °C for 24 h. Immunogold labeling specificity was detected by displacing the anti-mCherry antibodies with rabbit preimmune serum. The method for TEM immunodetection was performed as previously described [23].
For scanning electron microscopy (SEM) observation, 2-mm root segments from the wheat, barley, oat, maize, and rice seedlings treated with different strains of S. sclerotiorum grown on MS medium for 15 days were used. All segments were fixed in a 0.4% (v/v) glutaraldehyde solution overnight at 4 °C. For SEM analysis, the sections were allowed to air-dry overnight in a desiccator at room temperature, sputter-coated with gold, and prepared for SEM analysis (EVO MA 10 Carl Zeiss SMT AG, Germany). Root segments from nontreated wheat were sampled and observed as a control comparison.
For confocal microscopy observations, 45 days after the root crown inoculation with mycelia of the DT-8VFRFP strain, stems from wheat plants growing in soil in the greenhouse were carefully washed with distilled water and embedded in Tissue-Tek O.C.T. compound medium (Sakura Finetek USA, Inc., Torrance, CA) at −23 °C overnight. Microtome sections (25 µm thick) were sliced using a freezing microtome (LEICA SP8, Germany). The stem microtome sections were then visualized with a LEICA confocal microscope using the 561-nm excitation wavelength for mCherry.
Re-isolation of S. sclerotiorum
To further probe whether S. sclerotiorum can go to aerial parts of wheat when inoculated on root of wheat seedling, plants were grown in sterile nutrient soil for 45 days. The samples for re-isolation were taken from the second segment near the base of wheat stem of eleven individual DT-8VFRFP-treated wheat plants, and were cut into 5 mm long segments, then surface-sterilized by dipping them in 70% EtOH for 2 min and then in 0.5% NaClO for 2 min, followed by rinsing three times with sterile distilled water. The sterilized stem segments were placed on hygromycin-amended PDA medium plates (50 µg/mL) and incubated at 20 °C for 8 days. Then, the emerging colonies were identified as S. sclerotiorum based on colony morphology, and PCR amplification [24].
PCR determination of S. sclerotiorum and mycoviruses
DNA samples of either wheat or fungi were extracted using a cetyltrimethylammonium bromide method. A primer pair XJJ21/XJJ222 (GTTGCTTTGGCGTGCTGCTC/CTGACATGGACTCAATACCAATCTG) was used for detection of S. sclerotiorum [24], and a primer pair pRep-F/pRep-R (GTCACCACCCAAACATTACAAAGAGCGTATTCC/ACGTCA GGTGC) was used for detection of viral DNA of SsHADV-1. A procedure described by Yu et al. [21] was used for PCR amplification.
Seed treatment of wheat with S. sclerotiorum and sowing
To determine whether S. sclerotiorum could promote growth and disease resistance of wheat in the greenhouse, wheat seeds were washed with tap water and surface-sterilized with 0.5% NaClO for 10 min, then washed three times with sterile water. Surface-sterilized seeds were soaked in sterile water for 4 h, and then collected and blotted dry. Meanwhile, fresh mycelium of strains DT-8 and DT-8VF were collected from PDA medium and then cultured in PDB medium in 250-mL flasks in a shaker at 20 °C for 4 days to obtain hyphal fragment suspensions; the number of viable fungal fragments was adjusted to 1.4 × 105 cfu/mL before inoculation of wheat seeds. The hyphal fragment suspensions were then used to inoculate the prepared seeds (100 mL of hyphal fragment suspension/kg wheat seeds) by thoroughly mixing the hyphal fragment suspension and wheat seeds for 6 h. The inoculated seeds were further dried using an electric fan for 12 h at room temperature. Wheat seeds soaked only in sterile water for 6 h and then dried with the same method were used as a control.
To test whether the treated seeds were successfully colonized by S. sclerotiorum, S. sclerotiorum-treated seeds were surface-sterilized with 0.5% NaClO for 10 min, rinsed three times with sterile distilled water, and then cut in half and placed on PDA plates and incubated at 20 °C for 7 days. All emerging colonies from S. sclerotiorum-treated seeds were confirmed as S. sclerotiorum based on colony morphology and PCR amplification. In addition, a small number of S. sclerotiorum-treated seeds were randomly picked and sown into soil taken from the field and grown in the greenhouse. The seedlings were then tested after 21 days for the presence of S. sclerotiorum by PCR amplification; all seedlings tested positive. Hence, the seed treatment was confirmed as an efficient method for inoculation of wheat seeds with S. sclerotiorum. Consequently, we used this method to treat wheat seeds for the rest of the study, and kept treated seeds under dry conditions at room temperature for up to 7 days before sowing.
Seeds treated either with strains DT-8VF or DT-8 were sown in pots in a greenhouse. In a laboratory test, treated seeds were allowed to germinate in a Petri dish on a layer of wet filter paper. Germinating seeds were then sown either into soil that was taken from the field or sterile nutrient soil. Twelve plants were grown in each pot. Wheat seedlings were maintained under greenhouse conditions at 20 °C. Strain DT-8-treated seeds were also tested in the field. Those seeds were sown as rows in the wheat field in exactly the same way as farmers normally do in each sowing season. Field management was conducted as per normal farmer practice, except that no fungicide was applied.
In order to investigate whether S. sclerotiorum colonization can spread to aerial parts of wheat plants originating from DT-8-treated seeds, ten plants from each plot were randomly sampled from the field. A total of 30 wheat plants from the strain DT-8-treated group and 20 plants from the nontreated group were sampled. First, the roots were rinsed thoroughly with tap water. Then the roots, flag leaves, and spikes were given three brief rinses in distilled water. Each wheat plant was given a number, from 1 to 30 for DT-8-inoculated plants, and from 1 to 20 for nontreated plants. DNA samples were extracted individually from the root, flag leaf, and spike of each plant and used to determine the presence of S. sclerotiorum and mycovirus.
Evaluating the growth of treated wheat in the greenhouse and field
To evaluate the growth of S. sclerotiorum-treated wheat in greenhouse experiments, plant height was measured at the seedling and anthesis stages, while flag leaf and spike lengths were evaluated only at the anthesis stage. There were 60 plants in each group treated with strain DT-8 and strain DT-8VF and 60 plants in the nontreated group. Determination of 1000-grain weight was repeated, four in each group. Measurement data for each group were calculated for statistical analysis.
To investigate whether other strains could promote wheat growth, strains Ep-1PNA367, AH98, SCH941, and T1-1-20 were used instead of strains DT-8 and DT-8VF. Wheat seeds were treated and then sown in a sterile mixture of vermiculite and perlite at the ratio of 3:1 in pots and placed in greenhouse, with ~50 seedlings in each pot. Seedling shoot fresh weight was measured at 25 days after planting. There were 30 plants in each treatment and control, and the average weight of ten seedlings was calculated.
For field tests, plant height and the length, width, and thickness of flag leaves at the early flowering stage were measured in the field at EZhou. Forty plants from each plot were randomly measured from a total of 120 plants in the DT-8-treated group and from 120 plants in the nontreated group. Measurement data for each group were calculated for statistical analysis. All the results were confirmed with independent lines and over two planting seasons.
Field experiments and wheat yield tests
To examine whether S. sclerotiorum treatment could enhance wheat yield under natural field conditions, DT-8-treated seeds were sown in a wheat field located at EZhou in late October of 2016 and harvested in mid-May of 2017. This experiment was repeated at EZhou, Jingmen City, and Xiangyang City in late October of 2017 and 2018 and harvested in mid-May of 2018 and 2019. Furthermore, seeds of the spring wheat cultivar Yongliang 4 were also treated with strain DT-8 and sown in mid-April of 2017 at Minqin and Tianzhu Counties in Gansu province and harvested in late July of 2017. All wheat was managed as per normal farmer practice, except that no fungicide was applied. The treatments with or without strain DT-8 were replicated four times at Tianzhu County and Jingmen City in 2017 and five times at other places and the wheat yield from 5 m2 was measured in each plot and used for statistical analysis.
Analysis of chlorophyll content and photosynthetic rate in flag leaves
To determine chlorophyll content, leaf tissues were harvested using a circular punch that yields 0.5-cm diameter leaf discs. There were four flag leaf replicates for each treatment. Chlorophyll was extracted from wheat flag leaves obtained from the field at EZhou using 95% (v/v) ethanol (analytically pure, Sinopharm Chemical Reagent Co., Ltd) and the extracted chlorophyll concentration was determined using a spectrophotometer (UV2102, Unico, Shanghai, China) [25].
For the photosynthetic rate, flag leaf samples were obtained from the field at EZhou. Each treatment had three flag leaf replicates. Photosynthetic rate determination was performed as previously described [26].
Assay of plant hormones
Five frozen flag leaf and spike replicates from each treatment (~100 mg for each flag leaf and spike sample) were ground to a fine powder in liquid nitrogen using a mortar and pestle. Each sample was weighed into a 1.5-mL tube, mixed with 750 μL of cold extraction buffer (methanol: water: acetic acid, 80:19:1, v/v/v) supplemented with internal standards, 10 ng of 2H6ABA, 10 ng of DHJA, and 5 μg of NAA, vigorously shaken on a shaking bed for 16 h at 4 °C in the dark, and then centrifuged at 12,000 rpm for 15 min at 4 °C. Supernatant was carefully transferred to a new 1.5-mL tube and pellets remixed with 400 μL of extraction buffer, shaken for 4 h at 4 °C, and centrifuged. The two supernatants were combined and filtered using a syringe-facilitated 13-mm diameter nylon filter with a pore size of 0.22 μm (Nylon 66; Jinteng Experiment Equipment Co., Ltd, Tianjing, China). The filtrate was dried by evaporation under nitrogen gas flow for ~5 h at room temperature and then dissolved in 200 μL of methanol. Aliquots of dissolved samples were further diluted 40 times using methanol for jasmonic acid (JA), abscisic acid (ABA), and indole-3-acetic acid (IAA) quantification. Liquid chromatography was carried out using an ultrafast liquid chromatography with an autosampler (Shimadzu Corporation, Kyoto, Japan). The method used for hormone determination was as previously described [27].
Inoculation of F. graminearum and rating of disease
Infection assays on flowering wheat spikes were performed as previously described [28]. At the early flowering stage, a conidial suspension of F. graminearum strain PH-1 was collected from 5-day-old cultures growing in carboxymethylcellulose medium, then filtered through three layers of lens-wiping paper and then mixed with 0.01% (v/v) Tween 20. Ten microliters of 1 × 105 conidia mL−1 conidial suspension was inoculated individually onto the fourth spikelet from the bottom using a micropipette. The inoculated wheat spikes were maintained at a relative humidity of 95% for 72 h. Symptomatic spikes were examined and images captured after 14 days.
For the greenhouse test, 15 spikes from each treatment were inoculated and the spikelet infection rate of each spike was calculated; then, the average spikelet infection rate for each treatment was calculated for statistical analysis. For the field test, ten spikes from each plot were inoculated from a total of 30 inoculated spikes in the strain DT-8-treated plots and 30 spikes in the nontreated plots. The spikelet infection rate for each plot was calculated and the average spikelet infection rate for strain DT-8 treatment and nontreated control were then calculated for statistical analysis. The field test was conducted twice, once in 2017 and repeated in 2018.
FHB survey in a natural, noninoculated field
To investigate natural FHB infection in strain DT-8-treated wheat, an FHB field survey was conducted in experimental fields located at EZhou City, Jingmen City, and Xiangyang City in 2018. The field survey protocol described by the National Agricultural Technology Extension Service Center of China was adopted for the FHB survey with minor modifications. A total of 500 spikes were sampled randomly in each plot; in total, 1500 spikes were collected from DT-8-treated plots, with the same number of spikes being collected from nontreated plots to calculate disease incidence (spikelet infection rate) and severity (disease index). The number of infected and noninfected spikelets on each spike was counted and the average spikelet infection rate for each plot was calculated and used for statistical analysis. To calculate the disease index, the infected spikes were divided into five grades, namely: grade 0, no spikelet was infected; grade 1, the spike was infected, but less than 25% of spikelets were infected; grade 2, more than 25%, but less than 50%, of the spikelets were infected; grade 3, more than 50%, but less than 75%, of the spikelets were infected; and grade 4, more than 75% of the spikelets were infected. Finally, the disease severity for each plot was calculated using a formula for disease index, DI = ∑(nX/4 N) × 100, where “X” is the scale value of each spike, “n” is the number of spikes in the category, and “N” is the total number of spikes assessed for each plot. The disease index for each group was used for statistical analysis.
Inoculation of M. oryzae on barley and rice
To probe if S. sclerotiorum could enhance resistance against the rice blast fungus (M. oryzae) in barley and rice, an inoculation test was carried out according to a method described by Kong et al. [29]. Conidia of M. oryzae were collected with sterile water from 4-day-old cultures growing on tomato–oat medium and then filtered through three layers of lens-wiping paper. Infection assays were performed in whole plant leaves by spray inoculation using an airbrush nebulizer compressor. Strain DT-8-treated seedlings, which were grown in a sterile mixture of vermiculite and perlite at the ratio of 3:1 in pots and placed in greenhouse for 9 days for barley at 20 °C and 20 days for rice at 28 °C, were sprayed with a conidial suspension [105 conidia/mL mixed with 0.02% (v/v) Tween 20], using 4 mL suspension for each pot. The plants were further incubated at 28 °C, 80% relative humidity, under a 16 h light/8 h darkness photoperiod. Then, we assessed the presence of S. sclerotiorum in aerial parts of DT-8-treated barley and rice by PCR amplification in greenhouse plants; 92% of samples were confirmed as positive for S. sclerotiorum. For barley, lesions of leaves with the same leaf age (bottom leaves) were examined and typical infected leaves were photographed with a digital camera at 5 days post inoculation (dpi); for rice, the leaves were examined and photographed at 7 dpi.
Detection of toxins (DON)
To assay point-inoculated spikelets from strain DT-8-treated and nontreated plants, each sample was placed in a 50-mL centrifuge tube and mixed with 400 μL of a mixed isotope internal standard. The solution was remixed with 20 mL of an acetonitrile–water solution after standing for 30 min, vigorously shaken on a shaking bed for 4 h at 4 °C, and then centrifuged at 10,000 rpm for 5 min. The supernatants were carefully transferred to a new tube. The supernatants were combined and filtered using a syringe-facilitated 13-mm diameter nylon filter with a pore size of 0.22 μm (Nylon 66; Jinteng Experiment Equipment Co., Ltd, Tianjing, China). The filtrate was dried by evaporation under the nitrogen gas flow for ~5 h at room temperature, and then dissolved in 200 μL of methanol. The method for DON determination was performed as described by the National Agricultural Technology Extension Service Center of China (GB5009.111-2016).
RNA sequencing and analysis
Sterilized wheat seeds were inoculated with strain DT-8, and then were sown in experimental fields located at EZhou City in 2017. Wheat flag leaves and spikes were collected from this field during the initial bloom stage, with three spikes and leaves being randomly sampled from each of the three replicate plots. The flag leaves and spikes of nontreated wheat plants in the same field were randomly taken from each of the three control replicate plots. The samples were immediately placed in liquid nitrogen and ground into powder. Total RNA samples were extracted with a TRIzol Plus RNA Purification Kit (Takara, Dalian, China) and treated with RNase-free DNase I (Takara, Dalian, China) according to the manufacturer’s instructions. The RNA quality was checked using a Nanodrop Spectrophotometer (Thermo Fisher Scientific Inc., Wilmington, DE, USA). Then, mRNA was enriched with magnetic beads Oligo (dT) (TransGen Biotech, Beijing, China). Subsequently, cDNA was synthesized using the mRNA as template. The cDNA fragments were linked with adapters, and suitable fragments were selected for PCR amplification. Agilent 2100 Bioanalyzer and ABI StepOnePlus Real-Time PCR Systems were used in the quantification and qualification of the sample library. Subsequently, the library was sequenced for raw data using an Illumina HiSeq X sequencer at BGI (The Beijing Genomics Institute, China). Then, adapters, low-quality sequences, and reads with high content of unknown base (N) reads were removed to obtain clean reads. The clean reads were then mapped to the wheat genome or S. sclerotiorum genome and the sequence results evaluated in terms of read quality, alignment, saturation, and the distribution of reads on reference genes [30]. Mismatches of no more than two bases were accepted in the alignment. Gene expression was calculated by the number of reads mapped to the reference genomes using the fragments per kilobase of transcript per million mapped reads method [31]. Subsequently, differentially expressed genes (DEGs) were selected with FDR  More