Phenology of Trichodesmium spp. blooms in the Great Barrier Reef lagoon, Australia, from the ESA-MERIS 10-year mission

Trichodesmium, a filamentous bloom-forming marine cyanobacterium, plays a key role in the biogeochemistry of oligotrophic ocean regions because of the ability to fix nitrogen. Naturally occurring in the Great Barrier Reef (GBR), the contribution of Trichodesmium to the nutrient budget may be of the same order as that entering the system via catchment runoff. However, the cyclicity of Trichodesmium in the GBR is poorly understood and sparsely documented because of the lack of sufficient observations. This study provides the first systematic analysis of Trichodesmium spatial and temporal occurrences in the GBR over the decade-long MERIS ocean color mission (2002–2012). Trichodesmium surface expressions were detected using the Maximum Chlorophyll Index (MCI) applied to MERIS satellite imagery of the GBR lagoonal waters. The MCI performed well (76%), albeit tested on a limited set of images (N = 25) coincident with field measurements. A north (Cape York) to south (Fitzroy) increase in the extent, frequency and timing of the surface expressions characterized the GBR, with surface expressions extending over several hundreds of kilometers. The two southernmost subregions Mackay and Fitzroy accounted for the most (70%) bloom events. The bloom timing of Trichodesmium varied from May in the north to November in the south, with wet season conditions less favorable to Trichodesmium aggregations. MODIS-Aqua Sea Surface Temperature (SST) datasets, wind speed and field measurements of nutrient concentrations were used in combination with MCI positive instances to assess the blooms’ driving factors. Low wind speed (<6 m.s-1) and SST > 24°C were associated with the largest surface aggregations. Generalized additive models (GAM) indicated an increase in bloom occurrences over the 10-year period with seasonal bloom patterns regionally distinct. Interannual variability in SST partially (14%) explained bloom occurrences, and other drivers, such as the subregion and the nutrient budget, likely regulate Trichodesmium surface aggregations in the GBR.

Introduction Cyanobacterium Trichodesmium sp. is a marine diazotroph [1] phytoplankter found in all tropical and subtropical oceans [2][3][4]. In recent years, Trichodesmium erythraeum, the species occurring in Australian waters, has also been reported in the Mediterranean Sea [5] in late summer-early autumn and in UK coastal waters during winter [6]. Trichodesmium is known for forming extensive algal blooms that can cover hundreds of kilometers of the ocean surface. Aggregations have sometimes been so large that the sailors of the HMS Endeavor mistakenly identified a new shoal when the vessel sailed through the Torres Strait (North of Cape York, Australia) in August 1770. These persistent surface expressions occur toward the end of the algal bloom phase [7] and can last several days until the bloom subsides. Composed of senescent cells that are darker than healthy cells due to chlorosis [8], these surface matts increase the water-leaving radiance signal in the red-near-infrared (NIR) spectral region (680-750 nm), also called "red edge" [9], allowing their detection from satellite sensors for the past 30 years [10][11][12].
Blooms of Trichodesmium colonies likely play key roles in the ecosystem because of their ability to fix atmospheric nitrogen [13,14], thereby contributing to new nitrogen inputs in oligotrophic waters [15][16][17][18]. Trichodesmium colonies also contribute to the phosphorus budget by the uptake of phosphorus for growth [19,20], in addition to providing substrate and shelter to various organisms ranging from bacteria to crustacean larvae [21]. In the oligotrophic waters of the Great Barrier Reef (GBR), with chlorophyll-a concentrations that range from 0.2 to 0.8 μg.L -1 (e.g., S1 Fig; [22]), nitrogen fixation by Trichodesmium blooms may be a major source of new nitrogen, particularly in the offshore parts of the shelf [19,23].
However, because of the lack of systematic studies of Trichodesmium bloom events in the GBR, the understanding of their temporal and spatial dynamics within the system remains limited [4]. Based on field observations or satellite datasets, only a few studies have described the climatology of Trichodesmium in this region. The first study (1933) provided the seasonal abundance of Trichodesmium off the Low Islands (North of Cairns; Fig 1) over one annual cycle [24]. A major burst (~180 trichomes.L -1 ) occurred in August and was followed by another, although lower (~145 trichomes.L -1 ) but more sustained, peak in November-December. Four decades later, taxonomic counts of Trichodesmium at three stations off Townsville visited in 1976 and 1977 also reported the highest number of cells in August, December and March [25] (Fig 1). A 1992 study reported that Trichodesmium surface aggregations occurred mostly between August and February, but the spatio-temporal distributions of the events were uncertain [14]. More recent satellite-based studies have reported on Trichodesmium blooms in the western Pacific [26,27], the north Pacific subtropical gyre [28] and in the Cairns subregion [29]. In the latter, the phenology of Trichodesmium events during the ESA-MERIS sensor's mission from 2002 to 2012 occurred predominantly in August.
Thus, the complete spatial and seasonal occurrence of Trichodesmium in the GBR as a whole remains unresolved until this study, because to our knowledge, no previous research on the climatologic distribution of Trichodesmium has systematically analyzed these events over a decade-long period. This study builds on our previous research for the Cairns subregion [29] and extends this analysis to the whole of the GBR spanning 15˚of latitude (9.5˚to 24.5˚S; Fig  1). Using 10 years of MERIS Reduced Resolution (RR) satellite data and the Maximum Chlorophyll Index (MCI) [30], we provide new insights into the spatial distribution and decadelong temporal dynamics of Trichodesmium blooms in the GBR region. By comparing timeseries of MCI-positive instances along the GBR lagoonal waters, we identify the spatial and seasonal patterns of Trichodesmium surface expressions across the GBR subregions. Complemented by datasets of remotely sensed sea surface temperature (SST), wind speed and in situ nutrients, we assess the primary environmental factors associated with Trichodesmium surface expressions within the GBR.

Materials and methods
No specific permissions were required for these locations and research activities. The Australian Institute of Marine Science provided the nutrient datasets. No endangered or protected species were involved in this work.

The study region: Seasonal and hydrological aspects
The GBR has an area of~350,000 km 2 , with 10% covered by coral reefs (>3,000 individual reefs) [31]. This study focuses on the GBR lagoon, the water body located between the shoreline and the coral reef matrix (Fig 1). The width of the lagoon varies from up to 250 km in the southern GBR to less than 25 km north of 18˚S, as the reef matrix is closer to shore there ( Fig  1).
The GBR has a monsoonal climate with most rainfall occurring during the wet season (summer, November to April, >1,500 mm/year), resulting in episodic large river runoffs that may decrease the lagoonal water salinity (<< 33 ppt) and increase turbidity (Secchi disk depth < 1 m; [32]) in the near-shore lagoon [31,33]. These conditions are observed close to large river catchments, such as the Fitzroy River or the Burdekin River, or in shallow, tidal embayments, e.g., Broad Sound (Fig 1). Rainfall patterns and intensity are also influenced by the El Niño and La Niña climate phases, with El Niño summers characterized by higher than average SSTs, low wind speeds and clear skies, whereas La Niña summers are characterized by lower than average SSTs, higher wind speeds, and more cloudy conditions [34]. During the dry season months (May-October, <300 mm/year), southeast trade winds dominate. Wind strength eases in November-December and changes direction, allowing the intrusion of clearer water masses from the Coral Sea into the central and southern sections of the GBR.
For this study, the GBR was divided into five subregions (Fig 1), geographically defined as per the Natural Resource Management (NRM) areas. Typically, NRMs are used as reporting regions in the GBR report cards and Marine Monitoring Program (MMP). The water bodies vary in width from north to south. Cape York (CPY; Fig 1) is the subregion closest to the reef matrix. Most of the land in this subregion is undeveloped, with extensive land areas dedicated to nature conservation [35]; thus, the water quality is considered relatively pristine [36]. The Normanby is the major river and flows into Princess Charlotte Bay. Most of the coastal catchments in CPY are subject to heavy tropical rainfall. Although composed of small rivers, the CPY river catchments have an annual discharge equal to one-fourth of the entire GBR [31]. The Cairns subregion (CRN; Fig 1) is characterized by wet tropical rainforest and cleared land with agriculture being the major land-use. The primary rivers influencing the CRN are the Mossman, the Daintree, the Johnstone and the Tully rivers. In the Burdekin subregion (BDK; Fig 1), agriculture (cattle grazing in particular) is the main (90%) land-use activity [35]. The BDK is one of the largest (130,000 km 2 ) catchments in the GBR, and its main river, the Burdekin, can generate extensive river plumes during the wet season [33,37]. The Mackay subregion (MCK; Fig 1) is the farthest from the reef matrix (the GBR lagoon is the widest there) and hosts most of the sugarcane industry. The Fitzroy subregion (FIT; Fig 1)  River is the largest river system and catchment (140,000 km 2 ) discharging into the GBR [38,39].

Nutrient and metadata dataset from in situ sampling
In situ sampling of physical and chemical variables was conducted in each of the five subregions over the period 2002-2013. Data are summarized in S1 Fig. All water samples were collected with Niskin bottles, salinity and temperature were determined with a CTD SeaBird 911, and chlorophyll-a concentrations were measured fluorometrically using a Turner Designs 10 AU fluorometer. Further description of the methods used for this analysis can be found in [40,41].

Satellite image processing
During its mission, the MERIS ocean color sensor onboard the Envisat satellite provided reduced resolution (RR; 1.2 Km) imagery acquisitions of the GBR every two days on average. MERIS full resolution (FR; 300 m) imagery was intermittently acquired over Australia. Because the aim of this study was to use the entire MERIS mission, only RR imagery was used for this analysis. A total of 4,681 MERIS RR scenes from the third reprocessing, covering parts of the GBR between April 29, 2002, and April 4, 2012 (at~10 a.m. local time; UTC+10), were downloaded from the Optical Data processor of the European Space Agency (ODESA) provided by ACRI-ST (www.odesa-info.eu/). Subsets bounded to the five geographic subregions of interest were created, resulting in a total of 7,234 MERIS subscenes (S2 Fig). Land mass, cloud and sun glint-contaminated pixels were masked using the MERIS Level 1 quality flags.
The MERIS Maximum Chlorophyll Index (MCI) [42] was used for the detection and mapping of surface algal bloom expressions. The MCI algorithm is primarily designed for the detection of algal blooms with very high Chl-a concentrations globally (>30 μg.L -1 ) [30, 42], but such Chl-a concentrations are very high in comparison to the concentration ranges found in the waters of our study area: the GBR lagoonal waters are typically oligotrophic (Chl < 1 μg. L -1 ) (e.g., [22,43], S1 Fig), with higher Chl-a concentrations measured sporadically during the wet season [40,44]. Previous studies show that MCI can be used to detect surface blooms of slick-forming algal species such as Sargassum [45][46][47] in the Gulf of Mexico and Atlantic Ocean, in addition to Trichodesmium [10,29,48,49]. Computed from three MERIS bands in the near-infrared, namely, 681, 709, and 753 nm [50,51], the MCI was defined as follows: where L λ are level 1 TOA radiances at wavelengths λ using MERIS bands 10 (753 nm), 9 (709 nm) and 8 (681 nm), and k is a cloud factor set at the value of 1.005 and used to correct the influence of thin clouds. The MCI algorithm was applied to each Level 1 scene using the VISAT BEAM software toolbox (now replaced by the Sentinel Application Platform (SNAP)) and the default cloud correction factor of 1.005 to reduce the effect of thin clouds [52]. Processing the images to Level 2 included an atmospheric correction step, which might either result in flagging Trichodesmium sp. pixels as erroneous or as a saturated product over brighter bloom areas because radiance values for these pixels would likely be outside the expected range of the MERIS atmospheric correction [42,53]. MCI radiances typically vary between -3 and~15 mW -2 .sr -1 .nm -1 , with positive values above an MCI background level indicating phytoplankton-laden pixels. The MCI product does not have refined flags; thus, quality control was insured by computing the MCI for pixels with Level 1 top of atmosphere (TOA) radiance at 865 nm as <15 mW -2 .sr -1 .nm -1 [30], which further filtered pixels contaminated by land, high sun glint, haze or thick clouds. Disadvantages of the MCI include its sensitivity to submerged reefs [48,54], because coral reefs contain zooxanthellae that may trigger a positive MCI response [42]. A mask was applied to separate the GBR lagoon from the Coral Sea, thereby covering the entire reef matrix and avoiding the inclusion of possible false positives from the coral reef signals (Fig 1). The total counts of positive MCI values for each scene, hereafter called MCI positive instances and noted as MCI PI , were defined as: where t is the MCI threshold value, b is the background MCI value and n is the number of valid MCI pixels (i.e., "valid" refers to pixels for which a meaningful value is obtained not flagged for cloud, glint or any other L1 quality flag). The value of b is always negative, and as a result, the term (t-b) is always positive. The threshold t and median ocean background MCI values b of +0.4 mW -2 .sr -1 .nm -1 and -0.4 mW -2 .sr -1 .nm -1 , respectively, were selected based on Gower et al. [48,55]. However, we found that the sun illumination effects and the sensor-sun geometry were highly variable across regions covering 14˚of latitude, and the Gower et al. background value of -0.4 mW -2 .sr -1 .nm -1 was deemed not representative for our study region. For this study, t was a constant set at 0 mW -2 .sr -1 .nm -1 , and b was computed for each scene by taking the mean MCI value of non-bloom-contaminated pixels in the area surrounding the blooms. The component (t-b) varied in accordance with the viewing or sun angle.
Remotely sensed datasets of sea surface temperatures (SST) were used to support our analysis and assess the relationship between SST and Trichodesmium occurrences in the GBR lagoon. SST datacubes were composed of the daily NASA-MODIS Aqua skin temperature datasets for the period spanning the MERIS acquisition dates starting from July 4, 2002 to April 7, 2012. A box delineating the boundaries of each subregion and the same masks used for the MCI PI data covering the land, the reef matrix and the Coral Sea were also applied to the SST. SST daily average values were extracted for each subregion for the lagoon waters only. Across the entire time-series for all regions, SST varied from 17.5˚C to 29.9˚C, with an overall median of 23.5˚C.

Algorithm validation
The performance of MCI PI at retrieving Trichodesmium bloom surface expressions was assessed using field observations collated from four sources: the GBR-wide marine monitoring program (MMP), the GBRMPA database, CSIRO field campaigns and those listed in [56]. These records provided only qualitative observations (i.e., absence, presence) of Trichodesmium bloom surface expressions. The extent of the surface aggregations was not available, because these records were not collected explicitly for satellite validation purposes. Only the location (latitude, longitude), datum and the qualitative observation (absence, presence) were indicated. Records labeled as "present" between April 2002 and June 2009 were used (N = 432). For the purpose of this validation exercise, the search for the listed Trichodesmium events was conducted following McKinna et al [56]: independently of the sea-truth observation time, MERIS images within ±2 days of the observation date were used. The location of the field observation was used as the center pixel of a 10 x 10 pixel box (equating to~100 km 2 ). The maximum MCI PI value within that pixel box was extracted and used to determine absence or presence of Trichodesmium aggregations. When several images showed positive results for a single field observation, only the image the closest in time and date to the observation was considered. A coral reef mask was applied to each selected scene.

Wind speed datasets
Daily average wind speed was estimated in each of the five subregions from the eReefs CSIRO Environmental Modelling Data Node (www.ereefs.org.au/). The GBR 4-km hydrodynamic model was used to estimate the daily eastward and northward wind components from which the wind speed and direction were computed. The components were modeled at 10 a.m. local time to match the overpass time of the MERIS satellite sensor. Only daily wind speed data from September 1, 2010, onward were available, which limited the dataset to a total of 581 days (i.e., September 1, 2010-April 4, 2012).

Data analysis and generalized additive model (GAM) selection
To analyze the seasonal patterns of Trichodesmium blooms in the GBR, MCI PI were extracted from all MERIS scenes for the five subregions (Fig 2). Only daily subscenes with a minimum of 10 MCI PI were counted as containing bloom events and used for the subsequent analysis, and the value of b (Eq 2) was computed on at least 50 valid pixels. Using the R package "wq" for water quality data analysis [57], daily MCI PI were aggregated into monthly means to construct a 121-month time-series for the five subregions, and decadal trend models were computed using semiparametric smoothing techniques. An overall average of 13 MERIS scenes per month was available for monthly aggregations for a specific subregion. MERIS scenes with missing pixels were highly seasonal and occurred mostly during the wet season, i.e. January-March, due to cloud cover. The southern regions of MCK (2%) and FIT (4%) were the least affected by missing data, whereas CRN (33%) was the most affected (S2 Fig).
Two generalized additive models (GAMs) were built to model the influence of predictors (SST, time of year, interannual variation, and region) on the occurrence of Trichodesmium blooms, based on our response variable (MCI PI ). Wind was not used as a predictor, because the wind dataset was limited to only two years of the study period. Models were based on daily data collected between July 4, 2002, and April 4, 2012. Days for which < 20% of a scene was visible (i.e., � 80% cloud cover) were excluded from the analysis. The resulting dataset was zero-inflated (45%) and was analyzed using the two-step hurdle approach: (1) explain presence/absence, and (2) when present, explain extent [58].
Firstly, to assess the probability of a Trichodesmium bloom occurring, a GAM was fitted to the binary data (bloom or no bloom, hereafter called Bernoulli) using the mgcv package in R and a binomial link function. A GAM approach was used because several key predictor variables were considered likely to have non-linear effects: SST, time of year (or seasonality, represented by Julian Date) and interannual effects (represented by year). To test for regional differences in seasonality, the factor "region" was included in the model as both a main effect and as a varying coefficient of Julian Date (similar to an interaction effect). Secondly, a beta GAM was used to assess the variables influencing the extent of Trichodesmium blooms, within the subset of data in which blooms had been observed to occur. The response variable "bloom extent" was represented as the number of bloom pixels divided by the number of valid pixels in each image. This approach controlled for differences in the number of valid pixels in each image (e.g., non-cloud pixels) and for differences in the size of each region. The resulting variable was a proportion that did not include any zeroes or ones and was therefore ideally suited to a beta distribution. The response variable was cube root transformed to increase the spread of data and improve model dispersion, while staying within the constraints of the beta distribution. The beta GAM was fitted using a theta of 30 and the identity link function. Smoothers were applied to both SST and time of year data for the Bernoulli and beta GAMs.

Results
A summary of the physical and chemical variables in the surface waters of each subregion for the period 2002-2013 is presented in S1 Fig. Distinct differences were found in averaged key environmental parameters collected in the five subregions, with the largest spatial variations for temperature, particulate organic carbon (POC), nitrogen (PN), phosphate (PP), dissolved inorganic nitrogen (DIN) and organic phosphate (DOP), total dissolved nitrogen (TDN) and phosphate (TDP). Most of the nine variables showed slight north to south gradients across the five subregions: temperature decreased from CPY to FIT by up to 3.5˚C associated with the change in latitude (tropical to subtropical), whereas Chl-a showed little change (s.d.~0.2 μg.L -1 ) across the subregions. Some of the nutrients (e.g., POC, DIN, and TDP), with the exception of PP, showed minor increases in concentrations across the subregions on a north to south gradient. DIN and PN showed a very minor increase from CPY to BDK, and notably, DOP showed a >80% increase from CPY to FIT. This latitudinal gradient in nitrogen and phosphorus is consistent with recent findings on GBR water quality [40,59], with greater reported loads of DIN and PN in the southern GBR.
McKinna et al. used a set of 13 NASA-MODIS high-resolution (250 m) images coincident with Trichodesmium events observed in the central GBR (16˚-23˚S) to validate a binary classification algorithm [56]. The events used for the satellite retrieval assessment in their study were recorded in January 2005, April and October 2007, July, September, October, and November 2008 and February and June 2009. The algorithm retrieval was robust (85% success), albeit tested on a limited number of scenes (N = 13). Our validation exercise used 25 field observations with matching MERIS images (Table 1). Observations were reported in all months except March and December, with more than half (52%) recorded between September and November. A total of 19 these records had positive MCI PI , equivalent to 76% success. A total of nine MERIS scenes retrieved 13 of the sightings listed in McKinna et al. [56]. The distribution of the site locations is shown in  The number of MERIS images varied across subregions, with the lowest number of scenes for CRN (N = 1,117) and FIT with the most (N = 1,788) ( Table 2). A strong seasonality in satellite data availability was apparent, with the number of valid observations sharply decreasing during the wet season months because of increased cloud cover (S2 Fig) [60]. The magnitude of this seasonal decline decreased from CPY in the wet tropics to the subtropical region of FIT in which the wet season cloud cover was less pronounced. The monthly time-series of MCI background value, b, for each of the five regions is shown in S3 Fig. The difference in sun elevation and satellite observation angle variability across the 14˚of latitude and the seasons resulted in a 25% difference in mean MCI background value between the northernmost (CPY; -0.8 mW -2 .sr -1 .nm -1 ) and the southernmost (FIT; -0.6 mW -2 .sr -1 .nm -1 ) regions. Highest MCI background values were observed around June (early dry season) for all regions.
The multi-annual times-series of monthly averaged MCI PI for the five subregions are shown in Fig 2. The two southernmost subregions MCK and FIT accounted for >60% of all MCI PI for the GBR (Fig 2). Identical patterns were observed when these time-series were presented as ratios between daily MCI PI and the corresponding number of satellite observations for each region (data not presented). A general increase in MCI PI occurred from north to south, with CPY and CRN with the lowest (N = 272) overall number of monthly maximum counts (Fig 2) in comparison with the southernmost regions MCK and FIT for which monthly counts reached up to >2,000. Years with the most MCI PI varied across subregions but all occurred from 2007 onward with the exception of BDK (Fig 2; Table 3). The decadal seasonal signals (Fig 3) and their phenology (Fig 2) explicitly demonstrated that the peak of surface expressions captured in the imagery had a north to south gradient. The difference in algal bloom sizes also had a north to south increasing gradient in algal bloom area coverage, which was aligned with the overall increase in nutrients from northern to southern regions (e.g., The occurrence and extent of Trichodesmium blooms were analyzed following the two-step hurdle approach to: (1) explain the occurrence of blooms and (2) where present, assess the variables influencing the extent of Trichodesmium blooms. The Bernoulli GAM modeling the probability of a Trichodesmium bloom occurrence was fitted to the Bernoulli model assumptions (normality of residuals and homoscedasticity) and displayed an R 2 of 14% (N = 3,732). All of the variables included in the model had significant effects on the probability of a Trichodesmium bloom occurring (S1 Table). The effect of SST was strongly significant and non-linear, with a steep increase in probability of a bloom with increasing temperature up to approximately 24˚C and a near-plateau in probability of a bloom at temperatures above 24˚C (not shown). From the GAM analysis, the probability of a bloom occurrence peaked near the 200 th day of the year (mid-July) in the northernmost regions and progressively later in the year at more southerly regions (S1 Table), corroborating the graphical analyses of the monthly patterns in MCI PI (Fig 3). The beta GAM modeling the extent of Trichodesmium blooms was fitted to the beta model assumptions [61]. The prediction model was not over-dispersed (<1) and displayed an R 2 of 18.5% (N = 2,106) (S2 Table). A non-linear effect of SST was observed, with a statistically nonsignificant peak in bloom extents at 23 to 24˚C (p = 0.319). The largest measured influence on the extent of Trichodesmium blooms was subregion, with a highly significant regional main effect, in addition to a changing pattern in seasonal effects by region (S2 Table). Year had a weakly positive effect, with increasing extent of Trichodesmium blooms over the 10 year study period observed in all regions (Fig 6).
Overall, surface expressions of Trichodesmium blooms occurred more frequently from approximately July-August at northerly latitudes to November-December at more southerly latitudes. In CPY, surface expressions occurred for the months of August-September, with a second peak in November (S4 Fig). These later bloom events were mostly located in or close to Princess Charlotte Bay (Fig 1). In CRN, surface expressions tended to appear between June and September (S5 Fig) [29]. In BDK, the surface bloom seasonal activity (S6 Fig) was similar to that in CRN, although blooms could last until December and tended to occur near the mouth of the Proserpine River. In contrast to CPY and CRN, the spatial distribution of Trichodesmium surface expressions in BDK was widespread across the entire subregion and showed a higher frequency in pixel counts (S5 Fig). MCK and FIT had the greatest frequency and extent of surface blooms on the GBR; peaks in these regions occurred from September to November (S7 Fig): surface expressions in these subregions were widespread both along the coastline and across the lagoon. FIT had the highest percentage of bloom coverage over the lagoon for all subregions considered, with bloom spatial extents that covered up to 80% of the region in October. MCK had the second largest coverage (>60%).
The average daily wind speed corresponding to events with MCI PI larger than 100 varied between 1.6 and >8 m.s -1 (N = 313) across all regions (Tables 4 and 5). Although no strong correlation was found between the two parameters, MCI PI increased when wind speeds dropped below 6 m.s -1 . The average wind speed was 1.66 m.s -1 lower during large bloom Table 3. Monthly and yearly time-series statistics: Maximum MCI PI computed for each subregion when monthly and yearly aggregations were considered; dates are indicated between brackets. Decadal MCI PI (last column) and associated percentage of the total (all regions considered) were based on monthly aggregations over the 10-year period. Spatio-temporal variability of Trichodesmium blooms in the Great Barrier Reef lagoon from MERIS imagery events (>2,000 MCI PI ; 5.27 m.s -1 ) than that during smaller (<500; 6.93 m.s -1 ) blooms (Tables  4 and 5). The wind direction recorded for these aggregations, independently of their size, was largely (>80%) from the east-southeast, which indicated that the floating aggregations were pushed toward the shore.

Discussion
On a global scale, the onset and termination phases of Trichodesmium blooms depend on many environmental factors [62,63], which include the combination of optimal salinity [64] and light [65,66] conditions, nutrient (iron and/or phosphorus) availability [67], the lack of predation (e.g., copepod grazing or viral lysis) [68,69], water column stability [70] and temperature [71,72]. In the Indian Ocean (off Zanzibar, 6˚S), regular Trichodesmium aggregations are likely influenced by the warming and shallowness of the mixed layer in the summer, and the blooms are spatially sustained by iron supply carried eastward from the upwelling region south of Madagascar [73]. In the Red Sea (20˚-25˚N), Trichodesmium peaked in July Spatio-temporal variability of Trichodesmium blooms in the Great Barrier Reef lagoon from MERIS imagery [74]. From our results, seasonal and spatial patterns of Trichodesmium aggregations were evident for the GBR. Previous studies identified the possible environmental triggers of these highly seasonal surface blooms as a combination of low wind speed, nutrient availability (required for a large bloom to occur and persist [75]) and a marginal role played by SST. Yet, it has also been shown that Trichodesmium abundance and seawater temperature may not have clear relationships [63,76,77], in contrast to other studies in which temperature plays a key role [15,62,72].  [2] used SeaWiFS satellite reflectances to map the spatio-temporal distributions of Trichodesmium events at a global scale over a six-year period (1998)(1999)(2000)(2001)(2002)(2003). The GBR was masked out because their study solely focused on open ocean regions (see Westberry and Siegel [2]). The influence of SOI on Trichodesmium blooms was also examined, and the authors reported a correlation of~25% with La Niña events corresponding to larger and more frequent events. The authors noted that this correlation was particularly evident at the start of their time-series (1998)(1999)(2000) because a strong El Niño-La Niña transition period occurred. A similar transition might be observed over the period 2009-2011 in our study region, which might partly explain the increased frequency of bloom events noted above. The 2011 La Niña year was characterized by a decrease in salinity, increased nutrients and Chl-a and warmed sea temperatures along the northeastern coast of Australia [78], which might have supported the Trichodesmium blooms observed in the northern and central GBR regions of CPY, CRN and BDK in those years. In the GBR, nitrogen (N) fixation by Trichodesmium and N inputs from river runoff are proposed as major sources of new N to the system [23, 79], but no study to date has measured the detailed temporal and spatial variability characterizing N fixation activity. From the GAM results, an increase in bloom frequency of Trichodesmium occurred over the study's 10-year period (2002)(2003)(2004)(2005)(2006)(2007)(2008)(2009)(2010)(2011)(2012), suggesting that their importance as a source of N has also increased over the same period.
In other systems, the extent and densities of Trichodesmium aggregates are controlled by nutrient availability, such as phosphorus [4]. Dissolved organic phosphorus (DOP) is the primary form of phosphorus in the GBR [80], and previous studies show that Trichodesmium can grow in enriched DOP environments, such as those found in the southern GBR (S1 Fig.) [81]. The concentrations of both Nitrogen and available Phosphorus are generally very low in the GBR (S1 Fig), with a median inorganic Nitrogen:Phosphorus ratio of < 2, which is considerably less than the canonical Redfield Nitrogen:Phosphorus ratio (16:1) measured in phytoplankton biomass [82]. This result suggests that the GBR is severely nitrogen limited and that the nitrogen-fixing Trichodesmium has a competitive advantage when compared to that of non-nitrogen-fixing phytoplankton.
The relative quantity of nutrients exported into the lagoon differs significantly between catchments [38] and therefore subregions (S1 Fig.), but the north-south difference in nutrient concentrations was far less pronounced than cross-shelf differences within a subregion [25]. The effect of nutrient fluxes on phytoplankton populations was evident in some of the regions, such as MCK [83], and whereas comparatively little nutrient export occurs during the dry season months, other physical forcing, such as wind and tidal mixing, can play a major role in the cycling of nutrients. The dry tropics catchments export greater quantities of dissolved nitrogen than do the wet tropics catchments, but in both systems, the concentrations are only sufficient for a few days of phytoplankton growth [84].
Dust storms are another important source of nutrient enrichment (Nitrogen, Phosphorus and iron in particular), shifting phytoplankton species composition from pico-cyanobacteria (Synechococcus, Prochlorococcus) [85,86] to larger-sized groups such as diatoms and Trichodesmium, which may have important implications for other levels of the trophic system. Several studies in other parts of the world highlight the controlling factor of dust storms over periodic Trichodesmium blooms and their correlation to colony abundance, acting as iron  Spatio-temporal variability of Trichodesmium blooms in the Great Barrier Reef lagoon from MERIS imagery fertilization that stimulates blooms [87,88]. Regular dust storms occur over the Australian east coast, and the most notable occurred on October 23-24, 2002, [86] and September 22-24, 2009, as large dust storms (500 km-wide dust plumes) that spread over 2,400 km and 3,500 km (covering the distance Sydney-Cape York), respectively. These storms resulted in phytoplankton blooms, including Trichodesmium, mostly in the southern GBR within a week of these events. Such occurrences are typical of the GBR, as it is an ocean region in which a rapid biological response to dust inputs is most likely to be observed [14,89]. To date, the knowledge about iron cycling in the GBR remains limited, and therefore, to ascertain its role in triggering Trichodesmium blooms is currently not possible.
Trichodesmium often forms blooms parallel to the direction of the wind and as drifters, tend to accumulate in convergence zones (e.g., [56]). Low wind speed is strongly associated with Trichodesmium aggregations because of the resulting water column stability [14]. In this study, we found an inverse relationship between wind speed and the size of the aggregations, leading us to hypothesize that a wind speed <6 m.s -1 is required to observe large (>2,000 MCI PI ) aggregations.
The GBR is the only large tide-dominated reef in the world, and the tidal energy exerts a strong influence on parts of the shelf and across the coral reef matrix [90], in particular near Broad Sound (MCK; Fig 1) and CPY in which large tidal ranges (up to 9 m) can occur [91]. Shallower water depths (25-50 m) are found in the northern GBR (from 18˚S northward), whereas the southern GBR is typically much deeper and can reach up to 80 m average depth in some areas, exceeding 100 m near the Capricorn Channel. Thus, in CPY, the increase in water column mixing could hamper the development of regular Trichodesmium blooms. In other regions, however, such as CRN and BDK, tidal flows do not significantly contribute to water transport [92]. Because wind-driven resuspension occurs mostly between May and October [93], a frontal separation may be created between the lagoon and the coastal zone, thereby trapping terrestrial material and nutrients in nearshore waters and favoring the appearance of Trichodesmium blooms. During wet season flood events, large volumes of freshwater runoff entering the GBR reduce light penetration and decrease salinity levels (along the coast in particular, <<33 ppt; [94,95]), and despite being a euryhaline genus, maximum growth for Trichodesmium tends to occur at salinities of 33-37 ppt. Thus, we hypothesize that the tropical wet season conditions are suboptimal for Trichodesmium growth in the GBR, in particular near large catchments such as FIT and BDK (Fig 1).
Several interacting environmental factors are likely involved in affecting the probability of the occurrence and extent of Trichodesmium blooms (e.g., surface mixing by wind or/and tides, nutrient inputs from river plumes or dust storms, and grazer density, among others) beyond the spatial, temporal, and thermal effects assessed in this study. However, despite the complex dynamics of Trichodesmium blooms, the combined effects of spatial, temporal, and SST factors alone may account for 14% of the probability of a bloom occurring and 18.5% of the bloom extent when a bloom occurs. Spatio-temporal variability of Trichodesmium blooms in the Great Barrier Reef lagoon from MERIS imagery

MCI PI limitations
MERIS MCI PI was an appropriate algorithm for the detection and monitoring of surface expressions in the study region. Compared with this study, the match-up analysis by [56] had higher detection success at 85% but used a much finer resolution of 250-m pixel size from MODIS-Aqua and therefore was more likely to detect small and filamentous patches. This MODIS-based analysis is in contrast to the present study's validation for which the MERIS RR 1.2-km pixel size was used (MERIS FR imagery was irregularly acquired over Australia and therefore was not useable as a consistent time-series). The MODIS algorithm by [56] also relied on an atmospheric correction option only available in an outdated SeaDAS version; therefore, it was not possible to extend that analysis on the entire MODIS time-series. The limited information, other than the location and datum, provided by the field observation dataset also likely affected the quality of the retrieval analysis. The validation exercise in the present study was used for a general assessment of this algorithm. We recommend that additional criteria such as the size of the aggregations and a confirmation of the genus (i.e., is it Trichodesmium sp.?) should be considered to complement a future Trichodesmium-specific field database. The MCI has been used in several previous studies for the detection, monitoring and analysis of the seasonal dynamics, spatial distribution, and coverage of bloom surface expressions in a large range of freshwater (e.g., [96,97]) and marine (e.g., [55]) environments. This index is versatile and with low penetration depth (<1 m) is fairly insensitive to bottom reflectance (e.g., [98]), which is a significant advantage in this relatively shallow coastal system in which Secchi depths can reach the seafloor of the GBR lagoon [32]. The MCI also has limitations. The index has a limited application for low (Chl-a<10 μg.L -1 ) phytoplankton biomass conditions, as observed in the GBR [99]. In this study, the MCI was not employed to detect phytoplankton blooms based on their Chl-a concentration but on their surface expressions for which the MCI is suitable. Our previous study in the Cairns subregion combined the use of the FLH, Chl-a and MCI MERIS products, which showed that MCI and Chl-a were mostly uncorrelated, thereby confirming that MCI was predominantly responding to surface bloom expressions in the GBR [31]. Another recently documented limitation is the sensitivity of the MCI to inorganic particles (suspended sediment) [100], but we are confident this factor did not affect our results because the surface bloom observations occurred well outside the runoff (wet) season, when large sediment plumes may be observed. The use of the MERIS MCI PI alone, however, is not sufficient to confirm the genus of the phytoplankton bloom. For example, in the southern GBR (FIT), few instances of the surface expressions identified in this work as Trichodesmium sp. surface expressions could in fact be Lyngbya majuscula, a benthic cyanobacterium, or Hincksia sordida, a filamentous brown algae, as both are known to form surface aggregations in Southeast Queensland coastal waters (MCK, FIT) [101,102]. Hence, as no other taxa are likely to lead to consistent and significant blooms in the GBR that would trigger MCI PI , we can confidently interpret the MCI PI patterns as Trichodesmium sp. surface expressions.
Finally, an inherent bias exists in ocean color satellite observations due to cloud cover that hampers the collection of observations, in particular during the wet season. This variability must be acknowledged when analyzing the frequency of observed clear sky, pixel-based blooms in tropical and subtropical regions. Although this study was primarily based on monthly composites, i.e., the aggregation of several scenes over the course of each month, the analysis of time-series composing >1,000 scenes for each subregion should compensate for this bias.

Concluding remarks
This study used MERIS MCI PI to assess the spatio-temporal dynamics of Trichodesmium surface expressions over the entire GBR lagoonal waters during the decade-long MERIS mission. To date, this study is the only one that evaluated both the spatial and seasonal distributions of Trichodesmium blooms over the entire GBR and for a decade-long period.
The results showed a north to south gradient in bloom sizes (increasing from CPY to FIT), in bloom timings, appearing later in the year, from CPY in July to FIT in October-November (Table 4), and in bloom frequencies, with the largest and most frequent surface blooms occurring in FIT, the southernmost subregion. A temperature of 24˚C and a wind speed <6 m.s -1 were associated with larger events. Nitrogen fixation by Trichodesmium remains an incompletely constrained component of the GBR nitrogen budget, particularly in the subsurface [20], and our results suggest an increased importance of Trichodesmium as a source of new nitrogen to the GBR.
Recent modeling efforts are attempting to predict the distribution and growth of Trichodesmium from physiological-and remote sensing-based models incorporated into larger hydrodynamic and biogeochemical models of the Great Barrier Reef lagoon through eReefs (https:// ereefs.org.au/ereefs) to inform environmental management of this region [103,104]. By analyzing the phenology of Trichodesmium surface bloom aggregations along the GBR over a decade from satellite imagery, the findings of our work directly benefit this ecosystem modeling effort. The recent launch of multispectral Sentinel-2 MSI (launched in June 2015) and Sentinel-3 OLCI (launched in February 2016) satellite sensors that have matching wavebands to those of the MERIS sensor will allow this research to continue.