Ecology

Global Health Program, Smithsonian Conservation Biology Institute, Washington DC, USAJames M. Hassell & Dawn ZimmermanDepartment of Epidemiology of Microbial Disease, Yale School of Public Health, New Haven, CT, USAJames M. Hassell & Dawn ZimmermanCentre for Biodiversity & Environment Research (CBER), Department of Genetics, Evolution and Environment, University College London, London, UKTim Newbold & Lydia H. V. FranklinosDepartment of Ecology & Evolutionary Biology, Princeton University, Princeton, NJ, USAAndrew P. DobsonSanta Fe Institute, Santa Fe, NM, USAAndrew P. DobsonWalter Reed Biosystematics Unit (WRBU), Smithsonian Institution Museum Support Center, Suitland, MD, USAYvonne-Marie LintonDepartment of Entomology, Smithsonian National Museum of Natural History, Washington DC, USAYvonne-Marie LintonWalter Reed Army Institute of Research (WRAIR), Silver Spring, MD, USAYvonne-Marie LintonMarine Disease Ecology Laboratory, Smithsonian Environmental Research Center, Edgewater, MD, USAKatrina M. Pagenkopp Lohan More

Plant growthPrimary root length and rosette size were estimated. Control root length, measured until 21 DAG, was lightly minor in ddc vs WT (Fig. 1A). Cd differentially inhibited root growth in the two samples: at 21 DAG, 25 and 50 µM Cd-treated roots were 1.2 and 2.2 fold shorter than control roots in ddc, while in the WT Cd-treated roots were 1.8 and 2.8 fold shorter than control ones (Fig. 1A). Consequently, at 21 DAG root of Cd-treated samples was longer in ddc vs WT, particularly at the lowest Cd concentration.Figure 1(A) Primary root length (B) Picture of rosette leaf series and (C) rosette leaf area (cm2) of WT and ddc plants of A. thaliana, germinated and grown for 21 DAG in long day condition: (i) on growth medium added with 25 or 50 µM Cd; (ii) on growth medium without Cd as control (Ctrl). Root length was monitored up to 21 days after germination (DAG) every two days from germination. The results represent the mean value (± SD) of three independent biological replicates (n = 45). Asterisks indicate significant pairwise differences using Student’s t-test (*P ≤ 0.05; ** P ≤ 0.01; *** P ≤ 0.001), performed between ddc vs WT subjected to the same treatment. Bars, 0.5 cm.Full size imageRosette size was estimated at 21 DAG, corresponding to the period necessary for its full development11, by evaluating leaf number and area. Control plants of both ddc and WT exhibited a complete leaf series, although most leaves resulted smaller in ddc (Fig. 1B,C). Cd affected rosette development reducing leaf number and area, less in ddc than WT, resulting into a higher leaf area and/or number in ddc under both Cd concentrations (Fig. 1B,C).Gene expression profileRNA-Seq analysis provided an overview of gene expression profile of Cd-treated and control plants of both ddc and WT. The following comparisons were performed: ddc vs WT under control (Ctrl) conditions (ddc vs WT-Ctrl) and 25 and 50 µM Cd treatment (ddc vs WT-25 µM Cd; ddc vs WT-50 µM Cd); 25/50 µM Cd-treated vs Ctrl in ddc (25 µM Cd vs Ctrl-ddc; 50 µM Cd vs Ctrl-ddc); 25/50 µM Cd-treated vs Ctrl in the WT (25 µM Cd vs Ctrl-WT; 50 µM Cd vs Ctrl-WT).After DEGs identification (see Supplementary Fig. S1 online) 14 of them were analysed through qRT-PCR to validate transcriptomic analysis (see Supplementary Fig. S2 online). Results were fully consistent with RNA-seq data. Gene Enrichment analysis was also performed, evidencing that Cd strongly impacted on transcriptome in both ddc and WT, but in a largely different way (see Supplementary Figs. S3–S9 online). Notwithstanding, a common aspect was that in both ddc and WT the genetic pathways (GPs) more impacted by Cd dealt with photosynthesis, stress responses and hormone biosynthesis and signalling.Expression pattern of genetic pathways related to hormonesIn view of hormones pivotal role in plant development and stress response and considering the assessed epigenetic control on their action and signalling12, in this work we analysed in depth how the expression pattern (EP) of hormone-related GPs was modulated in ddc vs WT under Cd stress. The most relevant differences are discussed.AuxinsUnder control conditions, GPs related to auxin biosynthesis showed comparable EP in ddc and WT and no DEGs were detected (Fig. 2A). 25 µM Cd induced significant changes only in ddc resulting into: i) TAA1 and YUC5 downregulation along indole-3-pyruvic acid (IPA) pathway; ii) CYP71A13 and NIT2 overexpression along indole-3-acetaldoxime (IAOX) auxiliary pathway, while CYP79B3 was downexpressed (Fig. 2A). Differently, 50 µM Cd induced similar changes in ddc and WT consisting in: (i) a downexpression of YUC2 along IPA pathway in both samples and YUC5 and YUC9 in ddc and WT, respectively; (ii) overexpression of CYP71A12, CYP71A13, NIT2, NIT4 and downexpression of CYP71A16 along IAOX pathway in ddc and WT (Fig. 2A).Figure 2Genes differentially expressed (DEGs) along the pathway of (A) auxin biosynthesis, auxin conjugation, (B) indole-3-acetyl-amino acid biosynthesis, (C) methyl-indole-3-acetate interconversion and (D) auxin signalling in ddc and WT plants identified through a transcriptomic approach. For each comparison, the log2(fold change) of the analysed DEGs was shown in orange and in blue for the upregulated and downregulated genes, respectively. Plants were grown for 21 DAG in long day condition: (i) on growth medium added with 25 or 50 µM Cd; (ii) on growth medium without Cd as control (Ctrl).Full size imageAuxin level and homeostasis also depend on its oxidative degradation, conjugation and methylation13. Under control conditions, GPs related to auxin conjugation and methylation showed comparable EPs in ddc and WT and no DEGs were detected (Fig. 2B,C), but were differentially impacted by Cd, mainly at 25 µM concentration. Namely, at 25 µM Cd several genes related to auxin conjugation (GH3.3, GH3.17, YDK1) and methylation (MES7, MES17) were downregulated in ddc, whereas in WT only MES18, involved in methyl-indole-3-acetate production, was downregulated (Fig. 2B,C). At 25 µM Cd most of the above genes were downregulated in ddc vs WT (Fig. 2B,C), while at 50 µM Cd only two genes, working in auxin methylation were differentially modulated in ddc and WT (Fig. 2B,C).Under control conditions, GP related to auxin signalling exhibited similar EP in ddc and WT (Fig. 2D). Differences were induced by Cd. Namely, at 25 µM Cd, AUX/IAA family genes, which acts in signalling repression, were globally downregulated in ddc, while in WT the repressor IAA34 was overexpressed (Fig. 2D). Differently, 50 µM Cd effects on ddc and WT were quite similar, dealing with AUX/IAA family genes downregulation (Fig. 2D). In ddc vs WT comparisons, at 25 µM Cd above genes were downregulated, while no differences were found at 50 µM Cd (Fig. 2D). Moreover, and somehow unexpectedly, following 50 µM Cd it was observed in both ddc and WT a downregulation of several SAURs members, belonging to a large family of auxin responsive genes, which in turn can also have an impact on auxin pathway14 (Fig. 2D). Interestingly, such effect was more pronounced in ddc than WT. However, it must be mentioned that, although most of them are induced by auxin, several other hormones and co-factors acts upstream SAUR genes, regulating their activity in response to both endogenous stimuli and environmental cues14.In summary, in both ddc and WT, Cd induced: i) a downregulation of IPA pathway, which is the main auxin biosynthetic pathway15 and a simultaneous upregulation of IAOX auxiliary biosynthetic pathway; ii) an enhancement of hormone signalling. However, in the WT such effects occurred only at 50 μM Cd. Moreover, in ddc Cd also induced a downregulation of GPs related to auxin conjugation.CytokininsIn all comparisons, GPs related to CKs biosynthesis showed similar EPs, unless for the downregulation in 50 µM Cd-treated WT vs Ctrl of IPT5, encoding rate-limiting enzyme along the pathway16 (Fig. 3A).Figure 3Genes differentially expressed (DEGs) along the pathway of (A) trans-zeatin biosynthesis, (B) CKs degradation, (C) CKs N7- and N9-glucoside biosynthesis, (D) CKs O-glycosylation and (E) CKs signalling in ddc and WT plants identified through a transcriptomic approach. For each comparison, the log2(fold change) of the analysed DEGs was shown in orange and in blue for the upregulated and downregulated genes, respectively. Plants were grown for 21 DAG in long day condition: (i) on growth medium added with 25 or 50 µM Cd; (ii) on growth medium without Cd as control (Ctrl).Full size imageMajor differences were observed for GPs related to CKs catabolism and conjugation, occurring through cleavage by oxidation and glycosylation, respectively16. Under control conditions, these GPs also exhibited similar EPs in ddc and WT (Fig. 3B-D). At both 25 and 50 µM Cd a downregulation of CKX5 and CKX6, encoding cytokinin-oxidases, occurred only in ddc (Fig. 3B). Differently, Cd impact on GPs related to CKs N-glycosylation was almost comparable in ddc and WT, resulting into the overexpression of two different genes working in N7- and N9-glycosylation pathways at 25 µM Cd, and one gene at the higher concentration (Fig. 3C).Note that at 25 µM Cd the above genes were both overexpressed in ddc vs WT, while no differences were found at 50 µM Cd (Fig. 3C). Cd impact on GP related to cytokinin O-glycosylation was major, especially in ddc, involving at 25 µM Cd the overexpression of seven genes along this pathway compared to Ctrl (Fig. 3D), and one gene in WT 25 µM Cd vs Ctrl (Fig. 3D). By contrast, at 50 μM Cd the expression pattern along this pathway was similar in ddc and WT, being characterised by the upregulation of the same seven genes above mentioned and the downregulation of AT5G38010 (Fig. 3D). Finally, at 25 μM Cd, five genes along these pathways resulted upregulated in ddc vs WT while a similar EP occurred in ddc and WT under 50 μM Cd treatment (Fig. 3D).Concerning GP involved in CKs signalling, under control conditions A-ARRs, encoding negative regulators of CKs signalling17, were downregulated and signalling was likely enhanced in ddc vs WT (Fig. 3E). 25 µM Cd induced a downregulation of ARR11 A-type ARRs and B-ARR family ARR10 transcription factors, which control primary plant response to CKs, only in ddc (Fig. 3E). Whereas, at 50 µM Cd both ddc and WT showed A-ARR downregulation, supposedly leading to pathway upregulation (Fig. 3E). At 25 µM Cd AHP1, encoding positive regulators of CKs signalling18, was downregulated in ddc vs WT, while ARR17 was overexpressed, suggesting that signalling was downregulated also in ddc vs WT (Fig. 3E). No differences occurred between ddc and WT at 50 µM Cd (Fig. 3E).In summary, transcriptomic analysis evidenced that GP related to the biosynthesis of trans-zeatin, the most relevant CK, was negatively affected by Cd only in the WT at 50 μM Cd. In response to Cd, GPs related to CKs inactivation were enhanced in both ddc and WT, but in ddc a downregulation of GP related to CKs cleavage also occurred. Finally, hormone signalling was differentially modulated by Cd in relation to both the sample (ddc vs WT) and heavy metal concentration, resulting into a downregulation at 25 µM Cd only in ddc and an enhancement in ddc and WT at 50 µM Cd.GibberellinsUnder control conditions, GPs related GAs biosynthesis showed similar EPs in ddc vs WT (Fig. 4A). 25 μM Cd induced in ddc: i) a downregulation of GA2 encoding the ent-kaurene synthase, a pivotal enzyme along the early GAs biosynthetic pathways to synthetize GA12; ii) a downregulation of GA4, a key gene of GAs biosynthesis, along which bioactive GAs are synthetized19 (Fig. 4A). No Cd-induced modulation was observed in the WT (Fig. 4A). On the contrary, at 50 μM Cd both ddc and WT showed a downregulation of GA5 (Fig. 4A). Finally, in ddc vs WT the only difference dealt with GA2 downregulation at 25 µM Cd (Fig. 4A).Figure 4Genes differentially expressed (DEGs) along the pathway of (A) GAs biosynthesis, (B) GAs inactivation and (C) GAs signalling in ddc and WT plants identified through a transcriptomic approach. For each comparison, the log2(fold change) of the analysed DEGs was shown in orange and in blue for the upregulated and downregulated genes, respectively. Plants were grown for 21 DAG in long day condition: (i) on growth medium added with 25 or 50 µM Cd; (ii) on growth medium without Cd as control (Ctrl).Full size imageGPs controlling GAs inactivation also showed a comparable transcriptional pattern in ddc and WT under control conditions, and no DEGs were detected (Fig. 4B). 25 μM Cd induced a downregulation of DAO2 and AOP1, encoding GA2ox enzymes, only in ddc (Fig. 4B). Accordingly, in ddc vs WT these genes were downregulated only at the lowest Cd concentration (Fig. 4B).Under control conditions, also GP related to GAs signalling was not differentially modulated in ddc vs WT (Fig. 4C). A downregulation of genes encoding DELLA proteins, which act as repressors20, was induced only in ddc by 25 µM Cd (Fig. 4C) and in both ddc and WT at 50 µM Cd (Fig. 4C), suggesting an enhancement of hormone signalling. Consequently, DELLA-codifying genes resulted downregulated also in ddc vs WT only at the lowest Cd concentration (Fig. 4C).In summary, in ddc the lowest Cd treatment negatively affected GPs related to GAs biosynthesis but, at the same time, hormone signalling resulted enhanced. In the WT similar effects were observed only at 50 µM Cd.Jasmonic acidUnder control conditions, six genes along the GP related to JA biosynthesis were downregulated in ddc vs WT (Fig. 5A). 25 µM Cd induced in ddc a downregulation of this GP, except for LOX4 upregulation (Fig. 5A) and a downregulation involving eight genes in WT (Fig. 5A). At 50 μM, Cd effects were limited to LOX5 downregulation and LOX4 and OPR1 upregulation in ddc and LOX5 and AOS downregulation in WT (Fig. 5A). No Cd–induced differences were found in ddc vs WT (Fig. 5A).Figure 5Genes differentially expressed (DEGs) along the pathway of (A) JA biosynthesis, (B) JA signalling, (C) ABA biosynthesis, (D) ABA degradation, (E) ABA glucose ester biosynthesis and (F) ABA signalling in ddc and WT plants identified through a transcriptomic approach. For each comparison, the log2(fold change) of the analysed DEGs was shown in orange and in blue for the upregulated and downregulated genes, respectively. Plants were grown for 21 DAG in long day condition: (i) on growth medium added with 25 or 50 µM Cd; (ii) on growth medium without Cd as control (Ctrl).Full size imageConcerning the JA signalling-related GP, in control conditions JAZ5 gene, encoding a protein acting as repressor21, was downregulated in ddc vs WT, highlighting a signalling enhancement (Fig. 5B). Interestingly, JAZ10 and JAZ9 were differentially impacted by 25 µM Cd in ddc and WT resulting upregulated and downregulated, respectively. (Fig. 5B). At 50 µM Cd, JAZs were overexpressed in both ddc and WT (Fig. 5B). No differences were detected in ddc vs WT exposed to Cd (Fig. 5B).Globally, Cd negatively impacted on the GP related to JA biosynthesis especially in WT. Under Cd treatment hormone signalling was downregulated more in ddc than in WT, whatever concentration was applied.Abscisic acidUnder control conditions, ABA biosynthesis-related GP showed comparable EP in ddc and WT, and no DEGs were detected (Fig. 5C). The only significant Cd effect dealt with NCED3 downregulation both in ddc and WT, regardless of applied concentration (Fig. 5C). Under control conditions, also the GPs related to ABA catabolism showed a comparable EP in ddc and WT and no DEGs were detected (Fig. 5D), but Cd differentially impacted on CYP genes, involved in phaseic acid degradative production22. Namely, at 25 μM Cd, CYP707A3 was downregulated only in ddc (Fig. 5D). Moreover, also CYP707A2 appeared downregulated in ddc vs WT (Fig. 5D). At 50 µM Cd it was observed an upregulation of both CYP707A2 and CYP707A4 in ddc, and of only CYP707A4 in WT (Fig. 5D).Under control conditions, GP related to ABA inactivation through glucose conjugation showed similar EP in ddc vs WT (Fig. 5E). 25 μM Cd determined AT4G15260 upregulation and UGT71C3 downregulation only in ddc (Fig. 5E). At 50 μM, Cd equally impacted on ddc and WT, resulting into UGT71C1 and UGT2 downregulation, AT5G49690, UGT71B5, UGT71B6 upregulation and, limited to WT, AT4G15260 upregulation (Fig. 5E). No differences were highlighted in ddc vs WT exposed to Cd (Fig. 5E).Under control conditions, the GP related to hormone signalling also presented a comparable EP in ddc and WT (Fig. 5F). In ddc, 25 µM Cd impact on this GP appeared rather complex, resulting in an upregulation of PYL3, encoding ABA receptor, and a downregulation of PP2Cs (PP2CA and HAI1) encoding negative regulators of ABA signalling23. Moreover, ABI5, codifying a key transcription factor in ABA signalling24 belonging to AREBs/ABFs family, was upregulated. However, SnRK2.7 gene, codifying a protein which activate the AREBs/ABFs transcription factors24, was downregulated. Based on the prominent role of SnRK2s in plant response to ABA, it is likely that at 25 µM Cd ABA signalling was downregulated in ddc (Fig. 5F). Instead, 25 µM Cd determined in WT the upregulation of PYL6 and the downregulation of PP2Cs, suggesting an enhancement of ABA signalling (Fig. 5F). At 50 µM Cd, PYL6 was upregulated and SnRK2.7 downregulated in both ddc and WT, while HAI1 was downregulated only in ddc (Fig. 5F). When comparing ddc vs WT, at 25 µM Cd only SnRK2.7 was downregulated, while no differences occurred at 50 μM Cd (Fig. 5F).In summary, Cd determined a slight downregulation of GP related to ABA biosynthesis in both ddc and WT regardless of its concentration. ABA catabolic pathway was lightly downregulated in ddc at 25 µM Cd but upregulated in both samples at 50 μM Cd. At the transcriptomic level, ABA signalling featured as enhanced in WT and downregulated in ddc regardless of Cd concentration.EthyleneAlong GP related to ethylene biosynthesis, under control conditions ACS8 and ACS11 were upregulated in ddc vs WT (Fig. 6A). 25 µM Cd determined ACS8 and ACO5 downregulation in ddc and ACS4 upregulation in the WT (Fig. 6A). 50 μM Cd induced ACS7 upregulation and ACO5 downregulation in both ddc and WT and ACS2 and ACS11 overexpression only in WT (Fig. 6A). Finally, the only Cd-induced difference in ddc vs WT dealt with ACO1 downregulation at 25 µM Cd (Fig. 6A).Figure 6Genes differentially expressed (DEGs) along the pathway of (A) ethylene biosynthesis, (B) ethylene signalling and (C) SA signalling in ddc and WT plants identified through a transcriptomic approach. For each comparison, the log2(fold change) of the analysed DEGs was shown in orange and in blue for the upregulated and downregulated genes, respectively. Plants were grown for 21 DAG in long day condition: (i) on growth medium added with 25 or 50 µM Cd; (ii) on growth medium without Cd as control (Ctrl).Full size imageGP related to ethylene signalling showed similar EP in ddc and WT both under control conditions (Fig. 6B) and at 25 µM Cd, except for the upregulation of ETR2, encoding ethylene receptor25 in WT (Fig. 6B). At 50 µM Cd, both ddc and WT exhibited ETR2 and ERF1 overexpression suggesting an upregulation of ethylene signalling (Fig. 6B); no differences occurred in ddc vs WT (Fig. 6B).In summary, in control conditions GP related to ethylene biosynthesis was upregulated in ddc vs WT. Cd determined a downregulation and upregulation of this GP in ddc and WT, respectively. Concerning hormone signalling, at the highest Cd concentration in both ddc and WT an upregulation of this GP occurred.Salicylic acidRegarding SA, only the GP related to signalling resulted differentially expressed. Under control conditions, the GP related to SA signalling showed similar EP in ddc and WT. However, PR1, a useful molecular marker for the systemic acquired resistance (SAR) in response to pathogens26, was downregulated in ddc vs WT (Fig. 6C). 25 μM Cd induced a downregulation of genes codifying TGA10 transcription factor only in the ddc (Fig. 6C). Whereas, 50 μM Cd induced a downregulation of PRB1 in both ddc and WT and of TGA8 only in ddc (Fig. 6C). In ddc vs WT, differences were found only at 25 µM Cd, with the downregulation of TGA10 and CAPE3 (Fig. 6C).Altogether, these results evidenced a Cd-induced downregulation of this GP, likely resulting in an impairment of hormone signalling in both WT and ddc, but in the latter this effect already occurred at the lowest Cd concentration.Phytohormone levelBased on the major effects induced by 25 µM Cd treatment, hormone quantification was carried out on plants exposed to this concentration, compared to untreated control plants.Under control conditions, IAA amount was higher in ddc than WT, although not significantly. After Cd treatment, a decreasing trend was observed only in WT, resulting into a significant lower level as compared to ddc (Fig. 7A).Figure 7(A) IAA, (B–G) CKs, (H–R) GAs, (S) JA, (T) ABA, (U) SA and (V) SAG amount in A. thaliana ddc mutant and WT plants grown in Ctrl conditions and treated with 25 µM Cd estimated by GC–MS. The results represent the mean value (± SD) of three independent biological replicates. Statistical analysis was performed by using two-way ANOVA with Tukey post hoc test (P ≤ 0.05) after Shapiro–Wilk normality test. Means with the same letter are not significantly different at P ≤ 0.05.Full size imageConcerning CKs, both biological active (tZ) and inactive conjugate (tZR, cZR, tZOG, cZOG, iPR) forms were analysed (Fig. 7B–G). Under control conditions, all analysed CKs were present in ddc, but CKs conjugate forms and above all O-glycosylated exhibited the highest levels (Fig. 7B–G). By contrast, in WT tZ was not detectable and all the other CKs forms exhibited a lower level compared to the mutant, which appeared significant for tZR and cZOG (Fig. 7B–G). Following Cd treatment, CKs levels increased in ddc, except for tZ decrease. In WT, the unique Cd effect dealt with tZR increase and tZOG decrease. Consequently, under Cd treatment the level of all CKs forms remained higher (from 0.25 to 3 times) in ddc than in WT (Fig. 7B–G).Concerning GAs, precursors (GA9, GA19, GA20), biologically active forms (GA1, GA3, GA4, GA7) and catabolites (GA8, GA34, GA29, GA51) were analysed (Fig. 7H–R). Under control conditions, GA19 amount was significantly higher in ddc vs WT, while the amount of GA20, the other in serie precursor of hydroxylated forms, was comparable between the samples (Fig. 7H,I). Following Cd treatment, GA19 amount significantly increased in WT while a slight downtrend occurred in ddc, leading to comparable values in the two samples. The same trend was observed, but at less extent, for GA20 (Fig. 7H,I). Consistently, under control conditions also the amount of the active hydroxylated forms GA1 and GA3 was higher in WT than in ddc (Fig. 7J,K). Following Cd treatment, a decrease of their amount was detected only in WT, globally leading to a higher level of these GAs in ddc mutant compared to the WT (Fig. 7J,K). In addition, in ddc mutant also the related catabolites GA8 and GA29 were globally lower than in WT, under both control conditions and Cd treatment (Fig. 7L,M).Differences were observed also for GA9, precursor of non-hydroxylated GAs: under Cd treatment its amount decreased in the WT and was instead induced in ddc mutant, resulting in a quite comparable value between the two samples (Fig. 7N). Consistently, the amount of active non-hydroxylated forms, GA4 and GA7, increased under Cd treatment only in ddc mutant; also in this case, at the end of heavy metal treatment, comparable values were detected in ddc and WT (Fig. 7O,P) In agreement with these results, following Cd treatment, the amount of catabolites GA51 and GA34 did not change in the WT, whereas in ddc it increased and decreased, respectively (Fig. 7Q,R).As evident in Fig. 7S–V, differences were reported also for JA, ABA, SA and its predominant inactive conjugate, SA 2-O-β-D-glucoside (SAG). Under control conditions both JA and ABA amount was significantly lower in ddc vs WT and significantly decreased following Cd treatment only in the WT (Fig. 7S,T). Notwithstanding, under such condition the ABA amount remained lower in ddc than in WT while JA values were comparable in the two samples due its light, but not significant, increase in ddc (Fig. 7S,T). By contrast, under control conditions both SA and SAG amounts were significantly higher in ddc than in WT (Fig. 7U,V). Following Cd treatment, their amounts significantly decreased more in ddc than in WT, leading to an opposite condition (Fig. 7U,V).Testing of the involvement of SUPPRESSOR OF DRM1 DRM2 CMT3 (SDC) gene in ddc response to CdFinally, we planned to inquire on the possible involvement of SDC gene in the response of ddc triple mutant to Cd exposure. Indeed, it has been reported that in ddc mutant the misexpression of such gene, which encodes a F-Box protein, is ultimately responsible of the developmental phenotypes of ddc, such as curled leaves and reduced growth, as evidenced by its reversion in the drm1 drm2 cmt3 sdc quadruple mutant27. Note that in the WT SDC is silenced, being methylated in all its sequence contexts because of the redundant action of DRM2 and CMT3 enzymes. By contrast, in ddc, where DRM2 and CMT3 expression is silenced, the loss of non-CG methylation in the promoter region of SDC F-box gene determines its overexpression27.According to the above mentioned data27, we firstly verified that under control conditions SDC resulted silent in the WT and overexpressed in ddc also in our transcriptomic analysis (data not shown, complete raw transcriptomic data are available at NCBI SRA under the BioProject accessionPRJNA641242;https://www.ncbi.nlm.nih.gov/Traces/study/?acc=PRJNA641242). Moreover, our data also showed that at the transcriptomic level SDC is not modulated by Cd since its expression level did not significantly change in ddc nor in WT whatever heavy metal concentration was applied.Thereafter, we tested the involvement of SDC in the growth response of ddc mutant under Cd exposure, by monitoring primary root length of Arabidopsis thaliana WT, ddc and sdc plants grown under the following conditions: (i) on growth medium without Cd as control (Ctrl) (ii) on a medium supplemented with 25/50 µM Cd; (iii) limited to the WT and sdc mutant plants, on a medium supplemented with 25/50 µM Cd plus 15 µM 5-Azacytidine (5-Aza), an inhibitor of DNA methylation applied in order to mimic the hypomethylated state of ddc mutant.Under control conditions, all three samples showed a similar root length. However, at 21 DAG, root was lightly shorter in ddc vs WT, while sdc displayed an intermediate length (Fig. 8A). Under both Cd treatments, roots were averagely longer in ddc than in WT. Again, sdc roots exhibited an intermediate length, more similar to WT than ddc (Fig. 8 B,C). Interestingly, WT and sdc plants treated with Cd plus 5-Aza had longer roots than the plants treated only with Cd, and quite comparable to ddc roots exposed to Cd (Fig. 8 D,E).Figure 8Primary root length of WT, sdc and ddc plants of A. thaliana, germinated and grown in long day condition (A) on growth medium without Cd as control (Ctrl), (B,C) on a medium supplemented with 25/50 µM Cd, (D,E) limited to WT and sdc plants, on a medium supplemented with 25/50 µM Cd plus 15 µM 5-Azacytidine (5-Aza). Root length was monitored up to 21 days after germination (DAG) every two days from germination. The results represent the mean value (± SD) of three independent biological replicates (n = 45). Statistical analysis was performed between samples at the same growth stage, by using two-way ANOVA with Tukey post hoc test (P ≤ 0.05) after Shapiro–Wilk normality test. Means with the same letter are not significantly different at P ≤ 0.05.Full size image More

Genotyping control and datasets creationsAfter filtering the initial raw dataset of 1,071 goats for poor quality genotype and related animals we obtained the dataset used for the haplotype sharing analysis, which included 980 animals and 48,396 SNPs. This dataset was further pruned for linkage disequilibrium (LD; r2  More

MDE hypothesisGradient features of species diversity of plant communities refer to regular changes in species diversity along a gradient of environmental factors at the community level12,43. The altitudinal gradient includes gradient effects of multiple environmental factors. It is therefore important to study altitudinal patterns of species diversity to reveal changes in biodiversity along environmental gradients. The width and range of species distribution along geographical gradients reflect the ecological adaptability, diffusivity, and evolutionary history of species44. To some extent, geographical distribution patterns of species diversity can be interpreted as outcomes of synthetic actions across altitudinal gradients resulting from eurychoric species with a greater distribution width and stenochoric species with a smaller distribution range along geographical gradients44. Hence, the MDE, environmental gradient, distribution area, human disturbance, and habitat heterogeneity all have effects on the vertical distribution patterns of species diversity45,46. According to the MDE hypothesis, there is overlap in the distribution range of species along altitudinal gradients; the highest overlap intensity occurs at intermediate elevations. There is relatively low overlap intensity at low and high elevation areas, and the peak values of species diversity occur at intermediate elevations.In this study, forest ecosystems on the Loess Plateau were separated into three communities: tree, shrub, and herb communities. The altitudinal patterns and factors that influence species diversity of mountainous vegetation were then determined at the plant community level related to the form and family of the plant species. We discovered that the family numbers of the herb, shrub, and tree communities reached their greatest values at intermediate, intermediate, and lower elevations, respectively. We also discovered that correlations of species diversity indices with elevations conformed to unimodal change patterns for herb, shrub, and tree communities, which presented their greatest values at higher, lower, and intermediate elevations, respectively. This showed that MDE is another important factor that affects the distribution patterns of species diversity along regional altitudinal gradients, in addition to temperature, precipitation, and the terrain.A large number of studies have already verified that MDE is a significant mechanism that influences the gradient patterns of species diversity. MDE not only functions along altitudinal gradients but also acts along latitudinal and temporal gradients28,29,30,47. However, the effects of MDE on species diversity patterns are highly controversial. Some studies considered MDE to be the main factor that results in unimodal patterns of species diversity47,48, whereas other studies affirmed that the effects of MDE are smaller in contrast to the functions of the distribution area, environmental gradient, and other factors29,30. Besides MDE, other factors may also lead to unimodal vertical gradient patterns of species diversity, such as plant productivity, human disturbance, and the regional climate45,49. Our research indicated that the relationships of species diversity conformed to unimodal change patterns along various elevations for mountainous herb, shrub, and tree communities in a semi-arid region of the temperate zone. It can be concluded that vertical patterns of species diversity with a unimodal type may be a more universal phenomenon, relative to monotonic decreasing or increasing patterns of species diversity with increasing elevations.Factors influencing MDE at the plant community levelIn this study, more forbs and grasses were found at higher elevations, whereas more sedges occurred at lower elevations. The responses of the importance values of tree species to the altitudinal gradient demonstrated the following variation patterns: evergreen coniferous trees had higher importance values than deciduous coniferous trees, followed by deciduous broad-leaved trees. This showed that MDE was influenced by species life form. The species diversity of different life forms responded differently to the environment, and plant species with different life forms presented different diversity patterns along altitudinal gradients50. In New Zealand, the number of mountainous plants species decreased with increasing elevations and the total species number of all plants also decreased significantly, while species diversity had no significant distribution trend in response to elevation when plants with different life forms were considered under different layers of plant communities51.MDE is a common pattern of species diversity of mountainous plants with changing elevations. Our study area, located on the Loess Plateau of China, belongs to a semi-arid mountainous region in the temperate zone where the maximum species diversity of the tree community occurred at intermediate elevations. This finding was in accordance with the MDE hypothesis. Research from the Kinabalu Mountains in Sabah, Malaysia, indicated an obvious MDE pattern of species number linked to elevations52. On the Haleakala Mountains, Hawaii, USA, the highest species diversity occurred at intermediate elevations, which was also in accordance with the MDE hypothesis53. The MDE hypothesis was also proved by studies conducted in the Yu Mountains of Taiwan and the Emei Mountains of Szechwan in China54. However, the MDE pattern of species diversity in tree communities is caused by precipitation, which is the highest at intermediate elevations52. This situation also occurred in the herb community.There are many factors that affect the distribution of herbaceous plants, so the variations in species diversity in the herb community with elevations are complex55. In this study, we found that the herb community exhibited higher species diversity at higher elevations; more forbs and grasses were distributed at higher elevations, whereas more sedges were distributed at lower elevations. In the Siskiyou Mountains in Oregon, USA, the species diversity of herbaceous plants had a significantly positive correlation with elevation. This correlation occurred mainly due to an increase in the number of grass species, which was the primary reason that radiation was enhanced by a drastic reduction in community coverage as a result of increasing elevations. Consequently, there was an increase in the species diversity of herbaceous plants24. A decrease in species diversity with increasing elevations is a more familiar pattern for herbaceous plants in temperate56 and tropical22 forests.We also discovered that the family numbers of herb and shrub species all showed unimodal change patterns with high values at their central elevations in this semi-arid region. This conformed to the MDE hypothesis as well. In arid temperate grasslands, species diversity indicated an MDE distribution pattern. For example, the species diversity of herbaceous plants presented an MDE pattern in drought areas of the Siskiyou Mountains57. However, in semi-humid mountains in the temperate zone, the species diversity of herbaceous plants was principally in control of the community structure, and community coverage did not respond uniformly to elevation. Studies in New Zealand showed that there were no evident distribution trends for species diversity of herbaceous plants along elevations51. In low bush communities of Chile, the species diversity of herbaceous plants declined with increasing elevations after longstanding succession, but it increased during the early stage of succession58. Therefore, relationships between the species diversity of herb species and elevations were not completely clear in semi-humid regions.The major factors that control the distribution areas of species differ among different families and genera, and thus the vertical distribution patterns of species diversity differ22. We concluded that the family number of the tree species had a maximum distribution at lower elevations, unlike herb and shrub species; meanwhile, the responses of the importance values to the altitudinal gradient in the tree community were also different among evergreen coniferous trees, deciduous coniferous trees, and deciduous broad-leaved trees. These differences may have been related to environmental factors. Due to various distribution patterns of environmental parameters with elevations, the distribution patterns of species diversity showed large changes along elevations59. For example, the distribution of fern and Melastomataceae species is principally related to humidity, that of Acanthaceae and Bromeliaceae species is correlated with temperature, and that of Araceae species is related to transpiration59. Research conducted in the Gongga Mountains, China, showed that the species diversity with different floral components exhibited different distribution patterns along elevations due to differences in the environment and species origin60. We also discovered that the importance values of dominant families in the shrub (Rosaceae) and tree (Pinaceae) communities exhibited changing patterns in contrast to MDE. In our study, only the family numbers in the herb and shrub communities, as well as species diversity in the tree community, conformed to the MDE hypothesis. Therefore, we concluded that the MDE hypothesis of species diversity of mountainous vegetation is influenced by the species life form and family of different plant communities in the temperate semi-arid region of China.Factors influencing plant species diversity at the environment levelThe altitudinal distribution patterns of the plant community diversity had greater discrepancies in mountainous regions and between different community types, which might be connected to regional environmental conditions, relative heights of mountains, and the geological landscape35. Concerning the altitudes of mountains, serious human disturbances (e.g., deforestation, grazing, and land-use conversion) had negative effects on biodiversity in low-altitude regions61. In high-altitude regions, the cold climate slowed down plant growth and soil development, while other harsh environments exceeded the tolerance limits for growth of the majority of species, such as by intense solar radiation or large temperature differences between day and night62. In the middle-altitude regions, the species diversity was relatively higher due to less human disturbances and the formation of transition zones of plant species differentiation between the low- and high-altitude regions62. Hence, the plant community diversity and its altitudinal gradient patterns in mountainous regions were largely influenced by regional climate and human disturbances.Comparisons of the diversity at different levels indicated that the responses of the plant community diversity to the environment were not the same for diverse gradations, and different species exhibited different gradient patterns owing to restrictions from environmental factors63,64. The primary factors leading to the altitudinal differentiation of diversity included the temperature, moisture, soil nutrients, and succession process65. In our study area, compared with Guancen Mountain and Guandi Mountain, Wulu Mountain at the lower latitude of the Lvliang Mountains had a lower altitude and was located in the continental monsoon subhumid climate region of the warm temperate zone, making it suitable for the growth of secondary forests and shrub vegetation. However, the vegetation growth in the herb layer was restricted in Wulu Mountain, making diversity in the herb layer the greatest on Guandi Mountain at the middle latitude of the Lvliang Mountains35. This showed that the plant species diversity in the east of the Loess Plateau changed with the altitude, while being affected by complicated habitat conditions such as latitude and human disturbances. This characterization of the study area was correlated with the unimodal patterns observed. In this study, we observed that the family numbers of the herb and shrub communities presented unimodal patterns across an altitudinal gradient; the importance values of dominant families also presented unimodal patterns in the shrub and tree communities; and the species diversity indices of the herb, shrub, and tree communities conformed to unimodal change patterns following an altitudinal gradient as well.In our recent studies on the species diversity of herbaceous communities in the Lvliang Mountains66, we found that the results calculated for β-diversity using different indices revealed the highest value for the Cody index and the lowest value for the Bray–Curtis index at altitudes between 1900 and 2000 m, indicating that areas located between 1900 and 2000 m form a transition zone in which the herbaceous community undergoes a rapid process of species renewal and changes in its composition. The results for γ-diversity indicated a pattern of unimodal variation in relation to altitude. Changes in altitude gradient had highly significant impacts on changes in temperature and humidity, indicating that various environmental factors (notably humidity and temperature) and human disturbances had combined effects on changes in the values of the α-diversity indices.At present, it is widely believed that the formation of herbs in different life forms was principally impacted by precipitation, whereas in areas with similar rainfall, water, heat, and light conditions need to be considered. These conditions chiefly included average annual precipitation, accumulated temperature, and illumination time67. In this research, we observed that herb and tree species in different life forms showed different trends with altitudinal gradients in the Lvliang Mountains. At higher elevations, forbs and grasses grew well, whereas sedges grew well at lower elevations. The responses of different tree life forms to the altitudinal gradient were greater for evergreen coniferous tree species than for deciduous coniferous tree species and deciduous broad-leaved tree species. From the north to the south in the Lvliang Mountains, increases in the average annual precipitation increased the number of species and components of the annual herbs, while the hydrothermal matching requirements of Guandi Mountain at the middle latitude were preferred for annual herb growth in comparison with Guancen Mountain at the higher latitude67. However, considering whole mountains, the Lvliang Mountains located in the continental monsoon climate region of the warm temperate zone had four distinctive seasons with drought and wind in the spring, and a quick rise in air temperatures, and larger diurnal temperature differences35. These conditions conformed to the habitat features of herbs and trees. Hence, the hydrothermal distribution status affected by latitude and human disturbances determined the altitudinal distribution patterns of plant community diversity in the Lvliang Mountains.The MDE at different elevationsIn this study, we discussed the MDE of mountainous vegetation on the Loess Plateau with an elevational range from 1324 to 2745 m, including tree, shrub and herb community. This range was a very large elevation range for a case study, but a very short range in comparison to global elevation ranges, which extended from the sea level to well over 8000 m (though the highest locations did not have any vegetation). Therefore, owing to this limitation, the results we obtained in this research were suitable for lower elevation mountains in semi-arid areas.The MDE was changed with different elevations and vegetation layers. In studies on the Daiyun Mountains with an elevation from 900 to 1600 m, the phylogenetic diversity and species diversity of tree community indicated an intermediate high expansion pattern along elevations and their peak values all appeared at the elevation of 1200 m68. This conclusion conformed to the MDE pattern. In our studies on the Lvliang Mountains with an elevation from 1459 to 2610 m, higher species diversity of tree community was observed at intermediate elevations with a peak value being at the 2000 m, which conformed to the MDE pattern either. Therefore, at smaller elevations less than 2600 m, species diversity of tree conformed to the MDE pattern.When an elevation reached 2700 m on the Lvliang Mountains, the vegetation types changed to shrub and grass, and only their family numbers followed the MDE pattern across an altitudinal gradient. Slimily, the species richness of shrub and herb layer showed an obvious “lateral pattern” on an elevational gradient from 2950 to 3750 m on the Three River Headwater, both reaching the maximum value at the 3150 m; while with the rise of altitude, α diversity of shrub layer and herb layer showed a “wave”-shaped changed trend, reaching the lowest value at the 3550 m69. It illustrated that species diversity of shrub and grass did not conform to the MDE pattern completely at medium elevations from 2700 to 3700 m.As for an elevation extending from 3000 to 4400 m on the Gongga Mountains, the vegetation type was alpine meadow, and the species richness index presented an obvious unimodal pattern with a peak value at the 3850 m, which accorded with the MDE pattern70. Similarly, studies on an alpine meadow on the Gannan revealed that the number of richness, Shannon-Weiner index and phylogenetic diversity of plant community all showed a “humped-back” relationship with the increase of altitude from 3000 to 4000 m, and reached the maximum value at the 3800 m71,72. Thereby, at greater elevations more than 3800 m, species diversity of alpine meadow conformed to the MDE pattern.However, when an elevation exceeded 5000 m, research object on species diversity were not vegetation but animals along an altitudinal gradient. For example, on the Himalaya Mountains with an elevation from 3755 to 5016 m, the ant species richness illustrated a “unimodal curve” pattern along the rise of altitude, and the Shannon–Wiener index and Fisher α index of ant community commonly expressed the “Multi-Domain Effect” phenomenon73. Another research on mammalian richness was also conducted on the Himalaya Mountains. It concluded that most of elevational species richness patterns were hump-shaped from 100 to 6000 m on the Himalayas Mountains74. As a result, the MDE pattern was also extremely common in animal communities. More

The catenary sequence of soilsThe catena of Andosols down a slope near a hot spring in the Valley of Geysers was subdivided into four thermal zones (Fig. 1a–e), which are described below.Figure 1Location of study area (a–c), soil pits along the catena (d, e) and photos of soil pits (f–i). (a b) Study area location. Soil map is from133 (open access) with additions and corrections by I.N. Semenkov based on the map of soil temperature at a depth of 15 cm in the Valley of Geysers131 and the soil names from44 using CorelDraw X7 software (https://www.coreldraw.com/). (c) Top view of the left side of the Geysernaya River with the body of a catastrophic landslide, a visitor center (in the left lower part) and the location of the transect studied (1–15). (d) The location of transect studied. (e) A schematic profile of the catena studied with numbered soil pits. The main soil pits selected for comprehensive analyses (see section ‘Soil analyses’) are in red. (f) Non-heated Eutrosilic Silandic Andosols (Arenic, Cutanic) on pyroclastic material (pit no 16, Zone I), within levelled parts of the interfluve, under tall-herb meadow communities with local patches of Erman’s birch woods. (g) Slightly heated Eutrosilic Aluandic Andosols (Cutanic, Loamic, Natric) in the upper part of the catena, on hydrothermally altered sandy-loamy pyroclastic material (pit no 12, Zone II), on slightly heated slopes, under tall-herb meadows. (h) Moderately heated Eutrosilic Gleyic Aluandic Andosols (Loamic, Reductic, Protosalic, Hyperthionic) in the middle part of the catena, on hydrothermally altered clayey pyroclastic material (pit no 9.1, Zone III), under different moss and ‘microzonal’ communities. (i) Hot Gleyic Aluandic Andosols (Clayic, Reductic, Salic, Hyperthionic) in the lower part of the catena, on hydrothermal clays (pit no 4, Zone IV), on most heated bare slopes.Full size imageZone I. Non-heated Eutrosilic Silandic Andosols (Arenic, Cutanic) under Kamchatka’s tall herb communities and fragmented Erman’s birch woodsNon-heated Andosols with temperatures of 50 with very small amounts of PM5–50 ( More

Soil properties in different rice farming systemsFive treatments were designed in the three selected rice fields, including (1) rice monoculture field (RM); (2) planting area in the 1st year of rice–fish field (OP); (3) aquaculture area in the 1st year of rice–fish field (OA); (4) planting area in the 5th year of rice–fish field (FP); (5) aquaculture area in the 5th year of rice–fish field (FA). The soil properties of the five treatments were shown in Table 1. The highest soil available nitrogen (AN) content was observed in FP and was significantly higher than that in RM, OP and OA. The highest soil available phosphorus (AP) content was observed in RM and was significantly higher than that in the other 4 treatments. The highest soil available potassium (AK) content was measured in the 1st year of rice–fish field (OP and OA), followed by RM and the 5th year of rice–fish field (FP and FA), and significant differences were observed among different rice fields. The highest soil organic matter (OM) content appeared in the 5th year of rice–fish field (FP and FA), and was only significantly higher than that in OA. In addition, the soil pH in the 1st year of rice–fish field (OP and OA) was significantly lower than that in RM and the 5th year of rice–fish field (FP and FA). In summary, significant differences of soil properties were observed among the different rice farming systems.Table 1 Soil properties in different rice systems and areas.Full size tableSoil bacterial community compositionA total of 1,346,468 sequences were obtained by 16S rRNA MiSeq sequencing analysis after basal quality control (reads containing ambiguous bases were discarded; only overlapping sequences longer than 10 bp were assembled; Operational taxonomic units (OTUs) were clustered with 97% similarity). These sequences were classified as 46 phyla, 800 genera and 5335 OTUs. As shown in Fig. 1, the dominant bacterial phyla across different treatments were Proteobacteria (26.06–29.41%) and Chloroflexi (20.07–27.99%), followed by Actinobacteria (7.22–20.87%), Acidobacteria (11.36–14.46%) and Nitrospirae (3.11–8.50%). Since the implementation of rice–fish farming regime, the soil bacterial community composition has greatly changed. For example, Actinobacteria abundance decreased from 20.87% in RM to 7.22% in FA, while Nitrospirae abundance greatly increased from 3.11% in RM to 8.50% in FA. Between different areas in a same rice–fish field (i.e. OP vs OA or FP vs FA), the bacterial community composition were similar. The PCoA analysis on OTU level also showed that different areas within the same rice–fish field had high similarity in bacterial community composition. In contrast, the bacterial community composition differed distinctly among different rice farming systems (Fig. 2). Bacterial alpha diversity indices, as evaluated by Shannon, Simpson, ACE and Chao1, were shown in Table 2. Student’s t-test was adopted to evaluate the difference among treatments. The results showed that the alpha indices of FP were significantly lower than other treatments, except for Simpson index.Figure 1The average relative abundances on phylum level of soil bacterial communities in different rice systems and areas.Full size imageFigure 2PCoA analysis on OTU level based on bray_curtis distance algorithm (significance among treatments were conducted with ANOSIM test, R = 0.4294, P = 0.0010).Full size imageTable 2 Alpha diversity indices of soil bacterial in different rice systems and areas.Full size tableBased on the Kruskal–Wallis test, the statistical differences among treatments were evaluated in the abundances of the top 15 phyla. The results showed that 5 phyla, including Actinobacteria, Nitrospirae, Bacteroidetes, Unclassified_k_norank and SBR1093 were observed significant differences among treatments, and the most significant phylum was Nitrospirae (Fig. 3). In order to trace the source of the significant differences, the Wilcoxon tests were conducted between every two rice cultivation patterns separately (Fig. 4). The results indicated that the significant differences were mainly derived from the comparison between RM and F_group (FP & FA), as well as the comparison between the O_group (OP & OA) and F_group. In the comparison between the RM and O_group, only the phylum Gemmatimonadetes was observed to have a significant difference. Furthermore, we also compared the differences of the top 15 phyla between planting area (P_group) and aquaculture area (A_group) within rice–fish fields, and the results showed no phyla observed with significant differences in the abundances.Figure 3The differences with significance of the top 15 phyla in different rice systems and areas (* indicates 0.01  More

1.Cózar, A. et al. Plastic debris in the open ocean. Proc. Nat. Acad. Sci. USA 111, 10239–10244. https://doi.org/10.1073/pnas.1314705111 (2014).ADS 
CAS 
Article 
PubMed 

Google Scholar 
2.Lavers, J. L., Dicks, L., Dicks, M. R. & Finger, A. Significant plastic accumulation on the Cocos (Keeling) Islands, Australia. Sci. Rep. 9, 7102. https://doi.org/10.1038/s41598-019-43375-4 (2019).ADS 
CAS 
Article 
PubMed 
PubMed Central 

Google Scholar 
3.Cózar, A. et al. The arctic ocean as a dead end for floating plastics in the north atlantic branch of the thermohaline circulation. Sci. Adv. https://doi.org/10.1126/sciadv.1600582 (2017).Article 
PubMed 
PubMed Central 

Google Scholar 
4.Peeken, I. et al. Arctic sea ice is an important temporal sink and means of transport for microplastic. Nat. Commun. 9, 1505. https://doi.org/10.1038/s41467-018-03825-5 (2018).ADS 
CAS 
Article 
PubMed 
PubMed Central 

Google Scholar 
5.Woodall, L. C. et al. The deep sea is a major sink for microplastic debris. R. Soc. Open Sci. 1, 140317 (2014).ADS 
Article 

Google Scholar 
6.Chiba, S. et al. Human footprint in the abyss: 30 year records of deep-sea plastic debris. Mar. Policy 96, 204–212. https://doi.org/10.1016/j.marpol.2018.03.022 (2018).Article 

Google Scholar 
7.Bergmann, M., Tekman, M. & Gutow, L. Sea change for plastic pollution. Nature 544, 297 (2017).ADS 
CAS 
Article 

Google Scholar 
8.Jambeck, J. R. et al. Plastic waste inputs from land into the ocean. Science 347, 768–771. https://doi.org/10.1126/science.1260352 (2015).ADS 
CAS 
Article 
PubMed 

Google Scholar 
9.Gall, S. C. & Thompson, R. C. The impact of debris on marine life. Mar. Pollut. Bull. 92, 170–179. https://doi.org/10.1016/j.marpolbul.2014.12.041 (2015).CAS 
Article 
PubMed 

Google Scholar 
10.Camphuysen, C. J. Northern Gannets Morus bassanus found dead in the Netherlands, 1970–2000. Atlantic Seabirds 3, 15–30 (2001).
Google Scholar 
11.Gregory, M. R. Environmental implications of plastic debris in marine settings–entanglement, ingestion, smothering, hangers-on, hitch-hiking and alien invasions. Phil. Trans. R. Soc. B 364, 2013–2025 (2009).Article 

Google Scholar 
12.Ryan, P. G. The effects of ingested plastic on seabirds: Correlations between plastic load and body condition. Environ. Pollut. 46, 119–125 (1987).CAS 
Article 

Google Scholar 
13.Ryan, P. G. Effects of ingested plastic on seabird feeding: Evidence from chickens. Mar. Pollut. Bull. 19, 125–128 (1988).Article 

Google Scholar 
14.Pierce, K. E., Harris, R. J., Larned, L. S. & Pokras, M. A. Obstruction and starvation associated with plastic ingestion in a Northern Gannet Morus bassanus and a greater shearwater Puffinus gravis. Mar. Ornithol. 32, 187–189 (2004).
Google Scholar 
15.Ryan, P. G., Connell, A. D. & Gardner, B. D. Plastic ingestion and PCBs in seabirds: Is there a relationship?. Mar. Pollut. Bull. 19, 174–176 (1988).CAS 
Article 

Google Scholar 
16.Lavers, J. L., Bond, A. L. & Hutton, I. Plastic ingestion by Flesh-footed Shearwaters (Puffinus carneipes): Implications for chick body condition and the accumulation of plastic-derived chemicals. Environ. Pollut. 187, 124–129. https://doi.org/10.1016/j.envpol.2013.12.020 (2014).CAS 
Article 
PubMed 

Google Scholar 
17.Tanaka, K. et al. In vivo accumulation of plastic-derived chemicals into seabird tissues. Curr. Biol. 30, 723-728.e3. https://doi.org/10.1016/j.cub.2019.12.037 (2020).CAS 
Article 
PubMed 

Google Scholar 
18.Teuten, E. L. et al. Transport and release of chemicals from plastics to the environment and to wildlife. Phil. Trans. R. Soc. B 364, 2027–2045 (2009).CAS 
Article 

Google Scholar 
19.Tanaka, K., van Franeker, J. A., Deguchi, T. & Takada, H. Piece-by-piece analysis of additives and manufacturing byproducts in plastics ingested by seabirds: Implication for risk of exposure to seabirds. Mar. Pollut. Bull. 145, 36–41. https://doi.org/10.1016/j.marpolbul.2019.05.028 (2019).CAS 
Article 
PubMed 

Google Scholar 
20.Thiel, M. & Gutow, L. The ecology of rafting in the marine environment. I. The floating substrata. Oceanogr. Mar. Biol. Annu. Rev. 42, 181–264 (2005).
Google Scholar 
21.Kiessling, T., Gutow, L. & Thiel, M. Marine litter as habitat and dispersal vector. In: Bergmann M, Gutow L, Klages M, editors. Marine Anthropogenic Litter. p. 141–80 (2015).22.Day, R. H. & Shaw, D. G. Patterns of abundance of pelagic plastic and tar in the North Pacific Ocean, 1976–1985. Mar. Pollut. Bull. 18, 311–316 (1987).CAS 
Article 

Google Scholar 
23.Pichel, W. G. et al. Marine debris collects within the North Pacific Subtropical Convergence Zone. Mar. Pollut. Bull. 54, 1207–1211 (2007).CAS 
Article 

Google Scholar 
24.Yamashita, R. & Tanimura, A. Floating plastic in the Kuroshio Current area, western North Pacific Ocean. Mar. Pollut. Bull. 54, 485–488 (2007).CAS 
Article 

Google Scholar 
25.Titmus, A. J. & Hyrenbach, K. D. Habitat associations of floating debris and marine birds in the North East Pacific Ocean at coarse and meso spatial scales. Mar. Pollut. Bull. 62, 2496–2506 (2011).CAS 
Article 

Google Scholar 
26.Goldstein, M. C., Titmus, A. J. & Ford, M. Scales of spatial heterogeneity of plastic marine debris in the northeast pacific ocean. PLoS ONE 8, e80020 (2013).ADS 
Article 

Google Scholar 
27.Eriksen, M. et al. Plastic pollution in the world’s oceans: More than 5 trillion plastic pieces weighing over 250,000 tons afloat at sea. PLoS ONE 9, e111913 (2014).ADS 
Article 

Google Scholar 
28.IUCN. The IUCN Red List of Threatened Species. Version 2020–2. https://www.iucnredlist.org (2020).29.Lavers, J. L. & Bond, A. L. Ingested plastic as a route for trace metals in Laysan Albatross (Phoebastria immutabilis) and Bonin Petrel (Pterodroma hypoleuca) from Midway Atoll. Mar. Pollut. Bull. 110, 493–500. https://doi.org/10.1016/j.marpolbul.2016.06.001 (2016).CAS 
Article 
PubMed 

Google Scholar 
30.Roman, L., Hardesty, B. D., Hindell, M. A. & Wilcox, C. A quantitative analysis linking seabird mortality and marine debris ingestion. Sci. Rep. 9, 3202. https://doi.org/10.1038/s41598-018-36585-9 (2019).ADS 
CAS 
Article 
PubMed 
PubMed Central 

Google Scholar 
31.Jouventin, P. & Weimerskirch, H. Satellite tracking of wandering albatrosses. Nature 343, 746–748 (1990).ADS 
Article 

Google Scholar 
32.Kappes, M. A. et al. Hawaiian albatrosses track interannual variability of marine habitats in the North Pacific. Prog. Oceanogr. 86, 246–260 (2010).ADS 
Article 

Google Scholar 
33.Sakamoto, K. Q., Takahashi, A., Iwata, T. & Trathan, P. N. From the eye of the albatrosses: A bird-borne camera shows an association between albatrosses and a killer whale in the Southern Ocean. PLoS ONE 4, e7322 (2009).ADS 
Article 

Google Scholar 
34.Fukuoka, T. et al. The feeding habit of sea turtles influences their reaction to artificial marine debris. Sci. Rep. 6, 28015. https://doi.org/10.1038/srep28015 (2016).ADS 
CAS 
Article 
PubMed 
PubMed Central 

Google Scholar 
35.Nishizawa, B. et al. Albatross-borne loggers show feeding on deep-sea squids: Implications for the study of squid distributions. Mar. Ecol. Prog. Ser. 592, 257–265 (2018).ADS 
Article 

Google Scholar 
36.Hunt, G. L. Jr. & Schneider, D. Scale-dependent processes in the physical and biological environment of marine birds. In Seabirds: Feeding Ecology and Role in Marine Ecosystems (ed. Croxall, J. P.) 7–41 (Cambridge University Press, 1987).
Google Scholar 
37.Pinaud, D. & Weimerskirch, H. At-sea distribution and scale-dependent foraging behaviour of petrels and albatrosses: A comparative study. J. Anim. Ecol. 76, 9–19 (2007).Article 

Google Scholar 
38.Thiebot, J.-B., Nishizawa, B., Sato, F., Tomita, N. & Watanuki, Y. Albatross chicks reveal interactions of adults with artisanal longline fisheries within a short range. J. Ornithol. 159, 935–944 (2018).Article 

Google Scholar 
39.Froese, R. & Pauly, D. FishBase. World Wide Web electronic publication. www.fishbase.org, version (12/2019).40.Ryan, P. G. A simple technique for counting marine debris at sea reveals steep litter gradients between the Straits of Malacca and the Bay of Bengal. Mar. Pollut. Bull. 69, 128–136 (2013).CAS 
Article 

Google Scholar 
41.Mitani, Y. et al. Marine debris observed in the North Pacific during Oshoro-maru cruise in 2012. Bull. Fish. Sci. Hokkaido Univ. 64, 25–29 (2014).
Google Scholar 
42.Hyrenbach, K. D. et al. Plastic ingestion by Black-footed albatross from Kure Atoll, Hawai’i: linking chick loads and parental at-sea distributions. Mar. Ornithol. 45, 225–236 (2017).
Google Scholar 
43.Nevitt, G. A., Losekoot, M. & Weimerskirch, H. Evidence for olfactory search in wandering albatross, Diomedea Exulans. Proc. Nat. Acad. Sci. USA 105, 4576–4581 (2008).ADS 
CAS 
Article 

Google Scholar 
44.Savoca, M. S., Wohlfeil, M. E., Ebeler, S. E. & Nevitt, G. A. Marine plastic debris emits a keystone infochemical for olfactory foraging seabirds. Sci. Adv. 2, e1600395 (2016).ADS 
Article 

Google Scholar 
45.Santos, R. G., Andrades, R., Fardim, L. M. & Martins, A. S. Marine debris ingestion and Thayer’s law—The importance of plastic color. Environ. Pollut. 214, 585–588 (2016).CAS 
Article 

Google Scholar 
46.Castro, J. J., Santiago, J. A. & Santana-Ortega, A. T. A general theory on fish aggregation to floating objects: An alternative to the meeting point hypothesis. Rev. Fish Biol. Fish. 11, 255–277 (2002).Article 

Google Scholar 
47.Harrison, C. S., Hida, T. S. & Seki, M. P. Hawaiian seabird feeding ecology. Wildl. Monogr. 85, 1–71 (1983).
Google Scholar 
48.Hunte, W., Oxenford, H. A. & Mahon, R. Distribution and relative abundance of flyingfish (Exocoetidae) in the eastern Caribbean. II. Spawning substrata, eggs and larvae. Mar. Ecol. Prog. Ser. 117, 25–37 (1995).ADS 
Article 

Google Scholar 
49.Rapp, D. C., Youngren, S. M., Hartzell, P. & Hyrenbach, K. D. Community-wide patterns of plastic ingestion in seabirds breeding at French Frigate Shoals Northwestern Hawaiian Islands. Mar. Pollut. Bull. 123, 269–278 (2017).CAS 
Article 

Google Scholar 
50.Douglas, D. & Peucker, T. Algorithms for the reduction of the number of points required to represent a digitized line or its caricature. Cannadian Cartogr. 10, 112–122 (1973).Article 

Google Scholar 
51.Edelhoff, H., Signer, J. & Balkenhol, N. Path segmentation for beginners: an overview of current methods for detecting changes in animal movement patterns. Move. Ecol. 4, 21 (2016).Article 

Google Scholar 
52.R Core Team. R: A language and environment for statistical computing. R Foundation for Statistical Computing, Vienna, Austria. http://www.r-project.org/index.html (2020). More

1.Díaz, R. J. Overview of hypoxia around the world. J. Environ. Qual. 30, 275–281 (2001).PubMed 
Article 
PubMed Central 

Google Scholar 
2.Laffoley, D. & Baxter, J. M. (eds) Ocean deoxygenation: Everyone’s problem. Causes, impacts, consequences and solutions (IUCN, International Union for Conservation of Nature, 2019).
Google Scholar 
3.Booth, J. A. T. et al. Patterns and potential drivers of declining oxygen content along the southern California coast. Limnol. Oceanogr. 59, 1127–1138 (2014).ADS 
CAS 
Article 

Google Scholar 
4.Gilbert, D., Rabalais, N. N., Díaz, R. J. & Zhang, J. Evidence for greater oxygen decline rates in the coastal ocean than in the open ocean. Biogeosciences 7, 2283–2296 (2010).ADS 
CAS 
Article 

Google Scholar 
5.Altieri, A. H. & Gedan, K. B. Climate change and dead zones. Glob. Change Biol. 21, 1395–1406 (2015).ADS 
Article 

Google Scholar 
6.Breitburg, D. et al. Declining oxygen in the global ocean and coastal waters. Science (80-) 359, eaam7240 (2018).Article 
CAS 

Google Scholar 
7.Keeling, R. E., Körtzinger, A. & Gruber, N. Ocean deoxygenation in a warming world. Ann. Rev. Mar. Sci. 2, 199–229 (2010).PubMed 
Article 
PubMed Central 

Google Scholar 
8.Levin, L. A. & Breitburg, D. L. Linking coasts and seas to address ocean deoxygenation. Nat. Clim. Change 5, 401–403 (2015).ADS 
Article 

Google Scholar 
9.Rabalais, N. N., Turner, R. E., Díaz, R. J. & Justić, D. Global change and eutrophication of coastal waters. ICES J. Mar. Sci. 66, 1528–1537 (2009).Article 

Google Scholar 
10.Diaz, R. J. & Rosenberg, R. Spreading dead zones and consequences for marine ecosystems. Science (80-) 321, 926–929 (2008).ADS 
CAS 
Article 

Google Scholar 
11.Hofmann, A. F., Peltzer, E. T., Walz, P. M. & Brewer, P. G. Hypoxia by degrees: Establishing definitions for a changing ocean. Deep Res. Part I Oceanogr. Res. Pap. 58, 1212–1226 (2011).ADS 
CAS 
Article 

Google Scholar 
12.Rabalais, N. N. et al. Dynamics and distribution of natural and human-caused hypoxia. Biogeosciences 7, 585–619 (2010).ADS 
CAS 
Article 

Google Scholar 
13.Vaquer-Sunyer, R. & Duarte, C. M. Thresholds of hypoxia for marine biodiversity. Proc. Natl. Acad. Sci. 105, 15452–15457 (2008).ADS 
CAS 
PubMed 
Article 
PubMed Central 

Google Scholar 
14.Altieri, A. H. et al. Tropical dead zones and mass mortalities on coral reefs. Proc. Natl. Acad. Sci. U. S. A. 114, 3660–3665 (2017).CAS 
PubMed 
PubMed Central 
Article 

Google Scholar 
15.Grantham, B. A. et al. Upwelling-driven nearshore hypoxia signals ecosystem and oceanographic changes in the northeast Pacific. Nature 429, 749–754 (2004).ADS 
CAS 
PubMed 
Article 
PubMed Central 

Google Scholar 
16.Kim, T. W., Barry, J. P. & Micheli, F. The effects of intermittent exposure to low-pH and low-oxygen conditions on survival and growth of juvenile red abalone. Biogeosciences 10, 7255–7262 (2013).ADS 
Article 

Google Scholar 
17.Kolesar, S. E., Breitburg, D. L., Purcell, J. E. & Decker, M. B. Effects of hypoxia on Mnemiopsis leidyi, ichthyoplankton and copepods: Clearance rates and vertical habitat overlap. Mar. Ecol. Prog. Ser. 411, 173–188 (2010).ADS 
CAS 
Article 

Google Scholar 
18.Low, N. H. N. & Micheli, F. Lethal and functional thresholds of hypoxia in two key benthic grazers. Mar. Ecol. Prog. Ser. 594, 165–173 (2018).ADS 
CAS 
Article 

Google Scholar 
19.Thomas, P. & Saydur Rahman, M. Extensive reproductive disruption, ovarian masculinization and aromatase suppression in Atlantic croaker in the northern Gulf of Mexico hypoxic zone. Proc. R. Soc. B Biol. Sci. 279, 28–38 (2011).Article 
CAS 

Google Scholar 
20.Breitburg, D. Effects of hypoxia, and the balance between hypoxia and enrichment, on coastal fishes and fisheries. Estuaries 25, 767–781 (2002).Article 

Google Scholar 
21.Brown, J. H., Gillooly, J. F., Allen, A. P., Savage, V. M. & West, G. B. Toward a metabolic theory of ecology. Ecology 85, 1771–1789 (2004).Article 

Google Scholar 
22.Pörtner, H. O. & Knust, R. Climate change affects marine fishes through the oxygen limitation of thermal tolerance. Science (80-) 315, 95–97 (2007).ADS 
Article 
CAS 

Google Scholar 
23.Vaquer-Sunyer, R. & Duarte, C. M. Temperature effects on oxygen thresholds for hypoxia in marine benthic organisms. Glob. Change Biol. 17, 1788–1797 (2011).ADS 
Article 

Google Scholar 
24.Breitburg, D. L., Hondorp, D. W., Davias, L. A. & Diaz, R. J. Hypoxia, nitrogen, and fisheries: Integrating effects across local and global landscapes. Ann. Rev. Mar. Sci. 1, 329–349 (2009).PubMed 
Article 
PubMed Central 

Google Scholar 
25.Booth, J. A. T. et al. Natural intrusions of hypoxic, low pH water into nearshore marine environments on the California coast. Cont. Shelf. Res. 45, 108–115 (2012).ADS 
Article 

Google Scholar 
26.Walter, R. K., Woodson, C. B., Leary, P. R. & Monismith, S. G. Connecting wind-driven upwelling and offshore stratification to nearshore internal bores and oxygen variability. J. Geophys. Res. Ocean 119, 3517–3534 (2014).ADS 
Article 

Google Scholar 
27.Boch, C. A. et al. Local oceanographic variability influences the performance of juvenile abalone under climate change. Sci. Rep. 8, 1–12 (2018).CAS 
Article 

Google Scholar 
28.DiMarco, S. F., Chapman, P., Walker, N. & Hetland, R. D. Does local topography control hypoxia on the eastern Texas–Louisiana shelf?. J. Mar. Syst. 80, 25–35 (2010).Article 

Google Scholar 
29.Leary, P. R. et al. “Internal tide pools” prolong kelp forest hypoxic events. Limnol. Oceanogr. 62, 2864–2878 (2017).ADS 
CAS 
Article 

Google Scholar 
30.Walter, R. K., Brock Woodson, C., Arthur, R. S., Fringer, O. B. & Monismith, S. G. Nearshore internal bores and turbulent mixing in southern Monterey Bay. J. Geophys. Res. Ocean 117, 1–13 (2012).
Google Scholar 
31.Long, W. C. & Seitz, R. D. Trophic interactions under stress: Hypoxia enhances foraging in an estuarine food web. Mar. Ecol. Prog. Ser. 362, 59–68 (2008).ADS 
Article 

Google Scholar 
32.Kwiatkowski, L. & Orr, J. C. Diverging seasonal extremes for ocean acidification during the twenty-first centuryr. Nat. Clim. Chang. 8, 141–145 (2018).ADS 
CAS 
Article 

Google Scholar 
33.Safaie, A. et al. High frequency temperature variability reduces the risk of coral bleaching. Nat. Commun. 9, 1–12 (2018).Article 
CAS 

Google Scholar 
34.Woodson, C. B. The fate and impact of internal waves in nearshore ecosystems. Ann. Rev. Mar. Sci. 10, 421–441 (2018).CAS 
PubMed 
Article 

Google Scholar 
35.Woodson, C. B. et al. Harnessing marine microclimates for climate change adaptation and marine conservation. Conserv. Lett. 12(2), 1–9 (2018).
Google Scholar 
36.Micheli, F. et al. Evidence that marine reserves enhance resilience to climatic impacts. PLoS ONE 7, e40832 (2012).
ADS 
CAS 
PubMed 
PubMed Central 
Article 

Google Scholar 
37.Cox, K. W. California abalones, family haliotidae. Fish. Bull. 118 28–32 (1962).

Google Scholar 
38.Frieder, C. A., Nam, S. H., Martz, T. R. & Levin, L. A. High temporal and spatial variability of dissolved oxygen and pH in a nearshore California kelp forest. Biogeosciences 9, 3917–3930 (2012).ADS 
CAS 
Article 

Google Scholar 
39.Mayol, E., Ruiz-Halpern, S., Duarte, C. M., Castilla, J. C. & Pelegrí, J. L. Coupled CO2 and O2-driven compromises to marine life in summer along the Chilean sector of the Humboldt Current System. Biogeosciences 9, 1183–1194 (2012).ADS 
CAS 
Article 

Google Scholar 
40.Orellana-Cepeda, E., Granados-Machuca, C. & Serrano-Esquer, J. Ceratium furca: One possible cause of mass mortality of cultured Blue-Fin Tuna at Baja California, Mexico. Harmful Algae 2002, 514–516 (2004).
Google Scholar 
41.Bograd, S. J. et al. Oxygen declines and the shoaling of the hypoxic boundary in the California Current. Geophys. Res. Lett. 35, 1–6 (2008).Article 
CAS 

Google Scholar 
42.Bernardi, G., Findley, L. & Rocha-Olivares, A. Vicariance and dispersal across Baja California in disjunct marine fish populations. Evolution (N Y) 57, 1599–1609 (2003).
Google Scholar 
43.Haupt, A. J., Micheli, F. & Palumbi, S. R. Dispersal at a snail’s pace: Historical processes affect contemporary genetic structure in the exploited wavy top snail (Megastraea undosa). J. Hered. 104, 327–340 (2013).PubMed 
Article 
PubMed Central 

Google Scholar 
44.Al Najjar, M. W. Nearshore Processes of a Coastal Island: Physical Dynamics and Ecological Implications (Stanford University, 2019).
Google Scholar 
45.Hughes, B. B. et al. Climate mediates hypoxic stress on fish diversity and nursery function at the land-sea interface. Proc. Natl. Acad. Sci. U. S. A. 112, 8025–8030 (2015).ADS 
CAS 
PubMed 
PubMed Central 
Article 

Google Scholar 
46.Sydeman, W. J. et al. Climate change and wind intensification in coastal upwelling ecosystems. Science (80-) 345, 77–80 (2014).ADS 
CAS 
Article 

Google Scholar 
47.Fulton, S. et al. From fishing fish to fishing data: The role of Artisanal Fishers in Conservation and Resource Management in Mexico. In Viability and Sustainability of Small-Scale Fisheries in Latin America and The Caribbean (eds Salas, S. et al.) 151–175 (Springer International Publishing, 2019).
Google Scholar 
48.Chang, W., Cheng, J., Allaire, J. J., Xie, Y. & McPherson, J. shiny: Web Application Framework for R. R package version 1.4.0.2. https://cran.r-project.org/package=shiny (2020).49.R Core Team. R: A language and environment for statistical computing. R Foundation for Statistical Computing, Vienna, Austria. https://www.R-project.org/ (2020).50.Eerkes-Medrano, D., Menge, B. A., Sislak, C. & Langdon, C. J. Contrasting effects of hypoxic conditions on survivorship of planktonic larvae of rocky intertidal invertebrates. Mar. Ecol. Prog. Ser. 478, 139–151 (2013).ADS 
Article 

Google Scholar 
51.Low, N. H. N. & Micheli, F. Short- and long-term impacts of variable hypoxia exposures on kelp forest sea urchins. Sci. Rep. 10, 1–9 (2020).CAS 
Article 

Google Scholar 
52.Bograd, S. J. et al. Phenology of coastal upwelling in the California Current. Geophys. Res. Lett. 36, 1–5 (2009).Article 

Google Scholar 
53.Nam, S., Kim, H. J. & Send, U. Amplification of hypoxic and acidic events by la Nia conditions on the continental shelf off California. Geophys. Res. Lett. 38, 1–5 (2011).Article 
CAS 

Google Scholar 
54.Rogers-Bennett, L. et al. Dinoflagellate bloom coincides with marine invertebrate mortalities in Northern California. Harmful Algae News 46, 10–11 (2012).
Google Scholar 
55.Chan, F. et al. Persistent spatial structuring of coastal ocean acidification in the California Current System. Sci. Rep. 7, 1–8 (2017).Article 
CAS 

Google Scholar 
56.Montgomery, D. W., Simpson, S. D., Engelhard, G. H., Birchenough, S. N. R. & Wilson, R. W. Rising CO2 enhances hypoxia tolerance in a marine fish. Sci. Rep. 9, 1–10 (2019).CAS 
Article 

Google Scholar 
57.Boch, C. A. et al. Effects of current and future coastal upwelling conditions on the fertilization success of the red abalone (Haliotis rufescens). ICES J. Mar. Sci. 74, 1125–1134 (2017).Article 

Google Scholar 
58.Gobler, C. J. & Baumann, H. Hypoxia and acidification in marine ecosystems: Coupled dynamics and effects on
ocean life. Biol. Lett. 12, 20150976 (2016).PubMed 
PubMed Central 
Article 
CAS 

Google Scholar  More