Article

Coastal phytoplankton blooms expand and

intensify in the 21st century

https://doi.org/10.1038/s41586-023-05760-y Yanhui Dai1,9, Shangbo Yang1,9, Dan Zhao1, Chuanmin Hu2, Wang Xu3, Donald M. Anderson4,
Yun Li5, Xiao-Peng Song6, Daniel G. Boyce7, Luke Gibson1, Chunmiao Zheng1,8 & Lian Feng1 ✉
Received: 17 May 2022

Accepted: 25 January 2023

Phytoplankton blooms in coastal oceans can be beneficial to coastal fisheries
Published online: 1 March 2023
production and ecosystem function, but can also cause major environmental
Open access
problems1,2—yet detailed characterizations of bloom incidence and distribution are
Check for updates not available worldwide. Here we map daily marine coastal algal blooms between
2003 and 2020 using global satellite observations at 1-km spatial resolution. We found
that algal blooms occurred in 126 out of the 153 coastal countries examined. Globally,
the spatial extent (+13.2%) and frequency (+59.2%) of blooms increased significantly
(P < 0.05) over the study period, whereas blooms weakened in tropical and
subtropical areas of the Northern Hemisphere. We documented the relationship
between the bloom trends and ocean circulation, and identified the stimulatory
effects of recent increases in sea surface temperature. Our compilation of daily
mapped coastal phytoplankton blooms provides the basis for global assessments
of bloom risks and benefits, and for the formulation or evaluation of management
or policy actions.

Phytoplankton blooms are accumulations of microscopic algae in the ocean surface continuously since 1997 and have enabled bloom detec-
surface layer of fresh and marine water bodies. Although many blooms tion in many coastal regions17–19. However, the technical difficulties in
can occur naturally, nutrients linked to anthropogenic eutrophication dealing with complex optical features across different types of coastal
are expected to intensify their frequency globally2–4. Many algal blooms waters have so far prohibited their application globally20. To fill this
are beneficial, fixing carbon at the base of the food chain and support- knowledge gap, we developed a method to map global coastal algal
ing fisheries and ecosystems worldwide. However, proliferations of blooms and used this tool to examine satellite images between 2003
algae that cause harm (termed harmful algal blooms (HABs)) have and 2020, addressing three fundamental questions: (1) where and how
become a major environmental problem worldwide5–7. For instance, frequently global coastal oceans have been affected by phytoplankton
the toxins produced by some algal species can accumulate in the food blooms; (2) whether the blooms have expanded or intensified over the
web, causing closures of fisheries as well as illness or mortality of marine past two decades, both globally and regionally; and (3) the identity of
species and humans8–10. In other cases, the decay of a dense algal bloom the potential drivers.
can deplete oxygen in bottom waters, forming anoxic ‘dead zones’ that
can cause fish and invertebrate die-offs and ecosystem restructuring,
with serious consequences for the well-being of coastal communities1,11. Mapping global coastal phytoplankton blooms
Unfortunately, algal bloom frequency and distribution are projected We generated a satellite-based dataset of phytoplankton bloom occur-
to increase with future climate change12,13, with some changes causing rence to characterize the spatial and temporal patterns of algal blooms
adverse effects on aquatic ecosystems, fisheries and coastal resources. in coastal oceans globally. The dataset was derived using global, 1-km
Owing to substantial heterogeneity in space and time, algal blooms resolution daily observations from the Moderate Resolution Imag-
are challenging to characterize on a large scale5,14, and thus present ing Spectroradiometer (MODIS) onboard NASA’s Aqua satellite, and all
knowledge does not allow us to answer one of the most fundamen- 0.76 million images acquired by this satellite mission between 2003 and
tal questions: whether algal blooms have changed in recent decades 2020 were used. We developed an automated method to detect phyto-
on a global basis6,15,16. For example, although HAB events have been plankton blooms using MODIS images (Extended Data Fig. 1) (Methods).
compiled into the UNESCO (United Nations Educational, Scientific, In this study, we define a phytoplankton bloom as the accumulation of
and Cultural Organization) Intergovernmental Oceanographic Com- microscopic algae at the ocean surface that exhibits satellite-detectable
mission Harmful Algae Event Database (HAEDAT) globally since 1985, fluorescence signals21. However, whether a bloom produces toxins or is
bloom trends are difficult to resolve, owing to inconsistent sampling harmful to humans or the marine environment is not distinguishable
efforts and the diversity of the eco-environmental or socio-economic from satellite data. We delineated bloom-affected areas (that is, the areas
effects6. Alternatively, satellite data have been used to monitor the where algal blooms were detected), and enumerated the bloom count
1
School of Environmental Science and Engineering, Southern University of Science and Technology, Shenzhen, China. 2College of Marine Science, University of South Florida, St. Petersburg,
FL, USA. 3Shenzhen Ecological and Environmental Monitoring Center of Guangdong Province, Shenzhen, China. 4Woods Hole Oceanographic Institution, Woods Hole, MA, USA. 5School of
Marine Science and Policy, College of Earth, Ocean, and Environment, University of Delaware, Lewes, DE, USA. 6Department of Geographical Sciences, University of Maryland, College Park,
MD, USA. 7Bedford Institute of Oceanography, Fisheries and Oceans Canada, Dartmouth, Nova Scotia, Canada. 8EIT Institute for Advanced Study, Ningbo, China. 9These authors contributed
equally: Yanhui Dai, Shangbo Yang. ✉e-mail: fengl@sustech.edu.cn

280 | Nature | Vol 615 | 9 March 2023

 a

50

Bloom counts
30

10
1
b c Global total affected area:
30.3% 31.47 × 106 km2
10
60

Affected area (×106 km2)

21.5% 18.2%
Bloom counts

40
5
11.8%
9.0%
9.2%
20

0 0
SA AF EU NA AS AU Global EU NA AS SA AF AU

Fig. 1 | Global patterns of coastal phytoplankton blooms between 2003 and represents the median value, bottom and top bounds of boxes are first and
2020. a, The spatial distribution of annual mean bloom count based on daily third quartiles, and the whiskers show a maximum of 1.5 times the interquartile
satellite detections. b, Continental and global statistics for annual mean range. c, Continental statistics for the long-term annual mean of bloom-
bloom count (South America (SA), n = 3,846,125; Africa (AF), n = 2,516,225; affected areas (n = 18 years). The percentages show the corresponding
Europe (EU), n = 17,703,949; North America (NA), n = 10,034,286; Asia (AS), contributions to the global total. The bars represent s.d. Open circles are
n = 5,371,158; Australia (AU), n = 2,781,998 pixel observations). The centre line the affected areas during different years. Map created using Python 3.8.

at the 1-km pixel level (that is, the number of detected blooms per pixel) nutrient enrichment9,22–26. European LMEs showed mostly large propor-
(Fig. 1). We further estimated the bloom frequency (dimensionless) by tions of bloom-affected areas, and some also showed frequent bloom
integrating the bloom count and affected areas within 1° × 1° grid cells occurrences. By contrast, Asian LMEs exhibited mainly infrequent
(see Methods), and this metric was used to examine temporal dynam- blooms, given their large affected areas. We identified more bloom
ics in bloom intensity. Validation with independent satellite samples events in estuarine regions than along coasts in regions without major
selected via several visual inspection techniques showed an overall accu- river discharge (P < 0.05; Extended Data Fig. 8), highlighting the critical
racy level of more than 95% for our method, and comparisons using dis- role of terrestrial nutrient sources in coastal algal blooms3.
crete events in HAEDAT6 indicated that we successfully identified bloom
counts for 79.3% of the historical HAB events in that database (Extended
Data Figs. 2–6). We examined phytoplankton blooms in the exclusive Long-term trends
economic zones (EEZs) of 153 coastal countries and in 54 large marine The total global bloom-affected area has expanded by 3.97 million km2
ecosystems (LMEs) (Extended Data Fig. 7). Our study area encompasses (13.2%) between 2003 and 2020, equivalent to 0.14 million km2 yr−1
global continental shelves and outer margins of coastal currents, which (P < 0.05; Fig. 2). Furthermore, the number of countries with significant
offer the majority of marine resources available for human use (see Meth- bloom expansion was about 1.6 times those with a decreasing trend.
ods). Out of the 153 coastal countries examined, 126 were observed to The global median bloom frequency showed an increasing rate of 59.2%
have phytoplankton blooms (Fig. 1). The total bloom-affected area was (+2.19% yr−1, P < 0.05) over the observed period. Spatially, areas showing
31.47 million km2, equivalent to approximately 24.2% of the global land significant increasing trends (P < 0.05) in bloom frequency were 77.6%
area and 8.6% of the global ocean area, with a median bloom count of larger than those with the opposite trends (Fig. 2). Globally, our analysis
4.3 per year during the past 2 decades (Fig. 1b). Europe (9.52 million revealed overall consistent fluctuations between the bloom-affected
km2—30.3% of the total affected area) and North America (6.78 million area and bloom frequency between 2003 and 2020 (Fig. 2b). However,
km2—21.5% of the total affected area) contributed the largest bloom there was no significant relationship between bloom extent and fre-
areas. By contrast, the most frequent blooms were found around Africa quency in 23 countries and 10 LMEs over the past two decades, under-
and South America (median bloom counts of more than 6.3 per year). scoring the spatial and temporal variability of algal blooms and the
Australia experienced the lowest frequency (2.4 per year) and affected importance of continuous satellite monitoring.
area (2.84 million km2—9.0% of the total affected area) of blooms. The entire Southern Hemisphere was primarily characterized by
Phytoplankton blooms occurred frequently in the eastern boundary increased bloom frequency, although weakened blooms were also
current systems (that is, California, Benguela, Humboldt and Canary), sometimes found. In the Northern Hemisphere, the low latitude
northeastern USA, Latin America, the Baltic Sea, Northern Black Sea and (<30° N) coasts were mainly featured with strong bloom weakening
the Arabian Sea (Fig. 1a). Five LMEs were found with the most frequent (Fig. 2a), primarily in the California Current System and the Arabian
blooms (annual median bloom count over 15), including Patagonian Sea. Bloom strengthening was found in the northern Gulf of Mexico and
Shelf, Northeast US continental Shelf, the Baltic Sea, Gulf of California the East and South China Seas, albeit at smaller magnitudes. At higher
and Benguela Current (Extended Data Fig. 7). These hotspots are often latitudes, weakening blooms were detected mainly in the northeastern
reported as having a high incidence of algal blooms, some of which are North Atlantic and the Okhotsk Sea in the northwestern North Pacific.
HABs, driven by either coastal upwelling or pronounced anthropogenic Globally, the largest increases in bloom frequency were observed in six

Nature | Vol 615 | 9 March 2023 | 281

Article
a

50° N
Positive
Latitude (P < 0.05)
0° (P ≥ 0.05)
Negative
(P < 0.05)
50° S (P ≥ 0.05) Slope (104 yr −1)
100 0 100
Fraction (%) –1.0 –0.5 0 0.5 1.0

Affected area (×106 km2)

Bloom frequency (×103)

3 +0.14 × 106 km2 yr −1 (P = 0.001) 35

2 30
+2.19% yr −1 (P = 0.006)

1 25
2003 2006 2009 2012 2015 2018
Year

Fig. 2 | Trends of global coastal phytoplankton blooms between 2003 and b, Interannual variability and trends in annual median bloom frequency and total
2020. a, Spatial patterns of the trends in bloom frequency at a 1° × 1° grid scale. global bloom-affected area. The linear slopes and P-value (two-sided t-test) are
The latitudinal profiles show the fractions of grids with significant and indicated. The shading associated with the bloom frequency data represents
insignificant trends (positive or negative) along the east–west direction. an uncertainty level of 5% in bloom detection. Map created using Python 3.8.

major coastal current systems, including Oyashio (+6.31% yr−1), Alaska averaged over the growth window of algal blooms within a year (Meth-
(+5.22% yr−1), Canary (+4.28% yr−1), Malvinas (+3.02% yr−1), Gulf Stream ods and Extended Data Fig. 9)) in several high-latitude regions (>40° N),
(+2.42% yr−1) and Benguela (+2.30% yr−1) (Figs. 2a and 3). such as the Alaska Current (r = 0.44), the Oyashio Current (r = 0.48) and
the Baltic Sea (r = 0.41) (Fig. 3). These findings agree with previous stud-
ies, in which the bloom-favourable seasons in these temperate seas
Natural and anthropogenic effects have been extended under warmer temperatures27–29. However, this
Increases in sea surface temperature (SST) can stimulate bloom occur- temperature-based mechanism did not apply to regions with inconsistent
rence. We found significant positive correlations (P < 0.05) between the trends between SST and bloom frequency, particularly for the substantial
annual mean bloom frequency and the coincident SST (SST data were bloom weakening in the tropical and subtropical areas (Figs. 2a and 3b).

a b
8
7
2 6
1
3

5
4 10−6 °C m−1 decade−1 °C decade−1
–2 –1 0 1 2 –1.0 –0.5 0 0.5 1.0

c 1 California Current 2 Gulf Stream 3 Canary Current 4 Malvinas Current

12 10 28 8 24 16 6 10 24 20 11 15
S = –4.26 r = 0.61* r = –0.02 S = 4.28 r = 0.12

6 8 5 19 14 3 22 14 9 13
26 S = 2.42 8
Bloom frequency (104)

SST (10–6 ഒ per m)

S = 3.02
r = 0.81* r = –0.83* r = 0.84* r = –0.63* r = 0.83*
0 6 2 14 12 0 20 8 7 11
SST (ഒ)

2003 2011 2019 2003 2011 2019 2003 2011 2019 2003 2011 2019

5 Benguela Current 6 Oyashio Current 7 Alaska Current 8 Baltic Sea

40 16 21 10 16 14 4 11 9 11
r = 0.30
S = 5.22
8
11
25 r = 0.73* 13 5 14 12 2 9 6 7
r = 0.16
19
r = 0.58* 6
S = 2.30 S = 6.31 r = 0.48* r = 0.09 r = 0.44* S = 2.74 r = 0.41*
10 10 0 12 10 0 7 3 3 9
2003 2011 2019 2003 2011 2019 2003 2011 2019 2003 2011 2019
Year

Fig. 3 | Effects of climate change on phytoplankton blooms. a,b, Global and the correlation coefficient (r) between bloom frequency and the SST and
patterns of trends in SST gradient (a) and SST (b) from 2003 to 2020. c, Long- the SST gradient (∇SST) are shown. Asterisks indicate statistically significant
term changes in bloom frequency in the regions labelled in a and b, and their (P < 0.05) correlations. Maps created using ArcMap 10.4 and Python 3.8.
relationship to the SST and SST gradient. Linear slope (S) of bloom frequency

282 | Nature | Vol 615 | 9 March 2023

 Changes in climate can also affect ocean circulation, altering ocean HAEDAT. For example, in a recent global analysis of the HAEDAT events,
mixing and the transport of nutrients that drive the growth of marine Hallegraeff et al.6 report a dozen or more events per year for each of nine
phytoplankton and bloom formation30–32. We used the spatial SST gradi- regions over a 33-year study period, compared to the global median
ent (in °C m−1) as a proxy for the magnitude of oceanic mesoscale cur- bloom count of 4.3 in this study. There are several possible explana-
rents (the time-varying velocity of kinetic energy (also known as the eddy tions for this discrepancy, such as the many low-cell-concentration
kinetic energy (EKE))) by following the methods of a previous study33, HABs that are not detectable from space but that can still cause harm,
and examined its effects on algal blooms (Methods). The trend in the as well as sensor sensitivities and algorithm thresholds. Furthermore,
SST gradient appeared more spatially aligned to bloom frequency than our bloom count was averaged over all 1-km pixels within the EEZs,
SST. We found significant positive correlations (P < 0.05) between the whereas HAEDAT entries are based on discrete sampling regions. This
SST gradient and bloom frequency for various coastal current systems, underestimation does not, however, alter the trends and other conclu-
including the Canary (r = 0.84), Malvinas (r = 0.83), California (r = 0.81), sions of this study, as the metrics used here were constant across time
Benguela (r = 0.73), Gulf Stream (r = 0.61) and Oyashio (r = 0.58) currents. and space. Underestimates would have been uniform across regions
Trends in bloom frequency in subtropical eastern boundary upwelling globally. In this regard, it is of note that the study of Hallegraeff et al.6
regions (the California, Benguela and Canary currents) followed the found a significant link between the number of HAEDAT events over
changes in mesoscale currents (Fig. 3a,c). In the California Current Sys- time and the global expansion of aquaculture production, similar to
tem, the decrease in bloom frequency was probably due to the weakened findings in our study.
upwelling (represented by a reduced SST gradient and increased SST) The major contribution of our study is to provide a spatially and
and thus lower nutrient supply25. Conversely, the Canary and Benguela temporally consistent characterization of global coastal algal blooms
currents were characterized by strengthened upwelling and increased between 2003 and 2020. Globally, increasing trends in algal bloom
bloom frequency. The two western boundary current systems at high area and frequency are apparent. Regionally, however, trends were
latitudes (Malvinas and Oyashio)—although characterized by less pro- non-uniform owing to the compounded effects of changes in climate
nounced upwelling34—exhibited a similar mechanism to the subtropical (such as changes in SST and SST gradient and climate extremes),
eastern boundary regions. For the subtropical western boundary Gulf anthropogenic eutrophication and aquaculture development. Our
Stream current, the strengthened current jets (manifested as a larger SST daily mapping of bloom events offers valuable baseline information
gradient) brought more nutrients from the continental shelf35, trigger- to understand the mechanisms underlying the formation, mainte-
ing more algal blooms. Nevertheless, whether these changes in oceanic nance, and dissipation of algal blooms5,40. This could aid in developing
mesoscale activities were responses to wind, stratification, the shear of forecasting models (on either global or regional scales) that can help
ocean currents or other factors33 requires region-based investigations. minimize the consequences of harmful blooms, and can also help in
Global climate events, represented as the multivariate El Niño– policy decisions relating to the control of nutrient discharges and
Southern Oscillation index36 (MEI), also showed connections with coastal other HAB-stimulatory factors. Noting again that many blooms are
bloom frequency. The minimum MEI in 2010 (a strong La Niña year) beneficial, particularly in terms of their positive effects on ecosystems
was followed by a low bloom frequency in the following year, and the as well as on wild and farmed fisheries, the results here can also con-
largest MEI in 2015 (a strong El Niño year) was followed by the strongest tribute toward policies and management actions that sustain those
bloom frequency in 2016 (Fig. 2b and Extended Data Fig. 10a). beneficial blooms.
Changes in anthropogenic nutrient enrichment may have also con-
tributed to the trends in blooms37. For example, the decline in bloom
frequency in the Arabian Sea, without clear links to SST or SST gradient Online content
changes, could result from decreased fertilizer use in the surround- Any methods, additional references, Nature Portfolio reporting summa-
ing countries (such as Iran) (Extended Data Fig. 10). By contrast, the ries, source data, extended data, supplementary information, acknowl-
bloom strengthening in some Asian countries could be attributed to edgements, peer review information; details of author contributions
surges in fertilizer use38. We examined trends in fertilizer usage (either and competing interests; and statements of data and code availability
nitrogen or phosphorus) and bloom frequency and found high positive are available at https://doi.org/10.1038/s41586-023-05760-y.
correlations in China, Iran, Vietnam and the Philippines. Paradoxi-
cally, decreased fertilizer uses and increased bloom frequency were 1. Breitburg, D. et al. Declining oxygen in the global ocean and coastal waters. Science 359,
identified in some countries, suggesting that nutrient control efforts eaam7240 (2018).
2. Anderson, D. M. Turning back the harmful red tide. Nature 388, 513–514 (1997).
might have been counterbalanced by the stimulatory effects of climate 3. Beman, J. M., Arrigo, K. R. & Matson, P. A. Agricultural runoff fuels large phytoplankton
warming or other factors. Furthermore, the intensified aquaculture blooms in vulnerable areas of the ocean. Nature 434, 211–214 (2005).
industry in Finland, China, Algeria, Guinea, Vietnam, Argentina, Russia 4. Heisler, J. et al. Eutrophication and harmful algal blooms: a scientific consensus. Harmful
Algae 8, 3–13 (2008).
and Uruguay may also be linked to their increased bloom incidence, 5. Anderson, D. M., Cembella, A. D. & Hallegraeff, G. M. Progress in understanding harmful
as suggested by the significant positive correlations (r > 0.5, P < 0.05) algal blooms: paradigm shifts and new technologies for research, monitoring, and
between their aquaculture production and bloom frequency. A similar management. Annu. Rev. Mar. Sci. 4, 143–176 (2012).
6. Hallegraeff, G. M. et al. Perceived global increase in algal blooms is attributable to
relationship between aquaculture expansion and positive trends in HAB intensified monitoring and emerging bloom impacts. Commun. Earth Environ. 2, 117 (2021).
incidence was reported from an analysis of HAEDAT data6. However, 7. Smith, V. H. Eutrophication of freshwater and coastal marine ecosystems a global
analogous positive feedbacks for fertilizer or aquaculture were not problem. Environ. Sci. Pollut. Res. 10, 126–139 (2003).
8. Fleming, L. E. et al. Review of Florida red tide and human health effects. Harmful Algae
found in many other countries. Thus, an ecosystem model incorporat- 10, 224–233 (2011).
ing terrestrial and oceanic nutrient transport and nutrient–plankton 9. Richlen, M. L., Morton, S. L., Jamali, E. A., Rajan, A. & Anderson, D. M. The catastrophic
relationships of different species39 is required to quantify the contribu- 2008–2009 red tide in the Arabian Gulf region, with observations on the identification
and phylogeny of the fish-killing dinoflagellate Cochlodinium polykrikoides. Harmful
tions of natural and anthropogenic factors to algal blooms14. Algae 9, 163–172 (2010).
10. Hallegraeff, G. & Bolch, C. Unprecedented toxic algal blooms impact on Tasmanian
seafood industry. Microbiol. Aust. 37, 143–144 (2016).
11. Diaz, R. J. & Rosenberg, R. Spreading dead zones and consequences for marine ecosystems.
Future implications Science 321, 926–929 (2008).
We acknowledge that our criteria for a detectable bloom event is 12. Barton, A. D., Irwin, A. J., Finkel, Z. V. & Stock, C. A. Anthropogenic climate change drives
operationally defined by sensor sensitivities and other factors, and shift and shuffle in North Atlantic phytoplankton communities. Proc. Natl Acad. Sci. USA
113, 2964–2969 (2016).
that the bloom count metric used here may underestimate algal 13. Gobler, C. J. Climate change and harmful algal blooms: insights and perspective. Harmful
bloom incidence, particularly compared to harmful events entered in Algae 91, 101731 (2020).

Nature | Vol 615 | 9 March 2023 | 283

Article
14. Zohdi, E. & Abbaspour, M. Harmful algal blooms (red tide): a review of causes, impacts 31. Chelton, D. B., Gaube, P., Schlax, M. G., Early, J. J. & Samelson, R. M. The influence of
and approaches to monitoring and prediction. Int. J. Environ. Sci. Technol. 16, 1789–1806 nonlinear mesoscale eddies on near-surface oceanic chlorophyll. Science 334, 328–332
(2019). (2011).
15. Wells, M. L. et al. Future HAB science: directions and challenges in a changing climate. 32. Boyce, D. G., Petrie, B., Frank, K. T., Worm, B. & Leggett, W. C. Environmental structuring of
Harmful Algae 91, 101632 (2020). marine plankton phenology. Nat. Ecol. Evol. 1, 1484–1494 (2017).
16. Rabalais, N. N., Turner, R. E., Díaz, R. J. & Justić, D. Global change and eutrophication of 33. Martínez-Moreno, J. et al. Global changes in oceanic mesoscale currents over the satellite
coastal waters. ICES J. Mar. Sci. 66, 1528–1537 (2009). altimetry record. Nat. Clim. Change 11, 397–403 (2021).
17. Blondeau-Patissier, D., Gower, J. F., Dekker, A. G., Phinn, S. R. & Brando, V. E. A review of 34. Kämpf, J. & Chapman, P. in Upwelling Systems of the World 31–65 (Springer, 2016).
ocean color remote sensing methods and statistical techniques for the detection, 35. Lee, T. N., Yoder, J. A. & Atkinson, L. P. Gulf Stream frontal eddy influence on productivity
mapping and analysis of phytoplankton blooms in coastal and open oceans. Prog. of the southeast US continental shelf. J. Geophys. Res. 96, 22191–22205 (1991).
Oceanogr. 123, 123–144 (2014). 36. Wolter, K. & Timlin, M. S. El Niño/Southern Oscillation behaviour since 1871 as diagnosed
18. Wolny, J. L. et al. Current and future remote sensing of harmful algal blooms in the in an extended multivariate ENSO index (MEI. ext). Int. J. Climatol. 31, 1074–1087 (2011).
Chesapeake Bay to support the shellfish industry. Front. Mar. Sci. 7, 337 (2020). 37. Glibert, P. M. & Burford, M. A. Globally changing nutrient loads and harmful algal blooms:
19. Stumpf, R. P. et al. Monitoring Karenia brevis blooms in the Gulf of Mexico using satellite recent advances, new paradigms, and continuing challenges. Oceanography 30, 58–69
ocean color imagery and other data. Harmful Algae 2, 147–160 (2003). (2017).
20. Bernard, S., Kudela, R. M., Robertson Lain, L. & Pitcher, G. Observation of Harmful Algal 38. Lu, C. & Tian, H. Global nitrogen and phosphorus fertilizer use for agriculture production
Blooms with Ocean Colour Radiometry http://dx.doi.org/10.25607/OBP-1042 (IOCCG, 2021). in the past half century: shifted hot spots and nutrient imbalance. Earth Syst. Sci. Data 9,
21. Hu, C. et al. Red tide detection and tracing using MODIS fluorescence data: a regional 181–192 (2017).
example in SW Florida coastal waters. Remote Sens. Environ. 97, 311–321 (2005). 39. Falkowski, P. G., Barber, R. T. & Smetacek, V. Biogeochemical controls and feedbacks on
22. Andersen, J. H. et al. Long-term temporal and spatial trends in eutrophication status of ocean primary production. Science 281, 200–206 (1998).
the Baltic Sea. Biol. Rev. 92, 135–149 (2017). 40. Wells, M. L. et al. Harmful algal blooms and climate change: Learning from the past and
23. Gómez, F. & Boicenco, L. An annotated checklist of dinoflagellates in the Black Sea. present to forecast the future. Harmful Algae 49, 68–93 (2015).
Hydrobiologia 517, 43–59 (2004).
24. Townsend, D. W., Pettigrew, N. R. & Thomas, A. C. Offshore blooms of the red tide Publisher’s note Springer Nature remains neutral with regard to jurisdictional claims in
dinoflagellate, Alexandrium sp., in the Gulf of Maine. Cont. Shelf Res. 21, 347–369 (2001). published maps and institutional affiliations.
25. Pitcher, G. C., Figueiras, F. G., Hickey, B. M. & Moita, M. T. The physical oceanography of
upwelling systems and the development of harmful algal blooms. Prog. Oceanogr. 85,
5–32 (2010). Open Access This article is licensed under a Creative Commons Attribution
26. López-Cortés, D. J. et al. The state of knowledge of harmful algal blooms of Margalefidinium 4.0 International License, which permits use, sharing, adaptation, distribution
polykrikoides (a.k.a. Cochlodinium polykrikoides) in Latin America. Front. Mar. Sci. 6, 463 and reproduction in any medium or format, as long as you give appropriate
(2019). credit to the original author(s) and the source, provide a link to the Creative Commons licence,
27. Anderson, D. M. et al. Evidence for massive and recurrent toxic blooms of Alexandrium and indicate if changes were made. The images or other third party material in this article are
catenella in the Alaskan Arctic. Proc. Natl Acad. Sci. USA 118, e2107387118 (2021). included in the article’s Creative Commons licence, unless indicated otherwise in a credit line
28. Griffith, A. W., Doherty, O. M. & Gobler, C. J. Ocean warming along temperate western to the material. If material is not included in the article’s Creative Commons licence and your
boundaries of the Northern Hemisphere promotes an expansion of Cochlodinium intended use is not permitted by statutory regulation or exceeds the permitted use, you will
polykrikoides blooms. Proc. R. Soc. B 286, 20190340 (2019). need to obtain permission directly from the copyright holder. To view a copy of this licence,
29. Conley, D. J. Save the Baltic Sea. Nature 486, 463–464 (2012). visit http://creativecommons.org/licenses/by/4.0/.
30. Mahadevan, A., D’Asaro, E., Lee, C. & Perry, M. J. Eddy-driven stratification initiates North
Atlantic spring phytoplankton blooms. Science 337, 54–58 (2012). © The Author(s) 2023

284 | Nature | Vol 615 | 9 March 2023

Methods Oscillation activities and marine phytoplankton blooms. The dataset
is available from https://psl.noaa.gov/enso/mei/.
Data sources
MODIS on the Aqua satellite provides a global coverage within 1–2 Development of an automated bloom detection method
days. All images acquired by this satellite mission from January 2003 A recent study by the UNESCO Intergovernmental Oceanographic
to December 2020 were used in our study to detect global coastal Commission revealed that globally reported HAB events have
phytoplankton blooms, with a total of 0.76 million images. MODIS increased6. However, such an overall increasing trend was found to
Level-1A images were downloaded from the Ocean Biology Distributed be highly correlated with recently intensified sampling efforts6. Once
Active Archive Center (OB.DAAC) at NASA Goddard Space Flight Center this potential bias was accounted for by examining the ratio between
(GSFC), and were subsequently processed with SeaDAS software (ver- HAB events to the number of samplings5, there was no significant
sion 7.5) to obtain Rayleigh-corrected reflectance (Rrc (dimensionless), global trend in HAB incidence, though there were increases in certain
which was converted using the rhos (in sr−1) product (rhos × π) from regions. With synoptic, frequent, and large-scale observations, satel-
SeaDAS)41, remote sensing reflectance (Rrs (sr−1)) and quality control lite remote sensing has been extensively used to monitor algal blooms
flags (l2_flags). If a pixel was flagged by any of the following, it was then in oceanic environments17–19. For example, chlorophyll a (Chla) con-
removed from phytoplankton bloom detection: straylight, cloud, land, centrations, a proxy for phytoplankton biomass, has been provided as
high sunglint, high solar zenith angle and high sensor zenith angle a standard product by NASA since the proof-of-concept Coastal Zone
(https://oceancolor.gsfc.nasa.gov/atbd/ocl2flags/). MODIS level-3 Color Scanner (1978–1986) era43,44. The current default algorithm used
product for aerosol optical thicknesses (AOT) at 869 nm was also to retrieve Chla products is based on the high absorption of Chla at
obtained from OB.DAAC NASA GSFC (version R2018.0), which was the blue band45,46, which often shows high accuracy in the clear open
used to examine the impacts of aerosols on bloom trends. oceans but high uncertainties in coastal waters. This is because, in
We examined the algal blooms in the EEZs of 153 ocean-bordering productive and dynamic coastal oceans, the absorption of Chla in the
countries (excluding the EEZs in the Caspian Sea or around the blue band can be obscured by the presence of suspended sediments
Antarctic), 126 of which were found with at least one bloom in the past and/or coloured dissolved organic matter (CDOM)47. To address this
two decades. The EEZ dataset is available at https://www.marinere- problem, various regionalized Chla algorithms have been developed48.
gions.org/download_file.php?name=World_EEZ_v11_20191118.zip. Unfortunately, the concentrations of the water constituents (CDOM,
The EEZs are up to 200 nautical miles (or 370 km) away from coast- sediment and Chla) can vary substantially across different coastal
lines, which include all continental shelf areas and offer the majority of oceans. As a result, a universal Chla algorithm that can accurately
marine resources available for human use. Regional statistics of algal estimate Chla concentrations in global coastal oceans is not currently
blooms were also performed for LMEs. LMEs encompass global coastal available.
oceans and outer edges of coastal currents areas, which are defined Alternatively, many spectral indices have been developed to iden-
by various distinct features of the oceans, including hydrology, pro- tify phytoplankton blooms instead of quantifying their bloom bio-
ductivity, bathymetry and trophically dependent populations 42. mass, including the normalized fluorescence line height21 (nFLH),
Of the 66 LMEs identified globally, we excluded the Arctic and Ant- red tide index49 (RI), algal bloom index47 (ABI), red–blue difference
arctic regions and examined 54 LMEs. The boundaries of LMEs were (RBD)50, Karenia brevis bloom index50 (KBBI) and red tide detec-
obtained from https://www.sciencebase.gov/catalog/item/55c7772 tion index51 (RDI). In practice, the most important task for these
2e4b08400b1fd8244. index-based algorithms is to determine their optimal thresholds
We used HAEDAT to validate our satellite-detected phytoplankton for bloom classification. However, such optimal thresholds can be
blooms in terms of presence or absence. The HAEDAT dataset (http:// regional-or image-specific20, due to the complexity of optical fea-
haedat.iode.org) is a collection of records of HAB events, maintained tures in coastal waters and/or the contamination of unfavourable
under the UNESCO Intergovernmental Oceanographic Commission observational conditions (such as thick aerosols, thin clouds, and
and with data archives since 1985. For each HAB event, the HAEDAT so on), making it difficult to apply spectral-index-based algorithms
records its bloom period (ranging from days to months) and geoloca- at a global scale.
tion. We merged duplicate entries when both the recorded locations To circumvent the difficulty in determining unified thresholds for
and times of the HAEDAT events were very similar to one another, and a various spectral indices across global coastal oceans, an approach
total number of 2,609 HAEDAT events were ultimately selected between from a recent study to classify algal blooms in freshwater lakes52 was
2003 and 2020. adopted and modified here. In that study, the remotely sensed reflec-
We used the ¼° resolution National Oceanic and Atmospheric tance data in three visible bands (red, green and blue) were converted
Administration Optimum Interpolated SST (v. 2.1) data to examine the into two-dimensional colour space created by the Commission Interna-
potential simulating effects of warming on the global phytoplankton tionale del’éclairage (CIE), in which the position on the CIE chromaticity
trends. We also estimated the SST gradients following the method of diagram represented the colour perceived by human eyes (Extended
Martínez-Moreno33. As detailed in ref. 33, the SST gradient can be used Data Fig. 1a). As the algal blooms in freshwater lakes were manifested
as a proxy for the magnitude of oceanic mesoscale currents (EKE). as greenish colours, the reflectance of bloom-containing pixels was
We used the SST gradient to explore the effects of ocean circulation expected to be distributed in the green gamut of the CIE chromaticity
dynamics on algal blooms. diagram; the stronger the bloom, the closer the distance to the upper
Fertilizer uses and aquaculture production for different countries border of the diagram (the greener the water).
was used to examine the potential effects of nutrient enrichment Here, the colour of phytoplankton blooms in the coastal oceans
from humans on global phytoplankton bloom trends. Annual data can be greenish, yellowish, brownish, or even reddish53, owing to
between 2003 and 2019 on synthetic fertilizer use, including nitrogen the compositions of bloom species (diatoms or dinoflagellates) and
and phosphorus, are available from https://ourworldindata.org/fer- the concentrations of different water constituents. Furthermore, the
tilizers. Annual aquaculture production includes cultivated fish and Chla concentrations of the coastal blooms are typically lower than
crustaceans in marine and inland waters, and sea tanks, and the data those in inland waters, thus demanding more accurate classification
between 2003 and 2018 are available from https://ourworldindata. algorithms. Thus, the algorithm proposed by Hou et al.52 was modified
org/grapher/aquaculture-farmed-fish-production. when using the CIE chromaticity space for bloom detection in marine
The MEI, which combines various oceanic and atmospheric variables36, environments. Specifically, we used the following coordinate conver-
was used to examine the connections between El Niño–Southern sion formulas to obtain the xy coordinate values in the CIE colour space:
Article
x = X /(X + Y + Z ) We converted the Rrc data of 53,820 selected sample pixels into
y = Y /(X + Y + Z ) the xy coordinates in the CIE colour space (Extended Data Fig. 1a).
(1) As expected, these samples of bloom-containing pixels were located
X = 2.7689R + 1.7517G + 1.1302B
in the upper half of the chromaticity diagram (the green, yellow and
Y = 1.0000R + 4.5907G + 0.0601B
orange–red gamut) (Extended Data Fig. 1a). We determined the lower
Z = 0.0000R + 0.0565G + 5.5943B
boundary of these sample points in the chromaticity diagram, which
where R, G and B represent the Rrc at 748 nm, 678 nm (fluorescence represents the lightest colour and thus the weakest phytoplankton
band) and 667 nm in the MODIS Aqua data, respectively. By contrast, blooms; any point that falls above this boundary represents stronger
the R, G and B channels used in Hou et al.52 were the red, green and blue blooms. The method to determine the boundary was similar to Hou
bands. We used the fluorescence band for the G channel because, for a et al.52: we first binned the sample points according to the x value in
given region, the 678 nm signal increases monotonically with the Chla the chromaticity diagram and estimated the 1st percentile (Q1%) of
concentration for blooms of moderate intensity21, which is similar to the the corresponding Y for each bin; then, we fit the Q1% using two-order
response of greenness to freshwater algal blooms. Thus, the converted polynomial regression. Sensitivity analysis with Q0.3% (the three-sigma
y value in the CIE coordinate system represents the strength of the value) resulted in minor changes (<1%) in the resulting bloom areas for
fluorescence. In practice, for pixels with phytoplankton blooms, the single images. Notably, sample points were rarely located near white
converted colours in the chromaticity diagram will be located within points (x = 1/3 and y = 1/3, represent equal reflection from three RGB
the green, yellow or orange–red gamut (see Extended Data Fig. 1a); the bands) in the diagram, and we used two polynomial regressions to
stronger the fluorescence signal is, the closer the distance to the upper determine the boundaries for x values greater and less than 1/3, which
border of the CIE diagram (larger y value). By contrast, for bloom-free can be expressed as:
pixels without a fluorescence signal, their converted xy coordinates will
be located in the blue or purple gamut. Therefore, we can determine
1
y1 = 4.8093x 2 − 3.0958x + 0.8357 x < (3)
a lower boundary in the CIE two-dimensional coordinate system to 3
separate bloom and non-bloom pixels, similar to the method proposed
by Hou et al.52. 1
y2 = 4.9040x 2 − 3.5759x + 0.9862 x> (4)
We selected 53,820 bloom-containing pixels from the MODIS Rrc data 3
as training samples to determine the boundary of the CIE colour space.
These sample points were selected from nearshore waters worldwide Based on the above, if a pixel’s xy coordinate (converted from Rrc
where frequent phytoplankton blooms have been reported (Extended spectrum) satisfies the conditions of (x < 1/3 AND y > y1) or (x > 1/3 AND
Data Fig. 2); the algal species included various species of dinoflagellates y > y2), it is classified as a ‘bloom’ pixel.
and diatoms20. A total of 80 images was used, which were acquired from Depending on the local region and application purpose, the mean-
different seasons and across various bloom magnitudes, to ensure that ing of ‘phytoplankton bloom’ may differ. Here, for a global applica-
the samples used could almost exhaustively represent the different tion, the pixelwise bloom classification is based on the relationship
bloom conditions in the coastal oceans. (represented using the CIE colour space) between Rrc in the 667-, 678-
We combined the MODIS FLHRrc (fluorescence line height based on and 754-nm bands derived from visual interpretation of the 80 pairs
Rrc) and enhanced red–green–blue composite (ERGB) to delineate of FLHRrc and ERGB imagery. Instead of a simple threshold, we used a
bloom pixels manually. The FLHRrc image was calculated as: lower boundary of the sample points in the chromaticity diagram to
define a bloom. In simple words, a pixel is classified as a bloom if its
fluorescence signal is detectable (the associated xy coordinate in the
FLHRrc = Rrc678 × F678 − [Rrc667 × F667 + (Rrc748 × F748
(2) CIE colour space located above the lower boundary). Histogram of the
− Rrc667 × F667) × (678 − 667)/(748 − 667)] nFLH values from the 53,820 training pixels demonstrated the minimum
value of ~0.02 mW cm−2 μm−1 (Extended Data Fig. 1a), which is in line
where Rrc667, Rrc678 and Rrc748 are the Rrc at 667, 678 and 748 nm, respec- with the lower-bound signal of K. brevis blooms on the West Florida
tively, and F667, F678 and F748 are the corresponding extraterrestrial solar shelf21,47. Note that, such a minimum nFLH is determined from the global
irradiance. ERGB composite images were generated using Rrc of three training pixels, and it does not necessarily represent a unified lower
bands at 555 (R), 488 (G) and 443 nm (B). Although phytoplankton-rich bound for phytoplankton blooms across the entire globe, especially
and sediment-rich waters have high FLHRrc values, they appear as dark- considering that fluorescence efficiency may be a large variable across
ish and bright features in the ERGB images (Extended Data Fig. 3), different regions. Different regions may have different lower bounds
respectively21. In fact, visual examination with fluorescence signals of nFLH to define a bloom, and such variability is represented by the
and ERGB has been widely accepted as a practical way to deline- predefined boundary in the CIE chromaticity diagram in our study. Cor-
ate coastal algal blooms on a limited number of images21,54,55. Note respondingly, although the accuracy of Chla retrievals may have large
that the FLHRrc here was slightly different from the NASA standard uncertainties in coastal waters, the histogram of the 53,820 training
nFLH product56, as the latter is generated using Rrs (corrected for both pixels shows a lower bound of ~1 mg m−3 (Extended Data Fig. 1a). Simi-
Rayleigh and aerosol scattering) instead of Rrc (with residual effects of larly to nFLH, such a lower bound may not be applicable to all coastal
aerosols). However, when using the NASA standard algorithm to further regions, as different regions may have different lower bounds of Chla
perform aerosol scattering correction over Rrc, 20.7% of our selected for bloom definition.
bloom-containing pixels failed to obtain valid Rrs (without retrievals Although the MODIS cloud (generated by SeaDAS with Rrc869 < 0.027)
or flagged as low quality), especially for those with strong blooms (see and associated straylight flags can be used to exclude most clouds, we
examples in Extended Data Fig. 4). Likewise, we also found various found that residual errors from thin clouds or cloud shadows could
nearshore regions with invalid Rrs retrievals. By contrast, Rrc had valid affect the spectral shape and cause misclassification for bloom detec-
data for all selected samples and showed more coverage in nearshore tions. Thus, we designed two spectral indices to remove such effects:
coastal waters. The differences between Rrs and Rrc were because the
assumptions for the standard atmospheric correction algorithm do Index1 = nRrc488 − nRrc443 − (nRrc555 − nRrc443) × 0.5 (5)
not hold for bloom pixels or nearshore waters with complex optical
properties57. In fact, Rrc has been used as an alternative to Rrs in various Index2 = nRrc555 − nRrc469 − (nRrc645 − nRrc469) × 0.5 (6)
applications in complex waters58,59.
where Index1 and Index2 were used to remove cloud shadows and estimated the ratio between the number of satellite images with ‘bloom
clouds, respectively. The nRrc443, nRrc488 and nRrc555 in index1 are the detected’ against the number of valid images (see definition above)
normalized Rrc, obtained by normalizing Rrc488. Similar calculations during the bloom periods across the entire globe (Supplementary
were performed for index2. The purpose of normalizations is to elimi- Table 1). Moreover, we calculated the number of events with at least
nate the effect of the absolute magnitude of the reflectance, so that one successful satellite bloom detection (Ns), and then estimated the
the thresholds of these two indices are influenced by only the relative ratio between Ns and the total HAB events for each year. Results showed
magnitude (spectral shape). We determined thresholds for Index1 that substantial amounts (averaged at 51.2%) of satellite observations
(>0.12) and Index2 (<0.012) through trial-and-error and ensured that acquired during the HAB event periods were found with phytoplank-
the misclassifications caused by residual errors from clouds and cloud ton blooms by our algorithm. Overall, 79.3% of the global HAB events
shadows could be effectively removed. After applying the cloud/cloud were successfully identified by satellite. The discrepancies between
shadow and various other masks that are associated with l2_flags, we satellite and in situ observations could be explained by the following
obtained an annual mean valid pixel observation (Nvobs) of ~2.0 × 105 reasons: first, our study focused only on the phytoplankton blooms
for global 1° × 1° grid cells, and the fluctuation patterns and trends of that are resolvable by satellite fluorescence signals; other types of HAB
Nvobs, either annually or seasonally, are different from that of the global events in the HAEDAT dataset may not have been detectable by satellite
bloom frequency and affected areas (see Supplementary Fig. 1). observations, such as events with lower phytoplankton biomass but
high toxicity, occurrences at the subsurface layers, or fluorescence
Assessments of the algorithm performance signals overwhelmed by suspended sediments65–67. Second, although
In addition to phytoplankton blooms, macroalgal blooms (Sargassum the HAEDAT recorded HAB events could be sustained for long periods,
and Ulva) frequently occur in many coastal oceans60–63. To verify high biomass of surface algae may not have occurred every day within
whether our CIE-fluorescence algorithm could eliminate such impacts, this period due to the changes in stratification, precipitation, wind, ver-
we compared the spectra between micro-and macroalgal blooms (see tical migration of cells, and many other factors68. Third, the spatial scale
Extended Data Fig. 1b). We found that the spectral shapes of Sargassum of certain HAB events may have been too small to be identified using
and Ulva are substantially different from those of microalgae, particu- the 1-km resolution MODIS observations. Fourth, a reduced MODIS
larly for the three bands used for CIE coordinate conversion. The con- satellite observation frequency by the contaminations of clouds and
verted xy coordinates for macroalgae were located in the purple–red land adjacency effects69. Therefore, we believe the underestimations
gamut of the CIE diagram, which was far below the predefined bound- of satellite-detected blooms compared to the in situ reported HAB
ary (Extended Data Fig. 1). Moreover, our algorithm is not affected by events were mainly due to inconsistencies between the two observa-
highly turbid waters for the following two reasons: first, extremely high tions rather than uncertainties in our algorithm.
turbidity tends to saturate the MODIS ocean bands64, which can be eas- Because Rrc depends not only on water colour but also on aerosols
ily excluded; second, without a fluorescence peak, the reflectance of (type and concentration) and solar and viewing geometry, a sensitivity
unsaturated turbid waters, after conversion to CIE coordinates, will be analysis was used to determine whether such variables could impact
located below the predefined boundary (see example in Extended Data bloom detection. Aerosol reflectance (ρa) with different AOTs at 869
Fig. 1b). We also confirmed that the spectral shapes of coccolithophore nm was simulated using the NASA-recommended maritime aerosol
blooms are different from dinoflagellates and diatoms (see example in model (r75f02, with a relative humidity of 75% and a fine-mode frac-
Extended Data Fig. 1b), and thus they are excluded from our algorithm. tion of 2%). Then, ρa of each MODIS band was added to Rrc images, and
Three different types of validation methods were adopted to dem- the resulting bloom areas with and without added ρa were compared.
onstrate the reliability of the proposed CIE-fluorescence algorithm for Results showed that even with a change of 0.02 in AOT at 869 nm, the
phytoplankton bloom detection in global coastal oceans, including bloom areas showed minor changes (<2%) in the tested images; minor
visual inspections of the RGB, ERGB and FLHRrc images, verifications changes were also found when we used different aerosol models to
using independent manually delineated algal blooms, and comparisons conduct ρa simulations70. Note that 0.02 represents the high end of
with the reported HAB events from the HAEDAT dataset. the AOT intra-annual variability in coastal oceans (see Extended Data
First, we selected MODIS Aqua images from different locations where Fig. 5), and the associated interannual changes are much smaller.
coastal phytoplankton blooms have been recorded in the published Thus, the use of Rrc instead of the fully atmospherically corrected
literature. We visually compared the RGB, ERGB, and FLHRrc images, and reflectance Rrs could have limited impacts on our detected global
our algorithm detected bloom patterns (see examples in Extended Data bloom trend.
Fig. 3). Comparisons from various images worldwide showed that our We also tried various index-based algorithms developed in previ-
algorithm could successfully identify regions with high FLHRrc values ous studies. However, results showed that all these methods require
and brownish-to-darkish features on the ERGB images, which can be image-specific thresholds to accurately determine algal bloom bounda-
considered phytoplankton blooms. ries for different coastal regions (see Extended Data Fig. 6). By con-
Second, we delineated additional 15,466 bloom-containing pixels trast, although our CIE-fluorescence algorithm may lead to different
from 35 images covering different coastal areas, using the same visual bloom thresholds for different regions, it can identify bloom pixels
inspection and manual delineation method as for the training sample without adjusting the coefficients and, therefore, is more suitable for
pixels. Moreover, we also selected 14,149 bloom-free pixels (bright or global-scale bloom assessment efforts.
turquoise green colours on ERGB images or low FLHRrc values) within We acknowledge that our satellite-detected algal blooms represent
the same images as bloom-containing images. We applied our algo- only high amounts of phytoplankton biomass on the ocean surfaces
rithm to all these pixels, and compared the algorithm-identified and without distinguishing whether such blooms produce toxins or are
manually delineated results. Our CIE-fluorescence algorithm showed harmful to marine environments. Furthermore, with only limited
high values in both producer and user accuracies (92.04% and 98.63%) spectral information from MODIS, it is difficult to discriminate the
(Supplementary Table 1), and appeared effective at identifying bloom phytoplankton species of algal blooms; such information could help
pixels and excluding false negatives (blooms classified as non-blooms) to improve our understanding of the impacts of these phytoplankton
and false positives (non-blooms classified as blooms). blooms. However, we expect these challenges to be addressed soon with
Third, we validated the satellite-detected phytoplankton blooms the scheduled launch of the Plankton, Aerosol, Cloud, ocean Ecosystem
using in situ reported HAB events from the HAEDAT dataset. For each (PACE) mission by NASA in 2024, where the hyperspectral measure-
HAB event in the HAEDAT dataset, we obtained all MODIS images over ments over a broad spectrum (350–885 nm) will make species-level
the reported bloom period (from days to months). Within each year, we classifications possible71.
Article
was set from day 150 to day 270 within the year. We further found that
Exploring the patterns and trends of global coastal the TMBAA for cluster II showed small changes over the entire period
phytoplankton blooms (Extended Data Fig. 9d), indicating relatively stable phenological cycles
We applied the CIE-fluorescence algorithm to all MODIS Aqua level-2 Rrc for those phytoplankton blooms32,77.
images, and a total number of 0.76 million images between 2003 and
2020 were processed. We mapped all detected blooms into 1-km daily Reporting summary
scale level-3 composites. The number of bloom counts within a year Further information on research design is available in the Nature Port-
for each location can be easily enumerated, and the long-term annual folio Reporting Summary linked to this article.
mean values were then estimated (Fig. 1a). We further calculated the
total global bloom-affected area (the areas where algal blooms were
detected at least once) for each year and examined their changes over Data availability
time (Fig. 2b). The satellite-based dataset of global coastal algal bloom at 1-km resolu-
We defined bloom frequency (dimensionless) to represent the den- tion and the associated code are available at https://doi.org/10.5281/
sity of phytoplankton blooms for a year by integrating the bloom count zenodo.7359262. Source data are provided with this paper.
and bloom-affected areas within 1°×1° grid cells within that year, which
is expressed as: 41. Hu, C. A novel ocean color index to detect floating algae in the global oceans. Remote
Sens. Environ. 113, 2118–2129 (2009).
n 42. Sherman, K. Adaptive management institutions at the regional level: the case of large
n
Bloom frequency =
N ∑ Mi (7) 43.
marine ecosystems. Ocean Coast. Manag. 90, 38–49 (2014).
Gordon, H. R., Clark, D. K., Mueller, J. L. & Hovis, W. A. Phytoplankton pigments from the
i =1 Nimbus-7 coastal zone color scanner: comparisons with surface measurements. Science
210, 63–66 (1980).
where Mi is the enumerated bloom count for each 1-km resolution pixel 44. Moore, J. K. & Abbott, M. R. Phytoplankton chlorophyll distributions and primary
production in the Southern Ocean. J. Geophys. Res. 105, 28709–28722 (2000).
in a year within one 1° × 1° grid cell, and n represents the associated
45. Hu, C. et al. Improving satellite global chlorophyll a data products through algorithm
number of bloom-affected pixels in the same cell (the number of pixels refinement and data recovery. J. Geophys. Res. 124, 1524–1543 (2019).
with Mi > 0), and N is the total number of 1-km MODIS pixels in this grid 46. Hu, C., Lee, Z. & Franz, B. Chlorophyll a algorithms for oligotrophic oceans: a novel
approach based on three-band reflectance difference. J. Geophys. Res. 117, C01011
cell. We estimated the bloom frequency for each year between 2003 and
(2012).
2020, and determined the long-term trend over global EEZs through a 47. Hu, C. & Feng, L. Modified MODIS fluorescence line height data product to improve
linear least-squares regression (see Fig. 2a). image interpretation for red tide monitoring in the eastern Gulf of Mexico. J. Appl. Remote
Sens. 11, 012003 (2016).
Continental and country-level statistics were performed for bloom
48. Siswanto, E. et al. Empirical ocean-color algorithms to retrieve chlorophyll-a, total
count, bloom-affected areas, and bloom frequency (Fig. 1b,c and suspended matter, and colored dissolved organic matter absorption coefficient in the
Supplementary Table 2), using boundaries for the EEZs of different Yellow and East China Seas. J. Oceanogr. 67, 627–650 (2011).
49. Ahn, Y.-H. & Shanmugam, P. Detecting the red tide algal blooms from satellite ocean
ocean-bordering countries (see above). Similar statistics were also con-
color observations in optically complex Northeast-Asia coastal waters. Remote Sens.
ducted for 54 LMEs (Extended Data Fig. 7 and Supplementary Table 3). Environ. 103, 419–437 (2006).
50. Amin, R. et al. Novel optical techniques for detecting and classifying toxic dinoflagellate
Correlations with SST and SST gradient Karenia brevis blooms using satellite imagery. Opt. Express 17, 9126–9144 (2009).
51. Shen, F., Tang, R., Sun, X. & Liu, D. Simple methods for satellite identification of algal
To assess the impacts of climate change on long-term trends in coastal blooms and species using 10-year time series data from the East China Sea. Remote Sens.
phytoplankton blooms, we correlated the annual mean bloom fre- Environ. 235, 111484 (2019).
52. Hou, X. et al. Global mapping reveals increase in lacustrine algal blooms over the past
quency and the associated SST and SST gradient in various coastal decade. Nat. Geosci. 15, 130–134 (2022).
current systems for grid cells with significant changes in bloom fre- 53. Dierssen, H. M., Kudela, R. M., Ryan, J. P. & Zimmerman, R. C. Red and black tides:
quency (Fig. 3c). The SST and SST gradient were averaged over the quantitative analysis of water-leaving radiance and perceived color for phytoplankton,
colored dissolved organic matter, and suspended sediments. Limnol. Oceanogr. 51,
growth window within a year, assuming that the changes within the 2646–2659 (2006).
growth window, either in water temperatures or ocean circulations, 54. Zhao, J., Temimi, M. & Ghedira, H. Characterization of harmful algal blooms (HABs) in the
play more important roles in the bloom trends compared to other Arabian Gulf and the Sea of Oman using MERIS fluorescence data. ISPRS J. Photogramm.
Remote Sens. 101, 125–136 (2015).
seasons32. 55. Qi, L. et al. Noctiluca blooms in the East China Sea bounded by ocean fronts. Harmful
We determined the growth window of phytoplankton blooms for Algae 112, 102172 (2022).
each 1° × 1° grid cell (Extended Data Fig. 9a) using the following method: 56. Behrenfeld, M. J. et al. Satellite-detected fluorescence reveals global physiology of
ocean phytoplankton. Biogeosciences 6, 779–794 (2009).
first, we estimated the proportion of cumulative bloom-affected pixels 57. Gordon, H. R. Atmospheric correction of ocean color imagery in the Earth Observing
within the grid cells for a year. Second, a generalized additive model72 System era. J. Geophys. Res. 102, 17081–17106 (1997).
was used to determine the shape of the phenological curves (Extended 58. Feng, L., Hou, X., Li, J. & Zheng, Y. Exploring the potential of Rayleigh-corrected reflectance
in coastal and inland water applications: a simple aerosol correction method and its merits.
Data Fig. 9b), where a log link function and a cubic cyclic regression ISPRS J. Photogramm. Remote Sens. 146, 52–64 (2018).
spline smoother were applied73,74. Third, the timing of maximum 59. Feng, L., Hu, C., Han, X., Chen, X. & Qi, L. Long-term distribution patterns of chlorophyll-a
bloom-affected areas (TMBAA) was then determined by identifying concentration in China’s largest freshwater lake: MERIS full-resolution observations with a
practical approach. Remote Sens. 7, 275–299 (2015).
the inflection point on the bloom growth curve (Extended Data Fig. 9c). 60. Xiao, J. et al. An anomalous bi-macroalgal bloom caused by Ulva and Sargassum
To facilitate comparisons across Northern and Southern Hemispheres, seaweeds during spring to summer of 2017 in the western Yellow Sea, China. Harmful
the year in the Southern Hemisphere was shifted forward by 183 days Algae 93, 101760 (2020).
61. Teichberg, M. et al. Eutrophication and macroalgal blooms in temperate and tropical
(Extended Data Fig. 9c). We characterized the similarity of the bloom coastal waters: nutrient enrichment experiments with Ulva spp. Glob. Change Biol. 16,
growth curve between different grid cells and grouped them into three 2624–2637 (2010).
distinct clusters using a fuzzy c-means cluster analysis method75,76. 62. Viaroli, P. et al. Nutrient and iron limitation to Ulva blooms in a eutrophic coastal lagoon
(Sacca di Goro, Italy). Hydrobiologia 550, 57–71 (2005).
We found uniform distributions of the clusters over large geographic 63. Dierssen, H. M., Chlus, A. & Russell, B. Hyperspectral discrimination of floating mats of
areas. Cluster I is mainly distributed in mid-low latitudes (<45° N and seagrass wrack and the macroalgae Sargassum in coastal waters of Greater Florida Bay
<30° S), where the maximum bloom-affected areas were expected in the using airborne remote sensing. Remote Sens. Environ. 167, 247–258 (2015).
64. Hu, C. et al. Dynamic range and sensitivity requirements of satellite ocean color sensors:
early period of the year. Cluster II was mostly found in higher latitudes, learning from the past. Appl. Opt. 51, 6045–6062 (2012).
with bloom developments (quasi-) synchronized with increases in SST. 65. Trainer, V. L. et al. Pelagic harmful algal blooms and climate change: lessons from
Cluster III was detected along the coastlines, where the bloom-affected nature’s experiments with extremes. Harmful Algae 91, 101591 (2020).
66. Mardones, J. I. et al. Disentangling the environmental processes responsible for the
areas increase throughout the entire year. In practice, the growth win- world’s largest farmed fish-killing harmful algal bloom: Chile, 2016. Sci. Total Environ.
dow for clusters I and III was set as the entire year, and that for cluster II 766, 144383 (2021).
67. Gilerson, A. et al. Fluorescence component in the reflectance spectra from coastal Acknowledgements We thank NASA for providing global MODIS satellite images, and the
waters. II. Performance of retrieval algorithms. Opt. Express 16, 2446–2460 (2008). Intergovernmental Oceanographic Commission (IOC) of UNESCO for providing the HAEDAT
68. Lee, J. H., Harrison, P. J., Kuang, C. & Yin, K. in The Environment in Asia Pacific Harbours dataset. L.F. was supported by the National Natural Science Foundation of China (no.
(ed. Wolanski, E.) 187–206 (Springer, 2006). 41890852, 42271322 and 41971304). D.M.A. was supported by the Woods Hole Center for
69. Feng, L. & Hu, C. Cloud adjacency effects on top-of-atmosphere radiance and ocean Oceans and Human Health (National Science Foundation grant OCE-1840381 and National
color data products: a statistical assessment. Remote Sens. Environ. 174, 301–313 (2016). Institutes of Health grants NIEHS-1P01-ES028938-01).
70. Ahmad, Z. et al. New aerosol models for the retrieval of aerosol optical thickness and
normalized water-leaving radiances from the SeaWiFS and MODIS sensors over coastal Author contributions Y.D. and S.Y.: methodology, data processing and analyses, and writing.
regions and open oceans. Appl. Opt. 49, 5545–5560 (2010). L.F.: conceptualization, methodology, funding acquisition, supervision and writing. D.Z.: data
71. Werdell, P. J. et al. The Plankton, Aerosol, Cloud, Ocean Ecosystem mission: status, processing and analysis. C.H., W.X., D.M.A., Y.L., X.-P.S., D.G.B., L.G. and C.Z. participated in
science, advances. Bull. Am. Meteorol. Soc. 100, 1775–1794 (2019). interpreting the results and refining the manuscript.
72. Hastie, T. J. Generalized Additive Models (Routledge, 2017).
73. Wood, S. N. Generalized Additive Models: An Introduction with R (CRC press, 2017). Competing interests The authors declare no competing interests.
74. Macgregor, C. J. et al. Climate-induced phenology shifts linked to range expansions in
species with multiple reproductive cycles per year. Nat. Commun. 10, 4455 (2019). Additional information
75. Bezdek, J. C. Pattern Recognition with Fuzzy Objective Function Algorithms (Springer, Supplementary information The online version contains supplementary material available at
2013). https://doi.org/10.1038/s41586-023-05760-y.
76. Bi, S. et al. Optical classification of inland waters based on an improved fuzzy C-means Correspondence and requests for materials should be addressed to Lian Feng.
method. Opt. Express 27, 34838–34856 (2019). Peer review information Nature thanks Bryan Franz, Bingkun Luo and the other, anonymous,
77. Kheireddine, M., Mayot, N., Ouhssain, M. & Jones, B. H. Regionalization of the Red Sea reviewer(s) for their contribution to the peer review of this work. Peer reviewer reports are
based on phytoplankton phenology: a satellite analysis. J. Geophys. Res. 126, available.
e2021JC017486 (2021). Reprints and permissions information is available at http://www.nature.com/reprints.
Article

Extended Data Fig. 1 | Development of the CIE-fluorescence algorithm to selected spectra for phytoplankton blooms, macroalgal blooms (Ulva and
detect phytoplankton blooms using MODIS satellite imagery. (a). A1: The Sargassum), coccolithophore blooms, and sediment-rich turbid waters. The x-y
density plot of manually delineated bloom-containing pixels in the CIE numbers indicate their corresponding positions in the CIE coordinate system.
coordinate system (n = 53,820), and their distribution in the CIE color space The black rectangular boxes in the three lower panels highlight different
(box in A2). A3: Histograms of nFLH and Chla for the delineated pixels, obtained spectral shapes between phytoplankton blooms and other features near the
using NASA standard algorithms 47,57. (b) MODIS true color composites and fluorescence band. Maps created using ArcMap 10.4.
Extended Data Fig. 2 | MODIS-detected bloom count within certain years Database (HAEDAT) within the same year. The lower right panel shows the
for several coastal regions with frequently reported blooms. The MODIS locations of all the HAEDAT records that were used for algorithm validations in
observational year is annotated within each panel, and overlaid points indicate this study (Supplementary Table 1), which also demonstrates the increase in
in situ recorded harmful algal bloom events from the Harmful Algae Event sampling effort in the most recent years. Created using ArcMap 10.4.
Article

Extended Data Fig. 3 | Performance of the CIE-fluorescence algorithm for bloom area (green pixels) detected by the CIE-fluorescence algorithm. Created
phytoplankton bloom detection in 12 selected coastal oceans. From left using ArcMap 10.4.
to right are the RGB-true color composite, ERGB composite, FLHRrc, and the
Extended Data Fig. 4 | Examples showing disadvantages of using NASA amounts of invalid Rrs retrievals can be observed in the red-encircled areas in
standard R rs (i.e., with the removal of both Rayleigh and aerosol scattering) which severe blooms can be found. Additionally, nFLH shows high values at
in algal bloom detection. From left to right are the RGB composites, ERGB, cloud edges (yellow-encircled areas), making it challenging to use a simple
nFLH, and the bloom areas (green pixels) detected by the CIE-fluorescence threshold to classify blooms. However, such problems can be circumvented in
algorithm (based on Rrc, without the removal of aerosol scattering). Substantial our CIE-fluorescence algorithm. Created using ArcMap 10.4.
Article

Extended Data Fig. 5 | Sensitivity analysis of the impacts of aerosols on show the RGB composites, and the right three columns show the bloom areas
bloom detection. (a) Responses of bloom area (BA) to changes in aerosol under different AOTs. The percentages of BA changes are annotated in the
optical thickness (AOT). Aerosol reflectance (ρa) with AOTs of 0.01 and 0.02 at panels. (b) The standard deviation between the 12 monthly mean values of AOT
869-nm is simulated and added to the MODIS images, and the resulting bloom in global coastal waters (i.e., 66.7% of the intra-annual variability), and the
areas (green pixels) with and without added ρa are compared. The left columns histogram is shown in (c). Maps created using ArcMap 10.4.
Extended Data Fig. 6 | Comparison of different index-based algorithms in bloom areas (i.e., high nFLH values, which appear as bright and darkish features
algal bloom detection in various coastal regions. Image-specific thresholds on the ERGB images). The left panels are the bloom areas (green pixels)
(annotated within the panels) are required (labeled within the panels) for RI50, extracted using our CIE-fluorescence algorithm. The RGB-true color and ERGB
ABI (estimated with FLHRrc)48, RBD51, KBBI51, and RDI52 to delineate accurate composites are shown in Extended Data Fig. 3. Created using ArcMap 10.4.
Article

Extended Data Fig. 7 | Annual median bloom count and the proportion of ordered from the largest to the smallest. The LMEs are grouped by continent,
bloom-affected areas for large marine ecosystems (LMEs). (a) Annual and their names, numbers, and locations are shown in (a) and (b). Map created
median bloom count, (b) proportion of bloom-affected areas. The data are using Python 3.8.
Extended Data Fig. 8 | Comparison of bloom counts in the estuarine and
non-estuarine regions. Boxplots for long-term mean bloom count in the
estuarine (n = 13,622 pixel observations) and non-estuarine (n = 361,604 pixel
observations) regions. Comparison analysis was performed by two sided
Welch’s t-test (P < 0.001).Upper and lower bounds are first and third quartiles,
the bar in the middle represents the median value, and the whiskers show the
minimum and maximum values. Sixty-two estuarine zones from large rivers
were selected, and the boundary of each zone was manually delineated
according to high-resolution satellite images.
Article

Extended Data Fig. 9 | Clusters of different bloom growth paths. (a) The deviations are shown with dashed lines and gray shading. The proportions of
spatial distribution of different clusters. The fractions of different clusters different clusters in the global bloom-affected areas are annotated. (c) and (f)
across different latitudes are summarized. (b) The development of the The mean timing of the maximum bloom-affected areas (TMBAA) and the
maximum bloom-affected areas within a year within 1° × 1° grid cells, where associated standard deviations between 2003 and 2019. The whole year in the
all global grid cells are grouped into three distinct clusters according to the Southern Hemisphere is shifted forward by 183 days in (c). Maps created using
similarity of the bloom growth curve. The colored bond curves represent the Python 3.8.
mean values of all the grid cells, and their mean SST and associated standard
Extended Data Fig. 10 | Changes in climate extremes, global fertilizer uses, and La Niña events, respectively. The dots show annual mean values.
and fishery production over the past two decades. (a) Changes in the (b–c) Trends of nitrogen and phosphorus from 2003 to 2019 for different
bi-monthly Multivariate El Niño–Southern Oscillation (ENSO) index (MEI) countries. (d) Trends of fishery production from 2003 to 2018. Gray indicates
between 2002 and 2020. Positive and negative MEI values represent EI Niño no data. Maps created using ArcMap 10.4.
 nature portfolio | reporting summary
Corresponding author(s): Lian Feng
Last updated by author(s): Nov 23, 2022

Reporting Summary
Nature Portfolio wishes to improve the reproducibility of the work that we publish. This form provides structure for consistency and transparency
in reporting. For further information on Nature Portfolio policies, see our Editorial Policies and the Editorial Policy Checklist.

Statistics
For all statistical analyses, confirm that the following items are present in the figure legend, table legend, main text, or Methods section.
n/a Confirmed
The exact sample size (n) for each experimental group/condition, given as a discrete number and unit of measurement
A statement on whether measurements were taken from distinct samples or whether the same sample was measured repeatedly
The statistical test(s) used AND whether they are one- or two-sided
Only common tests should be described solely by name; describe more complex techniques in the Methods section.

A description of all covariates tested

A description of any assumptions or corrections, such as tests of normality and adjustment for multiple comparisons
A full description of the statistical parameters including central tendency (e.g. means) or other basic estimates (e.g. regression coefficient)
AND variation (e.g. standard deviation) or associated estimates of uncertainty (e.g. confidence intervals)

For null hypothesis testing, the test statistic (e.g. F, t, r) with confidence intervals, effect sizes, degrees of freedom and P value noted
Give P values as exact values whenever suitable.

For Bayesian analysis, information on the choice of priors and Markov chain Monte Carlo settings
For hierarchical and complex designs, identification of the appropriate level for tests and full reporting of outcomes
Estimates of effect sizes (e.g. Cohen's d, Pearson's r), indicating how they were calculated
Our web collection on statistics for biologists contains articles on many of the points above.

Software and code

Policy information about availability of computer code
Data collection The satellite data were obtained from the U.S. National Aeronautics and Space Administration (NASA) Goddard Space Flight Center (GSFC).

Data analysis SeaDAS (Version 7.5) were used to analyze the satellite images.
For manuscripts utilizing custom algorithms or software that are central to the research but not yet described in published literature, software must be made available to editors and
reviewers. We strongly encourage code deposition in a community repository (e.g. GitHub). See the Nature Portfolio guidelines for submitting code & software for further information.

Data
Policy information about availability of data
All manuscripts must include a data availability statement. This statement should provide the following information, where applicable:
March 2021

- Accession codes, unique identifiers, or web links for publicly available datasets
- A description of any restrictions on data availability
- For clinical datasets or third party data, please ensure that the statement adheres to our policy

The MODIS Aqua data can be obtained from the U.S. National Aeronautics and Space Administration (NASA) Goddard Space Flight Center (GSFC).
The in situ reported HAB data are available from events from http://haedat.iode.org.
The Exclusive economic zones (EEZs) dataset is available at https://www.marineregions.org/download_file.php?name=World_EEZ_v11_20191118.zip.
The boundaries of large marine ecosystems (LMEs) were obtained from https://www.sciencebase.gov/catalog/item/55c77722e4b08400b1fd8244.

1
Annual data between 2003 and 2019 on synthetic fertilizer use, including nitrogen and phosphorus, are available from https://ourworldindata.org/fertilizers.
Annual aquaculture production includes cultivated fish and crustaceans in marine and inland waters, and sea tanks, and the data between 2003 and 2018 are

nature portfolio | reporting summary

available from https://ourworldindata.org/grapher/aquaculture-farmed-fish-production.
The dataset is available from https://psl.noaa.gov/enso/mei/.

Human research participants

Policy information about studies involving human research participants and Sex and Gender in Research.

Reporting on sex and gender Use the terms sex (biological attribute) and gender (shaped by social and cultural circumstances) carefully in order to avoid
confusing both terms. Indicate if findings apply to only one sex or gender; describe whether sex and gender were considered in
study design whether sex and/or gender was determined based on self-reporting or assigned and methods used. Provide in the
source data disaggregated sex and gender data where this information has been collected, and consent has been obtained for
sharing of individual-level data; provide overall numbers in this Reporting Summary. Please state if this information has not
been collected. Report sex- and gender-based analyses where performed, justify reasons for lack of sex- and gender-based
analysis.

Population characteristics Describe the covariate-relevant population characteristics of the human research participants (e.g. age, genotypic
information, past and current diagnosis and treatment categories). If you filled out the behavioural & social sciences study
design questions and have nothing to add here, write "See above."

Recruitment Describe how participants were recruited. Outline any potential self-selection bias or other biases that may be present and
how these are likely to impact results.

Ethics oversight Identify the organization(s) that approved the study protocol.

Note that full information on the approval of the study protocol must also be provided in the manuscript.

Field-specific reporting
Please select the one below that is the best fit for your research. If you are not sure, read the appropriate sections before making your selection.

Life sciences Behavioural & social sciences Ecological, evolutionary & environmental sciences
For a reference copy of the document with all sections, see nature.com/documents/nr-reporting-summary-flat.pdf

Ecological, evolutionary & environmental sciences study design

All studies must disclose on these points even when the disclosure is negative.

Study description This study developed a novel method to map global coastal algal blooms and used this tool to examine satellite images between
2003 and 2020, addressing three fundamental questions: 1) where and how frequently have global coastal oceans been affected by
phytoplankton blooms? 2) have the blooms expanded or intensified over the past two decades, both globally and regionally? and 3)
what are the potential drivers?

Research sample Three separate samples were selected. 1) MODIS Aqua images were used to develop the phytoplankton bloom extraction algorithm,
2) MODIS Aqua images and were used to verify the reliability of the algorithm and the accuracy of the phytoplankton bloom
extraction results, and 3) in situ reported HAB events from the HAEDAT dataset were used to validate the accuracy of the
phytoplankton bloom extraction results.

Sampling strategy A total of 115 MODIS Aqua images were selected from the different locations where coastal phytoplankton blooms have been
recorded in the published literature, of which 80 were used for algorithm development and 35 were used for algorithm validation. A
total number of 2609 HAB events that occurred in the coastal area were selected from the HAEDAT dataset.

Data collection The HAEDAT dataset is a collection of records of harmful algal bloom (HAB) events , maintained under the UNESCO
Intergovernmental Oceanographic Commission and with data archives since 1985.

Timing and spatial scale The satellite data were acquired from different seasons and across various phytoplankton bloom magnitudes between 2003 and
2020, and HAB data from 2003 to 2020 in the HAEDAT dataset were used.
March 2021

Data exclusions No data were excluded from analysis.

Reproducibility Our results could easily be reproduced with existing datasets.

Randomization Excluding data affected by clouds, a total of 0.76 million MODIS Aqua images from 2003 to 2020 were used to extract phytoplankton
blooms in global coastal area.

Blinding Not applicable in our study.

2
 Did the study involve field work? Yes No

nature portfolio | reporting summary

Reporting for specific materials, systems and methods
We require information from authors about some types of materials, experimental systems and methods used in many studies. Here, indicate whether each material,
system or method listed is relevant to your study. If you are not sure if a list item applies to your research, read the appropriate section before selecting a response.

Materials & experimental systems Methods

n/a Involved in the study n/a Involved in the study
Antibodies ChIP-seq
Eukaryotic cell lines Flow cytometry
Palaeontology and archaeology MRI-based neuroimaging
Animals and other organisms
Clinical data
Dual use research of concern

March 2021