Global Warming Impacts Micro-Phytoplankton at a Long-Term Pacific Ocean Coastal Station

Understanding impacts of global warming on phytoplankton–the foundation of marine ecosystems–is critical to predicting changes in future biodiversity, ocean productivity, and ultimately fisheries production. Using phytoplankton community abundance and environmental data that span ∼90 years (1931–2019) from a long-term Pacific Ocean coastal station off Sydney, Australia, we examined the response of the phytoplankton community to long-term ocean warming using the Community Temperature Index (CTI), an index of the preferred temperature of a community. With warming of ∼1.8°C at the site since 1931, we found a significant increase in the CTI from 1931–1932 to 2009–2019, suggesting that the relative proportion of warm-water to cold-water species has increased. The CTI also showed a clear seasonal cycle, with highest values at the end of austral summer (February/March) and lowest at the end of winter (August/September), a pattern well supported by other studies at this location. The shift in CTI was a consequence of the decline in the relative abundance of the cool-affinity (optimal temperature = 18.7°C), chain-forming diatom Asterionellopsis glacialis (40% in 1931–1932 to 13% in 2009 onward), and a substantial increase in the warm-affinity (21.5°C), also chain-forming diatom Leptocylindrus danicus (20% in 1931–1932 to 57% in 2009 onward). L. danicus reproduces rapidly, forms resting spores under nutrient depletion, and displays a wide thermal range. Species such as L. danicus may provide a glimpse of the functional traits necessary to be a “winner” under climate change.


INTRODUCTION
Effects of anthropogenic carbon emissions on our oceans include unprecedented warming, acidification, declining oxygen concentrations, and changes in nutrient cycling (IPCC, 2019). These physical and chemical changes are shifting the distribution, phenology, abundance, composition and trophic interactions of phytoplankton (IPCC, 2019). This is likely to increasingly have ecosystem-wide consequences, as phytoplankton undertake 50% of the world's photosynthesis Field et al., 1998) and underpin ocean productivity (Falkowski et al., 2004;Doney et al., 2006;Richardson and Poloczanska, 2008).
The physico-chemical environment regulates the structure of microbial assemblages in the ocean, and temperature is one of the most important environmental variables, directly and indirectly shaping phytoplankton community composition and abundance (Richardson and Poloczanska, 2008;Trainer et al., 2020). Ocean warming is predicted to reduce phytoplankton productivity in nutrient-limited tropical and mid-latitude waters (via reduced mixing and nutrient availability), but enhance it at higher latitudes (Richardson and Schoeman, 2004;Hays et al., 2005;Doney, 2006;Doney et al., 2012). Open oceans are also predicted to become increasingly stratified under global warming, resulting in a decrease in transport of nutrients into the upper layers. Coastal environments, on the other hand, where the thermal contrast between land and coastal oceans will increase under global warming, may see an increase in winddriven upwelling leading to an increase in nutrient availability and phytoplankton abundance (Falkowski and Oliver, 2007;Deppeler and Davidson, 2017).
Analyzing changes in phytoplankton assemblages in relation to ocean warming is difficult however, as phytoplankton exhibit orders-of-magnitude variability over seasonal, inter-annual and inter-decadal time scales (Zingone et al., 2010;Barrera-Alba et al., 2019;Benway et al., 2019). Nevertheless, sustained ocean time series measurements are now sufficiently long to investigate impacts of climate change Parmesan et al., 2011;Trainer et al., 2020). The Port Hacking 100 m station (hereafter PH 100m ), one of the Pacific Ocean's longest time series stations, is 5 km from Sydney on the east coast of Australia (Figure 1). As part of the National Reference Station Network in Australia's Integrated Marine Observing System (Lynch et al., 2014;Eriksen et al., 2019), PH 100m is in a highly complex and energetic oceanographic region of the Tasman Sea. It is dominated by the East Australian Current (EAC), a western boundary current of the South Pacific sub-tropical gyre. The EAC travels southward (up to 2 m s −1 ) along the continental shelf break until it separates and bifurcates at ∼32 • S into eastward and southward components. Over the past 60 years, the EAC circulation pathway in the bifurcation zone has changed in response to climate change (Ridgway and Hill, 2009;Wu et al., 2012;Sloyan and O'Kane, 2015), with long-term warming, increases in salinity and nitrate, and a decline in silicate over the past 60 years at Maria Island (Thompson et al., 2009). The changes to the path of the EAC has contributed to ocean warming off southeast Australia ∼3-4 times the global average (Ridgway, 2007;Ridgway and Hill, 2012).
Over the ∼90 years of phytoplankton research at PH 100m , there have been many short-term phytoplankton investigations undertaken (Ajani et al., 2016b) and references therein), but few studies over long time periods. Most recently, using phytoplankton community composition data, the realized niches of phytoplankton were examined over ∼50 years to assess if they were fixed or if they shifted in response to changing environmental conditions at this station (Ajani et al., 2018). Using the Maximum Entropy Modeling framework (MaxEnt) (Phillips and Dudik, 2008), a commonly used algorithm used to model species distributions, it was revealed that the phytoplankton community tracked changes in environmental conditions (temperature, salinity, mixed layer depth, nitrate, and silicate) over decadal time periods at this station. While this study provided the first evidence that biotic communities may be adapting to changing conditions at PH 100m , there has been no examination of the changes in relative composition across this same period in relation to regional warming.
A challenge of analyzing long time series is that often the sampling and analysis methods change over time, making it difficult to attribute change to a particular driver . For example, estimates of phytoplankton abundance differ depending on whether samples are collected using nets or bottles (Sournia, 1978). However, a careful choice of data analysis method can minimize these biases . For instance, by comparing changes in relative rather than absolute abundances of species, a more robust assessment of differences in sampling periods attributed to climate impacts can be made (Fodrie et al., 2010;Brown et al., 2011). One such index is the Community Temperature Index (CTI) (Devictor et al., 2008;Stuart-Smith et al., 2015). It is the overall temperature preference of species in a sample, weighted by their abundance, so that more abundant species have greater influence. As it is based on relative rather than absolute composition within a sample, it is less impacted by the variations in sampling methods over time. The CTI has been used extensively to assess the preferred temperature of communities of microbes, butterflies, fish and benthic invertebrates (Devictor et al., 2008;Bates et al., 2013;Cheung et al., 2013;Zografou et al., 2014;Stuart-Smith et al., 2015), although it has never been used for phytoplankton.
Here we use the CTI to investigate changes in relative phytoplankton community composition and temperature spanning 90 years  from PH 100m to examine the response of the phytoplankton community to long-term ocean warming. We answered two key questions: (1) Has the phytoplankton community changed overtime in response to warming; (2) If so, which species are responsible for this change? Insights from these findings are important for identifying sensitive indicators for further monitoring, developing a better understanding of the functional traits of "winners and "losers" under climate change, and the possible flow on effects to other trophic levels dependent on phytoplankton composition and production.

MATERIALS AND METHODS
From 1965 to 2019, there have been four major phytoplankton sampling surveys at PH 100m (34 • 7 3.36 S, 151 • 13 5.52 E, Figure 1). These surveys were carried out during 1965-1966, 1997-1998, 1998-2009, and 2009-2019 and had different sampling frequencies, methods and references ( Table 1 and Supplementary Figure S1). A fifth sampling campaign was included in the analyses due to its proximity to PH 100m . This earlier dataset by Dakin and Colefax (1933) was collected ∼6 km east of North Head, Port Jackson, Sydney (max. water depth 55 m; <40 km from PH 100m ) from 1930 to 1931. For the remaining analysis, this survey is considered part of the PH 100m dataset (Table 1 and Figure 1).
All phytoplankton samples were collected weekly, fortnightly or monthly using either discrete or depth integrated bottle samples (Dorn/Niskin) or phytoplankton mesh nets (20 µm or 37 µm) and preserved before microscopic examination. Cells were identified to the lowest possible taxon using light microscopy and where identification to species level was not possible, cells were assigned to a genus (e.g., Chaetoceros spp., Thalassiosira spp.) according to relevant taxonomic guides (Ajani et al., 2016b).
We used the HadISST1 dataset 1 to assess the temperature change in the area. This is a monthly, 1 • gridded product based on in situ and satellite temperature measurements (Rayner et al., 2003). We calculated an annual mean temperature for the 1 • grid square around the Port Hacking site.
To investigate whether warming at PH 100m changed the phytoplankton community over the past 90 years, we compared the temperature preference of the phytoplankton community through time (across sampling surveys) using the CTI (Devictor et al., 2008;Stuart-Smith et al., 2015). The CTI is a semiquantitative index that is reasonably robust to different sampling devices. It is calculated as a community-weighted mean of the thermal preference of a community based on the individual temperature preferences of each species [species temperature index (STI)], and as such is reasonably robust to different water volumes captured. It is calculated as a community weighted mean of the thermal preference of a community based on the individual 1 https://www.metoffice.gov.uk/hadobs/hadisst/data/download.html temperature preferences of each species (STI). It is normally calculated on a per sample basis. We calculated the CTI based on species that have been found across multiple surveys because this indicates species that are more likely to be consistently identified by different taxonomists, and it is less biased by species only identified in one particular survey. Thus, an increase in the CTI indicates that the proportion of species with warmer preferences have increased and/or there has been a decrease in the proportion of species with cooler temperature preference.
For the CTI analysis, the first stage was to calculate the STI for the species in the sampling surveys at PH 100m . To do this we combined 329 phytoplankton samples collected across five sampling surveys from 1930 to 2019 and from these samples, we selected the species that had been consistently identified across the sampling periods, rejecting those that were not verifiable based on published taxonomic information. From this group of species, we then selected those that had been identified in the IMOS sampling period and at least one other survey at PH 100m . From these, 30 species were retained as consistently identified and were included in the analysis of CTI (Supplementary Tables S1, S2).
The STIs were calculated using Australia's Integrated Marine Observing Systems (IMOS) long-term plankton datasets from the National Reference Stations (∼750 samples) (Eriksen et al., 2019) and the Australian Continuous Plankton recorder survey (∼6000 samples) (Davies et al., 2016), The National Reference Station and Continuous Plankton Recorder surveys have analyzed phytoplankton in a consistent manner for 12 years and provide abundance data for phytoplankton across the full temperature range of Australian waters and beyond. The sea surface temperature for each sample over the past 12 years was sourced from the GHRSST L4 gridded product, which gives daily temperature over a 10 km resolution based on the time, latitude and longitude of the sample 2 . We extracted every abundance record for the 30 retained species from this dataset and fitted these, with the associated temperature for each sample, to a Kernel Density Estimation (KDE) function. This is a nonparametric smoothing technique used to estimate the probability density function of a random variable. The temperature at the maximum point of the Kernel Density curve, the temperature of maximum abundance, is defined as the STI for each species.
Using the STI information for each species (s), the CTI was then calculated as the mean temperature preference of the species in a sample (Devictor et al., 2008): To understand if CTI had changed over time we used a linear model with the CTI as the response and Period and Month as predictors (Supplementary Table S3). Period was modeled as a factor including a level for each survey. We adjusted for the seasonal cycle by including Month as a factor in the model. This compensates for the bias introduced by the frequency and timing of the samples. If more samples in a particular sampling period were collected in the warmer months we would expect the CTI to be higher than if more samples were collected in the colder months. Accounting for the month compensates for this variation so the model is testing the long term response of CTI. Although CTI is robust to differences between sampling methods we also tested for effects of the variation in sampling design over time. We added predictors for Method (with two levels: bottle, net) and Depth (with two levels: 0-50 m, 0-100 m) (see Table 1) to the linear model, but neither predictor turned out to be significant (p > 0.05) supporting the assumptions we have made in using this analysis technique. Both predictors were dropped from the final model.
The CTI represents the overall temperature preference of the community present, and has units of • C. Because the CTI is defined on temperature, we expect it to be sensitive to changes in temperature (warming) and not other environmental variables. However, we investigated the importance of other available environmental variables (i.e., nitrate and phosphate; silicate was not included because it was only available for the most recent two sampling periods) in our linear models and found they were unrelated to CTI (Supplementary Table S3). Note that our CTI analysis and the focus on temperature contrasts with an analysis of the general drivers of phytoplankton communities using multivariate analyses of raw species abundances and including multiple potential environmental drivers. Unfortunately, such an analysis was not possible here because the phytoplankton data over the 90 years was sampled with different methods and so the CTI, with its ability to use relative abundance data, was considered more appropriate.
All the analyses were performed in R 3 , using tidyverse, ggplot2, and visreg packages. The KDE for calculating STI was performed using the density() function in the stats package.
Finally, as a decline in cooler-water species and an increase in warmer-water species can both contribute to an increase in the CTI, we then investigated the species primarily responsible for the increase in CTI over time. We selected the five most abundant species from each survey used in the CTI calculation. As the dominant species were mostly consistent across surveys, this resulted in a total of eight species, which we investigated further.

Caveats
There are several caveats of combining disparate datasets over time such as we have done in the present study. Firstly, the data in the present study focuses on the most frequently observed, 3 https://www.r-project.org/ taxonomically verified species, i.e., the diatoms. Dinoflagellates are often poorly resolved using preservatives such as Lugol's Iodine and light microscopy. On the other hand, however, the phytoplankton taxa used to calculate the thermal preferences in our current study (predominantly the diatoms) might be a true reflection of a diatom-dominated community at this location (Ajani, 2014;Ajani et al., 2014b). We have also carefully curated the species lists included in these analyses to include only species that an experienced analyst could reliably and repeatedly identify over the various surveys. We assessed taxa stability (defining characteristics, the ease of discrimination with light microscopy, any drawings, diagrams, photographs or taxonomic descriptions that would have been available to all analysts) over the entire period and made expert consensus decisions on which taxa to retain. We traced nomenclatural changes through the World Register of Marine Species (WoRMS) and thus could reconstruct the community composition recorded in the past with contemporary nomenclature. This is important as it not possible to re-examine the material from previous surveys, and indeed that may introduce a further bias in the results. We discarded observations only made to genus or higher levels as they do not add to or improve the analysis.

RESULTS
Species temperature index across all taxa ranged from 16.65 • C (Noctiluca scintillans) to 24.73 • C (Tripos candelabrum) with an average of 20.3 • C. Kernel density plots of these eight dominant species (seven diatoms and one dinoflagellate) show the range of their thermal preferences (Figure 2). The STI for the most abundant species ranged from 17 to 23 • C, clustered around 20.8 • C, similar to the mean water temperature over the IMOS sampling period (21 • C, Figure 3). Species with narrower, taller peaks are more stenothermal (e.g., Asterionellopsis glacialis), whilst species with broader curves are more eurythermal (e.g., Proboscia alata).
The model output indicates an increase in CTI over time, with the significance values of each period relative to the 1930s increasing with time (1965-1966, p = 0.4; 1997-2009, p    with warming over the same period at PH 100m (Figure 3). However, the warming of 1.8 • C was substantially greater than the increase in CTI of 0.8 • C for the same period. A clear seasonal cycle in CTI was also observed, with highest values toward the end of the austral summer (February/March) and lowest at the end of winter (August/September) ( Figure 4B). The seasonal variation in CTI agrees with the trend of water temperatures across seasons at PH 100 , although it is not as pronounced, and reflects the changing thermal preference of the community across season. There were substantial changes in the relative proportion of several of the eight dominant species for each time period (Figure 5). The generally cooler water, chain-forming diatom A. glacialis (STI = 18.7 • C) substantially decreased in abundance overtime (40% in 1931-1932 to 13% in 2009 onward). There were also declines over time in the dinoflagellate Tripos furca (STI = 18.2 • C, 9-1.6%), and the diatom Guinardia striata (STI = 22.3 • C, 7.3-2.3%), while the diatom Neocalyptrella robusta (STI = 20.0 • C) was observed in the 1930s (12%) but infrequently since (close to 0%). By contrast, the warmer-water diatom species Leptocylindrus danicus (STI = 21.5 • C; 20% in 1931-1932 to 57% in 2009 onward), Cylindrotheca closterium (STI = 21.9 • C, 0-7%) and P. alata (STI = 22.7 • C, 1.8-8.2%) all substantially increased in abundance over time. The cool-affinity diatom Guinardia delicatula (STI = 19.8 • C, 1-7%) also increased in abundance overtime.

Increase in Community Temperature Index at Port Hacking
Despite the dynamic nature of PH 100m , with a phytoplankton community that exhibits orders of magnitude variation over daily, weekly, seasonal and inter-annual time scales, there is still a clear increase in CTI at this location. The most pronounced of this "short term" variability is a result of winddriven or current-driven upwelling which brings cold, nutrientrich (phosphate, silicate, and nitrate) water into the euphotic zone and causes diatom blooms . These slope water intrusions can last between 2 and 22 days, most commonly from September to February (Lee et al., , 2007Pritchard et al., 2003). Seasonal variation is also well-defined at PH 100m , with the water column well mixed between May and September (decreasing in recent years) and highly stratified from December to March (Grant and Kerr, 1970;Ajani et al., 2001). Notwithstanding these well-documented patterns in the phytoplankton community at this location, these studies are still considered "short term" in relation to climate change. Using disparate datasets over much longer sampling periods (including the present study), emerging evidence suggests that the phytoplankton community at this location may track the magnitude and direction of environmental change also occurring in this region (Ajani et al., 2018). Indeed, the overall shift to a warmer community at PH 100m is most likely a consequence of the strengthening of the East Australian Current along this coastline. The increasing southerly flow of the EAC, and significant warming documented in this region, have not only seen the poleward shift in certain phytoplankton species (Ajani et al., 2014a,b;Buchanan et al., 2014) and larger species such as invertebrates and fish (Ling et al., 2009;Figueira and Booth, 2010;Last et al., 2011), but is predicted to cause a substantial shift in the microbial assemblage structure and function along this coastline (Messer et al., 2020).

Smaller Increase in Community Temperature Index Than Warming
Interestingly, the shift we observed in the CTI (0.8 • C) is less than that observed in the regional temperature data (1.8 • C). This may be for a number of reasons. While the increase in CTI at PH 100m may be due to an increase in the abundance of warmwater species such as L. danicus and a decrease in cold-water species such as A. glacialis, it may also be affected (reduced) by a decrease in other warm-water species (e.g., G. striata) or an increase in cold-water species (e.g., G. delicatula). These latter two species however, were far less abundant than L. danicus or A. glacialis at PH 100m , and although contributing to the overall CTI, would have had significantly less impact than the more abundant species.
Furthermore, it may be a consequence of the shorter time scales in variation of SST data compared to a lagged response by the phytoplankton community. The phytoplankton data examined does not include the whole phytoplankton community, but rather a subset of consistently and unequivocally identified taxa. STIs cannot be calculated for rare species due to the need for a minimum number of observations for a statistically robust result. As our final data set only includes relatively abundant, larger species that can be identified using microscopy and that we were confident were consistently and correctly identified across surveys (∼20% of the phytoplankton community), our results should be interpreted with caution; changes in rarer, smaller (picoeukaryotes) or species with narrow temperature preferences may also significantly impact the marine food web under climate change.
Species diversity may also confound this signal with changes in individual species having less effect on the CTI at PH 100m due to it relatively high diverse community composition. Moreover, our study has only considered temperature, and to a lesser degree nitrate and phosphate (found to be unrelated to CTI), yet it is the combination and interaction of all physico-chemical parameters (including nutrients, light and CO 2 ) that will influence marine community structure. The Port Hacking coastal location has seen an increase in temperature, salinity, and nitrate and a decline in silicate recorded at this station over the past 60 years (Thompson et al., 2009). However, not all responses to physicochemical variables are the same, and although out of the scope of the present study, larger scale oceanic processes or extreme event such as heat waves, must also been taken into account (Trainer et al., 2020). Finally, microbial assemblages are highly diverse and complex, often with thousands of species networking and interacting on both short (e.g., predation) and long-term time scales (e.g., competition), and with an enormous range of behaviors (sinking, buoyancy and vertical migration), life forms (single cell, colonies, spines etc.) and reproductive strategies (sexual, asexual, cyst formation etc.) (Barton et al., 2016). Therefore, although the CTI is a valuable community level metric to evaluate how well communities are suited to their thermal environments (Flanagan et al., 2019), other eco-physiological traits and bioprocesses, life history strategies and phytoplankton assemblage characteristics need to be considered when discussing the response of marine communities to temperature changes.

Potential Winners and Losers Under Climate Change
The species that contributed most to the increase in CTI at PH 100m was the warm-water species L. danicus. This diatom is a major component of coastal upwelling communities' worldwide (Nanjappa et al., 2013;Karthik et al., 2017) and a dominant member of the phytoplankton community along the south east coast of Australia (Ajani, 2014;Ajani et al., 2014bAjani et al., , 2016a. L. danicus reproduces rapidly (both sexually and asexually) and forms highly silicified resting spores under nutrient depletion. In a recent study to quantify the trait variability within and among closely related Leptocylindrus species, thirty-five isolates were established from six locations spanning 2000 km of the south-eastern Australian coastline (Ajani et al. under review). Comparable to levels of genetic diversity found across similar seascapes, significant intraspecific trait variability was observed for 8 of 9 traits measured (growth rate, biovolume, C:N, silica deposition, silica incorporation rate, chl-a, and photosynthetic efficiency under dark adapted, growth irradiance, and high light adaptation), providing evidence that L. danicus may rapidly adapt to fluctuating conditions. As such, it is this rapid reproduction, resting spore formation, and extremely high morphological and physiological strain variability that could be the functional traits necessary to be a "winner" under climate change.
Conversely, the cool-affinity, chain-forming diatom A. glacialis, which had a narrower thermal range compared with many other species from PH 100m , substantially declined in abundance overtime at PH 100m . This species has been reported to form highly dense (9.32 × 10 7 cells L −1 ), monospecific blooms as a result of upwelling events (cool, nutrient rich waters into the euphotic zone), most notably in India (Mishra and Panigrahy, 1995;Sahu et al., 2016), and along the east coast of Australia from Ballina Beach to Coffs Harbor (200 km) (Ajani et al., 2011). The first records of Asterionellopsis glacialis (as Asterionella glacialis) were by Castracane (1886) who examined Antarctic samples from the HMS Challenger voyages and as the name suggests, this mostly neritic species is more abundant in cold to temperate coastal waters (Hasle and Syvertsen, 1997). It is an "ecologically important planktoner" (Kaczmarska et al., 2014) because of its inclusion in the top 21 species (including L. danicus) that contribute ∼90% of the global biomass produced by diatoms (Leblanc et al., 2012). Under ideal conditions, it blooms repeatedly over one growing season, often in nuisance proportions, and is tolerant of a wide range of temperature, irradiance, energy (estuaries, open ocean and surf-zones) and salinity conditions (Smayda, 1958;Abreu et al., 2003;Kaczmarska et al., 2014). Wood (1964a) noted both species as important at Port Hacking and specifically noted Asterionellopsis (described using the invalid name Asterionellopsis japonica throughout his works, following Cleve (Crosby and Wood, 1959) as an expected component of winter blooms in June and July. In contrast to L. danicus though, the narrower thermal range of A. glacialis combined with the sustained increase in surface water temperature over the study period at PH 100m , allows us to speculate that this taxa may further decline at this site under climate change.
Other diatom species in the "top eight" that increased in abundance over the study period were C. closterium, G. delicatula, and P. alata. The centric diatom G. delicatula (STI = 19.8 • C) was rarely observed in the early monitoring campaigns (DC, OSL, and Ajani), but is now one of the most frequently observed species in the IMOS monitoring (Supplementary Table S1). This species is known to form non-toxic high biomass blooms worldwide but was not noted as present at Port Hacking (Wood, 1964a) in his extensive surveys. Control of bloom dynamics of G. delicatula in the northern hemisphere is linked to infection by phagotrophic nanoflagellates (Drebes et al., 1996) and viruses (Arsenieff et al., 2019), but to date this has not been established for populations at Port Hacking. Increased frequency of occurrence may be linked to apparent ability to grow at lower Si:C ratios than other diatoms (Rousseau et al., 2002), however, incomplete nutrient records in our early monitoring campaigns mean that long-term changes in nutrient ratios cannot be easily tested. Rousseau et al. (2002) also noted increasingly earlier bloom timings linked to higher than average water temperatures in the North Sea, and the survival of larger overwintering populations and this may also explain the increase in abundance over time at Port Hacking.
The cosmopolitan benthic marine diatom C. closterium (STI = 21.9 • C) was noted as quantitatively unimportant in coastal NSW by Dakin and Colefax (1933), but the "most commonest diatom in the estuarine environment" by Crosby and Wood (1959). Its broad distribution globally and ecological success is likely enhanced by its ability to disperse easily from sediments into the water column, high genetic diversity and genotypes with a range of thermal niches (Stock et al., 2019). In their global comparison of strains of C. closterium, Stock et al. (2019) reported a clonal culture isolated from Port Hacking showed optimal growth above 20 • C, with wide variability in strain specific maximum growth rates.
P. alata (STI = 22.7 • C) is a common oceanic species with wide distribution globally (Dakin and Colefax, 1933), with one varietal form (var. indica) described as extremely stenothermal (Schöne, 1982), however, this observation is not supported by the broad temperature preference exhibited by this species at Port Hacking (Figure 2). Noted as one of the top 20 taxa in the OSL study (Grant and Kerr, 1970), P. alata was one of 15 common phytoplankton taxa identified to have shifted southward in its distribution along the east coast of Australia (Richardson et al., 2015) relative to its mean position prior to 1990. Much of the published data on growth rates and ecological niches for this species focus on polar regions, partly because of the occurrence of distinct winter and spring forms (Jordan et al., 1991), and the species prevalence in large diatom-dominated spring blooms. There is comparatively little information on its autecology and favored niches in temperate communities. Species in the "top eight" that decreased in abundance over the study period included G. striata, N. robusta, and T. furca. G. striata (STI = 22.3 • C) is rarely mentioned in early studies from Port Hacking, but is another species noted for its low silica requirements and ability to form blooms in silica-deplete conditions (Schapira et al., 2008). G. striata was one of 87 North Atlantic phytoplankton taxa modeled (as Rhizosolenia stolterfothii) for future changes in geographic distribution based on historical CPR data  and future ocean conditions (2051-2100) (Barton et al., 2016). Consistent longitudinal shifts reflected dynamical changes in multiple environmental drivers (temperature, nutrients, stratification/mixed layer depth, etc.) on top of expected latitudinal range shifts commonly attributed to temperature alone.
The only dinoflagellate present in the "top eight" was T. furca (STI = 18.2 • C) noted to be tolerant of a wide range of temperatures and observed in estuarine, neritic and oceanic waters, but absent from polar and sub-polar seas (Dakin and Colefax, 1933). This species can be up to 50% of the dinoflagellate community in nearshore/estuarine environments (Wood, 1964a), often present year-round but in notably lower concentrations in winter (Dakin and Colefax, 1933). T. furca decreased in abundance from 9 to 1.6% over the study period suggesting a preference for lower temperatures, consistent with its inclusion in the 15 species showing a southward movement off the East coast of Australia (Richardson et al., 2015). Finally, N. robusta (STI = 20.0 • C) was present in the earliest dataset (DC), but not detected in OSL, Ajani or Ajani2 datasets. The species was described as a "warm oligo-eurytherm" diatom and proposed as an indicator of warm water currents from tropical or subtropical areas (Baars, 1982). It has not been observed in "cold" waters (Hernandez-Becerril and Del Castillo, 1996).

CONCLUSION
Based on one of the longest time series of phytoplankton in the Southern Hemisphere (90 years), and using a communitylevel temperature index, we found a significant warming signature in the phytoplankton community overtime. The phytoplankton CTI also showed strong seasonality, with highest values observed toward the end of the austral summer (February/March) and lowest at the end of winter (August/September), a well-recognized temporal pattern for this location. Species whose abundance changed substantially over time, and resulted in the observed increase in the CTI, were Leptocylindrus danicus and Asterionellopsis glacialis among others, providing insights into the potential functional traits-broad/narrow thermal range, formation of resting spores, and high phenotypic/genotypic variability-that may determine the adaptive capacity or survivability of species under climate change.

DATA AVAILABILITY STATEMENT
Publicly available datasets were analyzed in this study. This data can be found here: https://portal.aodn.org.au/.

AUTHOR CONTRIBUTIONS
PA, CD, and AR were responsible for the conceptualization of this work. CD performed the analysis. All authors were responsible for the writing, reviewing, and editing of this manuscript.

ACKNOWLEDGMENTS
We acknowledge all contributors, past and present to the PH 100m time-series dataset especially Prof. Gustaaf Hallegraeff, Ms. Alex Coughlan, Dr. Tim Ingleton, and Dr. Andrew Allen. National Reference Station data are sourced from the Integrated Marine Observing System (IMOS). IMOS is a national collaborative research infrastructure, supported by the Australian Government. All data are available through the Australian Ocean Data Network (AODN) portal (https://portal.aodn.org.au/search). The data used for calculating the STIs can be sourced from the AusCPR and NRS data collections and the other datasets are included in the Australian Phytoplankton Database (Davies et al., 2016). The data collection titles at the AODN are as follows: IMOS-AusCPR: Phytoplankton Abundance; IMOS National Reference Station (NRS)-Phytoplankton Abundance and Biovolume; The Australian Phytoplankton Database (1844-ongoing)-Abundance and Biovolume.

SUPPLEMENTARY MATERIAL
The Supplementary Material for this article can be found online at: https://www.frontiersin.org/articles/10.3389/fmars. 2020.576011/full#supplementary-material Supplementary Figure 1 | Sampling frequency across each sampling survey from 1930 to 2019 at PH 100m , southeastern Australia.
Supplementary Table 1 | A list of the taxa, their Species Temperature Index (STI), and the number of times each was reported in a sample from the PH 100m coastal station from each sampling campaign, 1930-2019. LowerSST and UpperSST represent lower and upper quartile temperatures from Kernel Density Plots for each species. Note NA = no report for this species in this sampling period. Species with * indicate those that were the top eight most abundant across all surveys and included in Figure 5.
Supplementary Table 2 | (A) Species not used as consistent identification could not be verified across surveys. (B) Species not used as didn't occur in the IMOS survey plus one other survey. (C) Species not used as not enough data to calculate a STI, i.e., rare species not seen more than 30 times over the IMOS NRS and AusCPR surveys across Australia.
Supplementary Table 3 | Linear model outputs. Note: In the second model neither the NOX nor the seasonal cycle are significant, but we see an increase in relative significance between the periods. The nitrate data was from the long-term monitoring survey at Port Hacking but was not associated directly with the samples, therefore a climatology was used to estimate the nitrate at each sample time confounding the effect of season across the model. The first model includes a wider time range, maintains a good level of significance between periods and has a weakly significant seasonal cycle. This is the more robust model.