Interannual stability of phytoplankton community composition in the North-East Atlantic

As primary producers, phytoplankton play a pivotal role in the marine environment and are central to many biogeochemical processes. Changes to phytoplankton community composition could have major consequences for wider ecosystem functioning and may occur in response to climate change. Here we describe multi-decadal variability in phytoplankton community structure using taxonomic data from the Continuous Plankton Recorder collected in the North-East At lan tic from 1969−2013, using a total of 42 diatom and dinoflagellate taxa. We considered a range of characteristics of community structure, including taxonomic diversity and community stability and disorder, and how these characteristics change in response to sea surface temperature, mixed layer depth and the North Atlantic Oscillation. We found that phytoplankton community composition was largely stable on interannual timescales. A change in community composition occurred between 1985 and 1995 due to an increased dominance of 2 diatom taxa (Rhizosolenia styliformis and Thalassiosira spp.); however, after this period, the community returned to its previous composition. Further, a community disorder analysis found that phytoplankton compositional structure became more rigid in recent years, which may lead to an eventual community shift in the future. In contrast to previous studies that revealed relationships between total phytoplankton abundance or biomass and environmental forcing, we found that community structure had, at most, a very weak relationship with the environmental parameters tested. Changes to the physical environment may therefore have less influence at interannual timescales on phytoplankton community structure than previously thought.


INTRODUCTION
The stability of community taxonomic composition largely depends on the capacity for individual taxa to withstand perturbation, with external forcing capable of altering a community over multiple spatio-temporal scales (Beaugrand & Ibanez 2002, Edwards & Richardson 2004, Beaugrand et al. 2008. Because phytoplankton have many ecological and biogeochemical functions, the alteration of taxa within a community in response to global climate change may have the potential to alter global biogeochemical cycles, ecosys-tem functioning and predator−prey relationships, causing trophic mismatch (Edwards & Richardson 2004, Beaugrand et al. 2008, Eggers et al. 2014.
In a community of taxa with similar resource needs, the community composition can be maintained through trade-offs between intrinsic population growth capacity and competitive impact through growth rates and nutrient affinity (Dakos et al. 2009, Doncaster et al. 2016. During periods of exogenous forcing, ecological theory predicts an initial rise of dominant competitors, which may take 'keystone' roles that impose structural rigidity on the community (Doncaster et al. 2016). Ongoing exogenous forcing eventually disadvantages slow-replicating taxa, generally including the dominant competitors, tipping the community into a more fluid structure with prevalence of fast-replicating 'weedy' taxa (Wang et al. 2019). Analysis of temporal variations in community composition may thus provide an early warning signal of exogenous forcing affecting the community (Beaugrand & Reid 2003, Edwards & Richardson 2004, Beaugrand et al. 2008, Doncaster et al. 2016. Due to the diverse role that phytoplankton play in the marine environment, much research has been de voted to fluctuations in the abundance, biomass and community composition of these organisms (Dutkie wicz et al. 2001, Beaugrand & Reid 2003, Sommer & Lengfellner 2008. Current work concerning longterm variability in phytoplankton populations has predominantly focussed on the abundance and biomass of the entire phytoplankton population or large functional groups within it, such as diatoms and dinoflagellates (Edwards 2001, Henson et al. 2009, Barton et al. 2015). For example, seasonal changes in North Atlantic phytoplankton have been well documented, and a significant fraction of the variability in functional groups has been attributed to fluctuations in the physical environment , Dandonneau et al. 2004, McQuatters-Gollop et al. 2007a). The physical environment has also been implicated as the dominant factor driving population fluctuations over interannual time scales (Beaugrand & Reid 2003, Leterme et al. 2005; however, this remains under debate, with some studies concluding that predator− prey and other ecological interactions drive the variability (Benincà et al. 2008, Dakos et al. 2009, Boeing 2016. The variability of phytoplankton on longer time scales (multi-decadal, geological) is less studied, especially in regard to shifts in community composition (Cermeño et al. 2008).
Sea surface temperature (SST) directly impacts phytoplankton growth and metabolic rates (Eppley 1972) and indirectly impacts phytoplankton abun-dance by altering stratification and thus light and nutrient availability (Edwards & Richardson 2004). Interannual variability in bulk phytoplankton biomass across the North Atlantic basin has been associated with increased SST, particularly in the sub-polar gyre region (Colebrook 1982, Beaugrand & Reid 2003, Beaugrand et al. 2008, Barton et al. 2015). Further, dinoflagellates have been shown to outcompete diatoms (on average) in warmer years (Boyd et al. 2010, Henson et al. 2012, Zhai et al. 2013. Temporal variability in the mixed layer depth (MLD) also influences the magnitude and composition of phytoplankton communities. In a model study, longer periods of stratification and shallower mixed layers increased diatom, coccolithophore and prasinophyte biomass compared with years of frequent deep mixing (Litchman et al. 2006). Observational studies have also linked shallower mixed layers to the dominance of smaller protist taxa over haptophytes (van de Poll et al. 2013, Brun et al. 2015.
Multi-decadal variability in North Atlantic phytoplankton abundances has also been linked to a basin-wide climatic oscillation, the North Atlantic Oscillation (NAO). The NAO index is a measure of the difference in sea level pressure between the Azores and Iceland, which results in shifts to atmospheric circulation causing broad changes to the physical ocean environment (Hurrell 1995). The phase of the NAO index is correlated with changes in the relative abundance of diatoms, which increase in the North-East Atlantic during positive NAO phases, typical of cooler conditions with deeper mixed layers (Henson et al. 2012). Conversely, dinoflagellates dominate in warmer conditions during negative NAO periods in the sub-polar region of the North At lan tic (Edwards 2001, Henson et al. 2012). Long-term positive trends in phytoplankton biomass have also been associated with long-term trends in the NAO index, suggesting that atmospheric variability may be a key driving factor in phytoplankton variability on interannual to decadal timescales (Barton et al. 2003, Henson et al. 2012, Zhai et al. 2013.
Investigations into multi-decadal-scale variability in North Atlantic phytoplankton using in situ observations have tended to focus on bulk phytoplankton properties such as total biomass or productivity (Edwards 2001, Beaugrand et al. 2008, while studies exploring community composition generally focus on functional groups (Litchman et al. 2006, Boyd et al. 2010, Barton et al. 2016. In this study, we analysed multi-decadal variability in the North Atlantic diatom and dinoflagellate community composition using data from the Continuous Plankton Recorder (CPR) survey from 1969−2013. The CPR survey provides a unique multi-decadal view of plankton across the North Atlantic basin that re solves individual taxa, and has followed a consistent sampling methodology since 1948 for phytoplankton and 1958 for zooplankton. We consider several ecological metrics describing the observed diatom and dinoflagellate community, which is hereafter is referred to as the 'phytoplankton community', and relate the temporal variability in these metrics to 3 environmental drivers: SST, MLD and the NAO index. The aim of this study was to evaluate the stability of phytoplankton composition in the North-East Atlantic, an environment that has changed substantially over the past half century, by testing the following hypotheses: (1) there was no significant, persistent change in the diversity and structure of the phytoplankton community in the North-East Atlantic between 1969 andand (2) there is no relationship between phytoplankton community structure and SST, MLD or NAO index in the North-East Atlantic.

Phytoplankton community composition
The area of study is defined by 4 CPR Standard Areas (B5, C5, D5 and E5, Fig. 1). This region en com passes sections of the temperate and sub-polar regions of the North Atlantic.
Information on phytoplankton taxa abundance was obtained from the CPR survey, part of the Marine Biological Association (the CPR survey was previously managed by the Sir Alister Hardy Foundation for Ocean Science), Plymouth. The CPR device is towed behind a 'ship of opportunity', usually at speeds of 15−20 knots at a depth of ~7 m (Richardson et al. 2006). Plankton were captured on a moving silk layer (mesh size 270 μm), preserved by the application of a second silk layer and stored in a tank containing 4% formaldehyde. Each CPR sample consists of a 15 cm × 10.5 cm section of silk, representing 10 nautical miles of tow. Due to the scale and spatial coverage of the CPR survey, phytoplankton cells are recorded as present or absent across 20 microscopic fields of view, from which abundances are calculated in accordance with the number of counts (Richardson et al. 2006). Due to the large mesh size of the CPR device, smaller phytoplankton taxa are likely under-sampled, although organisms as small as coccolithophores are regularly recorded (Rivero-Calle et al. 2015). Larger phytoplankton, such as armoured dinoflagellates and chain-forming diatoms, are better sampled (Edwards et al. 2013). Data from the CPR survey are therefore regarded as semiquantitative, and are best used at monthly or annual temporal resolution (Richardson et al. 2006). Here we focus on dinoflagellate and diatom taxa composition at the annual scale and refer to the sub-set of the total community that the CPR samples as the 'phytoplankton community'. For additional information on the CPR methodology and history, see Barnard et al. (2004), Reid et al. (2003) and Richardson et al. (2006).
A total of 30 281 samples used in the analysis were collected between 1969 and 2013 across 4 standard CPR areas of the North-East Atlantic (Fig. 1 . Data within the remaining years were linearly interpolated to create an annual mean abundance for each phytoplankton taxon (Colebrook 1975), unless 2 or more consecutive years were missing data. Phytoplankton abundances were then fourth-root transformed to reduce the influence of the most dominant taxa (Clarke et al. 2008), and all taxa with a frequency of <1% for all CPR samples were removed from the dataset. By removing taxa that contributed to <1% of the community, smaller and poorly sampled taxa such as coccolithophores were removed from the CPR data set. Our statistical analysis of community composition therefore focusses on diatoms and dinoflagellates, organised in a community matrix of up to 45 years by 42 taxa for each region. Phytoplankton colour index (PCI) data from the CPR survey were also included in our analysis. PCI is a visual assessment of the 'greenness' of the silk ('nil', 'very pale green', 'pale green' and 'green') that is subsequently assigned numerical values based on acetone extracts using spectrophotometric methods (Richardson et al. 2006). PCI data have been previously verified using satellite data (Raitsos et al. 2005(Raitsos et al. , 2014. The PCI is used as a semi-quantitative assessment of total community biomass, which accounts for diatoms and dinoflagellates, as well as other microand nano-flagellates that were destroyed or too damaged for taxa identification while nevertheless contributing to the retained pigments (Richardson et al. 2006).

Environmental data
Vertical profiles of ocean temperature and salinity were obtained from the Hadley Centre's EN4 qualitycontrolled dataset (at 1° resolution, downloaded from www.metoffice.gov.uk/hadobs/en4/) (Ingleby & Huddleston 2007, Good et al. 2013). Density was calcu-lated from monthly temperature and salinity values, and MLD was defined as the depth at which a density difference of > 0.125 kg m −3 from the surface was observed (Chiswell 2011, Chiswell et al. 2013. Annual averages were calculated for both SST and MLD values to correspond with the data from the CPR survey. The NAO index time series was obtained from the UCAR Climate Analysis Section (Hurrell et al. 2003, https: // climatedataguide.ucar.edu/climate-data/ hurrell-north-atlantic-oscillation-nao-index-pcbased). The NAO index used was the principal component based time series of surface pressure differences between the Azores and Iceland (Hurrell 1995). For the present study, we used an average of the NAO index for the months December, January, February and March of each year.
Nutrients were not used as an environmental factor in this study, because, although they were considered in the design, data were not available over the appropriate time and space scales to compare with the CPR data for all regions. As such, only physical parameters of the environment were retained for analysis.

Statistical analysis
Different statistical analyses were used to explore changes in community composition, relationships to environmental forcing and the stability and potential for change within the phytoplankton community over time. Analysis was undertaken using R version 3.3.2 (R Core Team 2016) and MatLab version R2019b (The MathWorks).

Characterizing long-term changes
To characterize temporal patterns and determine the presence of temporal changes in the SST and PCI time series, a combination of trend detection and change-point analysis was used. We used an approach that compares 8 models that are commonly used to represent environmental time series: (1) a constant mean (no change), (2) a constant mean with background autocorrelation, (3) a long-term trend, (4) a long-term trend with background autocorrelation, (5) change points in the mean, (6) change points in the mean with background autocorrelation, (7) change points in the slope of the trend and (8) change points in the slope of the trend with background autocorrelation ). The models with change points were fitted using the pruned exact linear time algorithm developed by Killick et al. (2012). To select the best model to characterize a time series, we used Akaike's information criterion (AIC). The smallest AIC value indicates the best model fit and the AIC difference between models is a measure of how significant the differences are . The R package 'EnvCpt'  automates the fit and selection between models. Results of EnvCpt analysis were used to identify change points and trends in the physical environment that may be linked to changes in the phytoplankton community.

Analysis of community composition
To evaluate changes in community composition over time, Bray-Curtis dissimilarity coefficients (C n ) were calculated for each of the 4 CPR regions of the North Atlantic (Bray & Curtis 1957, Anderson et al. 2011, Legendre & Legendre 2012. This technique calculates a diversity measure that accounts for both the presence and abundance of taxa: where n a is the total number of individuals in year a, n b is the total number of individuals in year b, and 2w is the sum of the lowest value of abundances for taxa found in both year a and year b (Bray & Curtis 1957). Bray-Curtis analysis matches abundances at a taxonby-taxon level, comparing each year to all other years in pairwise testing. Resulting values of C n range between 0 (2 identical communities) and 1 (2 entirely different communities) and take account of both the presence and abundance of taxa within a sample.

Relationship between phytoplankton composition and the physical environment
For each environmental variable, a dissimilarity matrix was produced in order to apply a Mantel test (Mantel 1967) that evaluates the relationship between the phytoplankton community structure and the environment. The Mantel test is a pairwise correlation to test the null hypothesis that there is no relationship between 2 dissimilarity matrices (i.e. community data and an environmental parameter): (2) where C n1 and C n2 are 2 dissimilarity matrices with N samples, in which a and b are row and column in-dices, the first being the community dissimilarity matrix and the second an environmental parameter. The resulting correlation coefficient, Z M , ranges be tween −1 (a strong negative correlation) and 1 (a strong positive correlation), with 0 as no correlation and a pvalue indicating the significance of the relationship.

Analysis of taxon dominance
To further examine the structure of the phytoplankton community, the dominance of each taxon was calculated as the percentage contribution of an individual taxon to the total abundance. Taxon dominance allows for examination of interannual changes in specific taxa, as opposed to the whole community quantified in Bray-Curtis analysis.

Analysis of community stability
To assess the stability and potential for change within the phytoplankton community, 2 metrics were applied to each of the 4 standard areas of the North-East Atlantic. The first metric is a measure of resilience defined as the rate of return of a community to its mean state after a disturbance (Pimm 1984). Communities with a high resilience to change demonstrate fast return times, whereas communities with lower resilience take longer to return to their mean state (Dakos et al. 2012, De Keersmaecker et al. 2014. In this study, the temporal relationship between the observed years is used as a measure of resilience, and is expressed through the first-order autocorrelation, defined as when consecutive errors are correlated. For the taxonomic time series, ts A , with N observations from time t 1 to t N , the resilience of a community is calculated as: ( 3) where ts A is the ecological time series (here the abundance of each diatom and dinoflagellate), and ts A (i) is the time series at time i (adapted from Dakos et al. 2012, De Keersmaecker et al. 2014). The resilience metric ranges between 0 and 1; the lower the value of ρ 1 , the more resilient a community is to disturbance and the faster it will return to its mean state after a disturbance.
The second test of resilience is a test of ecological variance, used to assess the stability of a community. Ecological variance is defined as the total variability of a community and provides a more general meas- ure of ecological stability than measures such as resilience (De Keersmaecker et al. 2014). The coefficient of variance (V ) is calculated as: ( 4) where ts A is the ecological time series and S is its standard deviation. Variance was normalised to the average of the time series (ts A 3) in order to aid comparison between sampled regions. Communities were deemed to have a relatively large sensitivity to change in the environment if they take a longer time to return to equilibrium after a disturbance and de monstrate a higher variance than less sensitive communities.

Analysis of community disorder
Community disorder (D 0 ) provides a method of evaluating changes in functional composition (Thé bault & Fontaine 2010, Doncaster et al. 2016. When community composition is disordered through time, the structure of the community is influenced by an unpredictable turnover of taxa. We analysed disorder following methods detailed by Doncaster et al. (2016), using the R function 'nestedtemp' in the 'vegan' library (Oksanen & Carvalho 2013, R Core Team 2016. Each disorder value was calculated from an m × n matrix of community incidence comprising m = 10 consecutive years with n total taxa present in 1 or more of the 10 years. The matrix was shuffled to order the years (in rows) by descending abundance from top to bottom and the taxa (in columns) by descending prevalence from left to right. A perfectly ordered matrix would have all presences in the top left, and all absences in the bottom right. The resulting divide between these 2 extremes is known as the 'threshold of perfect fill', and is calculated as: for matrix coordinates (x, y) with c selected to create a curve that integrates to the same area as the proportion of presences (Oksanen & Carvalho 2013).
Any observed presence to the right of the threshold, or absence to the left, is classified as an anomalous observation, and disorder is quantified as a function of the sum across all observations of these anomalies. Values of disorder range between 0 for a perfectly ordered community to 100 for a completely disordered community (Rodriguez-Girones & Santamaria 2006). In a perfectly ordered community (D 0 = 0), the year in the bottom row of the matrix has a taxa composition that is a perfectly ordered subset of the composition in the year above it in the matrix, and so on through all of the years. This means that community composition in any year is determined by an ordered accumulation and/ or attrition of taxa from previous years. In contrast, complete disorder (D 0 = 100) means that the taxa composition in one year provides no information about composition in the year for the rows above or below in the matrix. While disorder measures the unpredictability of community composition through time, biodiversity quantifies the number of species. Disorder and diversity values were calculated and correlated in order to inform on the relative frequency of occurrence of commonly abundant species. For each year of available data, disorder values were tested for correlation with both taxa richness, defined as the number of taxa present, and with Hill's N 2 biodiversity parameter, which considers the richness and evenness of the community (Hill 1973): (6) where A is the proportional abundance of taxon α in a community of n taxa, with abundance given by the square-root percentage presence (using the R package 'vegan', Oksa nen & Carvalho 2013). Both taxa richness (R) and biodiversity (N 2 ) were calculated at the last time point in the matrix used for each 10 yr period disorder calculation. Correlations were calculated using Pearson's correlation on paired first differences in biodiversity (or taxa richness), and disorder values were calculated for consecutive years over the time series (e.g. 1970−1969, 1971−1970, etc.). It must be noted that this analysis required a time series without gaps. Consequently, data from 45 consecutive years (1969−2013) were used in region B5, 32 yr (1969− 2000) in region C5, 20 yr (1969−1988) in region D5 and 45 yr (1969−2013) in region E5. Other years from C5, D5 and E5 that could not be interpolated were discarded.
The sequential correlation coefficients between richness (R) and disorder (D 0 ) provide information on the relative occurrence of different types of taxa. A prolonged period of negative coefficients indicates a dominance of more persistent taxa over ephemeral taxa, as upward fluctuations in biodiversity are associated with an ordered accumulation of taxa and downward fluctuations with disordered turnover. In communities of taxa that compete for similar resources, the most persistent taxa are likely to function as keystones, which shape community composi- tion through strong inter-specific interactions (Doncaster et al. 2016). A switch from negative coefficients to a period of positive coefficients indicates a change to an ordered sequence of losses and disordered gains, consistent with exogenous forcing and a prevalence of fast-replicating weedy taxa (Doncaster et al. 2016).

Temporal change in the physical environment and PCI
All physical parameters displayed interannual and decadal variability (Fig. 2). In particular, SST showed an increasing trend, with all regions displaying an eventual change to a warmer regime (Fig. 2a). However MLD remained between 35 and 65 m (Fig. 2b), and the NAO varied widely, spanning from its lowest values in 1969 and 2010 of −2.5 to its largest value in 1990 of + 2.5 (Fig. 2c). The EnvCpt analysis conducted on the SST time series indicated that region B5 had an overall increasing trend with a change point in 1978, whereas regions C5, D5 and E5 all showed a change point in the mean during 1994 that shifted the environment to a warmer regime.
The time series of the PCI anomalies demonstrated an overall increasing trend in all regions (Fig. 3). Ana lysis of the PCI time series using the 8 models tested via the EnvCpt framework revealed structural changes for each region throughout the time series. The best fit model for the PCI time series within regions B5 and C5 indicated a change in the mean of the time series, with change points in 1997 for region B5 and 1984 for region C5 (Fig. 3). Small AIC values for this model suggest a good fit in the model choice (Fig. 4). The region D5 time series was best modelled as mean with autocorrelation (model 2, constant mean with autocorrelation); however, the AIC value for this region was similar to model 6, which indicated a change point. The most robust model for the time series in region E5 was model 7, which indicated a trend with a change point in 1978 (Fig. 4).
Examination of the contributions of different taxa to the community throughout the study period re vealed few temporal changes, with the clear exceptions of the dominance of R. styli- Year 1 9 6 9 1 9 7 3 1 9 7 7 1 9 8 1 1 9 8 5 1 9 8 9 1 9 9 3 1 9 9 7 2 0 0 1 2 0 0 5 2 0 0 9 2 0 1 3 Ecological resilience (1 − ρ 1 ) and variance metrics (V ) applied to the biological data represent the stability of the phytoplankton community in the North-East Atlantic and its ability to return to community equilibrium after perturbation. Ecological resilience values (1 − ρ 1 ) were predominantly close to 0 for all 4 regions, indicating a capacity for communities to return quickly to equilibrium following a disturbance in composition ( Table 2). The phytoplankton community composition of region B5 was found to be the most stable, with a resilience value of only 0.0001 ( Table 2). The resilience value for region D5 was ex cep tionally high compared to all other regions (0.4618) and may be the slowest to return to equilibrium following a perturbation. The ecological variability experienced across the North-East Atlantic was very similar in all regions ( Table 2).
The most abrupt change in disorder (D 0 ) occurred in region B5 around 1995, with a greater than 2-fold increase in disorder at a time of rising biodiversity (N 2 ) (Fig. 6) and taxa richness (R, not shown here). Region E5 had 2 increases in disorder, the first in 1977 with return to low values of disorder before a second increase in 1997 (Fig. 6). The decline of disorder values after each dramatic increase demonstrates that the community consistently returns towards a stable state (Fig. 6). Such a return further supports the results of the ecological resilience test that demonstrated the ability of phytoplankton to return quickly to equilibrium following perturbation (Table 2).
Sequential correlations were calculated on the first differences of disorder with biodiversity and with species richness, to seek insight into the role of community composition in population dynamics (Fig. 7). Regions B5 and E5 showed a switch in the sign of sequential coefficients of correlation on first differences from 1997 onwards, to a more negative correlation between disorder and taxa richness for B5 (Fig. 7a), and to a more positive correlation between disorder and biodiversity for E5 (Fig. 7d). These switches show recent samples in B5 undergoing a sustained period of ordered accumulations and disordered losses of taxa, consistent with endogenous competition; in contrast, recent samples in E5 undergo a sustained period of more ordered losses of taxa and disordered accumulations, consistent with exogenous forcing. Samples C5 and D5 show little change in correlation coefficient, except in the most recent samples which calculate coefficients from the fewest replicates (Fig. 7b,c).

Interannual variability of phytoplankton community structure
This study provides new insight into the consequences of changes to the physical environment for phytoplankton community structure in the North-East Atlantic at a taxon level. Bray-Curtis analysis demonstrated that the community remained relatively stable between 1969 and 2013, regardless of the variability observed in the physical environment (Figs. 3 & 4). In region B5 during 1988, and region C5 in 1993, the community composition temporarily veered away from its equilibrium due to a substantial increase in the abundance of 2 centric diatom taxa, Rhizosolenia styliformis in region B5 and Thalassiosira spp. in region C5 (Fig. 5). However, we found no links with the physical environment to explain these years of dramatic community change.
The lack of long-term trends in taxa composition de spite the observed increase in SST could potentially be explained by the fact that the majority of the taxa in this region have large thermal ranges, many of which span the annual and interannual averages of SST for each region, as well as the large seasonal variability found in this region (Irwin et al. 2012, Barton et al. 2013. The community could thus be well adapted to withstand the changes in annual mean SST observed in the North Atlantic, leading to a weak coupling between community structure and the phys ical environment over interannual time scales. The dominant community members, such as R. styli formis and Thalassiosira spp., have cosmopolitan distributions spanning large areas of the North Atlantic and thus are expected to have a large thermal range, which is further evidence that broad temperature tolerance is important for phytoplankton community stability ( Barnard et al. 2004, Barton et al. 2016. For example, the dinoflagellate composition proved to be highly stable in terms of taxa dominance, with Ceratium fusus dominant amongst our observed taxa throughout the time series for areas B5, C5 and D5. C. fusus has a wide biogeographic distribution ( Barnard et al. 2004, Beaugrand 2004, large thermal range (Irwin et al. 2012, Barton et al. 2013) and is mixotrophic, which potentially allows this taxon to remain dominant under changing SST conditions through inter-specific competition with other taxa ( Barnard et al. 2004). Similarly, the observed diatoms were dominated by the genus Thalassiosira, which also has a broad biogeographical range ( Barnard et al. 2004, Beaugrand 2004, Richardson et al. 2006. The disorder analysis provided a unique assess ment of phytoplankton community composition in terms of structural rigidity and predictability over the study period. Changes in disorder demonstrated fluctuations in the relative dominance of persistent taxa over ephemeral taxa in the community, with the most abrupt change occurring in region B5 around 1995, with an increase in disorder at a time of rising biodiversity (Fig. 6). Supporting prior dominance analysis, C. fusus and Thalassiosira spp. may be 'keystone' taxa to the community structure, i.e. taxa that dominate the community regardless of the physical environment by occupying wide niches. Sequential correlations suggested a high turnover of ephe meral taxa in more recent years, in favour of more persistent and potentially competitively dominant taxa. As no large or persistent shifts were de tected in the community over the studied time period, ecological theory suggests that the structural rigidity of the community may increase over time with an increasing prevalence of competitively dominant taxa, re ducing richness, 52 Fig. 6. Results of disorder analysis for regions B5, C5, D5 and E5. Each plot displays the time series of disorder (D 0 , red) and biodiversity (N 2 , blue). Each red line joins the disorder value to its corresponding biodiversity value (vertical part), and stretches back through the 10 time steps used for calculation of disorder (horizontal part). Note the different scales on each y-axis and may therefore render it more prone to collapse from external forcing (Doncaster et al. 2016).
The analysis of PCI across the North-East Atlantic demonstrates an increase in the mean that contrasts with the stability suggested by the Bray-Curtis community analysis (Figs. 3 & 5). Previous studies have linked changes in phytoplankton community biomass, through the use of the PCI from the CPR survey and satellite derived chlorophyll a measurements, to changes in the physical environment over inter annual timescales (Raitsos et al. 2005, McQuatters-Gollop et al. 2007a, Martinez et al. 2016. As a proxy for biomass, PCI allows insight into a greater proportion of the 'true' in situ phytoplankton community that was not fully taxonomically resolved in this study. The initial data manipulations that were conducted (as well as the mesh size of the CPR device and limitations of light microscopy) exclude many smaller and potentially less cosmopolitan taxa. Previous studies have re ported a general decline in PCI until the late 1980s and then an increase beyond the 1990s, which tend to agree with the findings of our study (Fig. 3) (Edwards 2001, Beaugrand & Reid 2003, Raitsos et al. 2014, Martinez et al. 2016. Future work should exa mine a larger proportion of the in situ community for taxa-based population changes, including taxa smaller than those captured by the CPR and taxa outside the diatom and dinoflagellate functional groups, and consider larger geographical regions.

Implications for the bottom-up control of community structure
The relationship between phytoplankton and the physical environment on interannual timescales has been widely debated (Dakos et al. 2009, Prowe et al. 2012, Barton et al. 2015; however, the weak correlations between the physical environment and community composition in this study were unexpected and largely contrasted with the outcomes of prior studies. Beaugrand & Reid (2003) suggested that an increase in global temperatures would lead to a marked change in the organisation of phytoplankton communities within the North Atlantic, and such links have been demonstrated in studies showing that dinoflagellates dominate over diatoms in warmer SST re gimes (Leterme et al. 2005, Henson et al. 2012, van de Poll et al. 2013. The influence of the NAO on the marine environment across the North Atlantic has also been well documented (Barton et al. 53 Year  Drink water et al. 2003, Henson et al. 2012. The ef fect on phytoplankton is suggested to be that diatoms dominate in positive NAO years when conditions create a cooler climate, with stronger wind stress and deeper mixed layers, whereas during negative NAO periods (where the reverse conditions occur), dinoflagellates dominate (Edwards 2001, Henson et al. 2009. Increased temperatures were suggested to also cause a reduction in the overall cell size of the phytoplankton community and reduce the ecological diversity in warmer regions (Edwards et al. 2013). However, our Mantel test re sults found that for this region of the North-East Atlantic between 1969 and 2013, the phytoplankton community remained relatively stable with little change in composition, despite an increase in SST. Furthermore, when changes did occur, such as in region B5 during 1986, the community soon returned to its original state (Fig. 5). Our analysis of the PCI time series revealed structural changes for each region throughout the studied time series, in agreement with the interannual analysis of Barton et al. (2015). In their study, the authors suggested that physical changes in the environment in the North-East Atlantic (surface wind speed, heat flux, SST, MLD, stratification and turbulent kinetic energy) were strongly related to phytoplankton on seasonal timescales but had little influence on the long-term (interannual) structure of phytoplankton communities.
As interannual fluctuations in community structure were not statistically linked to changes in the physical environment, we speculate that the variability of phyto plankton populations on interannual time scales may be driven by other factors, e.g. chaotic variability (Dakos et al. 2009, Barton et al. 2015. Dakos et al. (2009) utilised a model to describe variability of phytoplankton over seasonal and interannual time scales without the requirement for external forcing. Interannual variability within the phytoplankton population was found to be a naturally occurring property of multi-taxa communities, with chaotic fluctuations occurring frequently without environmental perturbation (Dakos et al. 2009). Population chaos has also been demonstrated in mesocosm experiments in which community fluctuations occur without external forcing (Benincà et al. 2008). The results of the Mantel test determined that there was no relationship between the phytoplankton community structure and changes to the physical environment. However, variability in the community may have resulted from the chaotic nature of a multi-taxa phytoplankton community, and this could be the driver of community fluctuations on interannual scales (Dakos et al. 2009).

Stability of phytoplankton in the North-East Atlantic
The apparent stability of the phytoplankton community in the North-East Atlantic identified by our analysis suggests that the community is well suited to the environmental perturbations occurring within the study period. When multiple taxa are well adapted to withstand environmental perturbations, we can expect a community to remain stable throughout time as demonstrated by the North Atlantic community exa mined here (Colebrook 1982, Leterme et al. 2005, Beaugrand et al. 2008. Taxa thermal ranges and biogeographical distributions reported in the literature show that the majority of taxa enumerated in the CPR survey have a large thermal tolerance that spans across the total SST variability from region B5 to E5, likely reducing their susceptibility to population change due to variable SST (Irwin et al. 2012, Barton et al. 2013. Coupled with this, the prevalence of other sources of interannual variability, including the NAO, may contribute to the broad distribution of taxa across the North Atlantic and the resilience of the community to environmental change. The robustness of the phytoplankton community to long-term change in the physical environment has also been demonstrated in the fossil record, assessed with Bray-Curtis analysis (Cermeño et al. 2010), where the biodiversity and community structure of marine diatoms were found to recover after significant environmental perturbations. Cermeño et al. (2010) suggested that the high dispersal of marine phytoplankton, resulting in large biogeographical distributions, ensured the survival and recovery of the community (Cermeño & Falkowski 2009, Cermeño et al. 2010. The observed stability and predictability of the community contrasts with the results from some model studies, which predict that continued global warming will alter the phytoplankton community composition (Barton et al. 2010, Bopp et al. 2013, Dutkiewicz et al. 2015. For example, Barton et al. (2016) demonstrated that phytoplankton communities appear poised for a biogeographical shift and community re-shuffle in the dominance of key members in the second half of the 21 st century. The shift was predicted to result in decreased taxa richness in sub-tropical regions and an increase in sub-polar regions, as only approximately 35% of local taxa survive in future environmental conditions (Barton et al. 2016). The Bray-Curtis analysis and stability metrics shown here demonstrate that even with changes to the physical environment, the community was able to return to its original state after perturbation, suggesting that larger temporal and spatial studies may be required to detect significant shifts in community composition linked to changes in the physical environment, and may become more evident in the future with the continued pressure of global climate change on the marine environment.
The findings here also contrast with results from the North Sea region, where a regime shift in the 1980s resulted in an increase in PCI values which was linked to increased SST (McQuatters-Gollop et al. 2007b, Raitsos et al. 2014. The North Sea is a relatively closed basin and is shallower than the North-East Atlantic, which may increase the impact of changes to the physical environment compared to a more open ocean basin. The regime shift described in the North Sea was also detected in PCI, which encompasses information from a greater proportion of the in situ community. Changes in the phytoplankton community may therefore be occurring in the smaller size fraction not captured in our study (McQuatters-Gollop et al. 2007b, Raitsos et al. 2014, and the relationship with the physical environment may be more complex; we therefore suggest that non-linear relationships should be examined in a future study. While a regime shift was not detected in the North-East Atlantic, the results of the disorder analysis highlighted the potential for future shifts in the community structure, as some areas of the present community appear to become increasingly rigid towards the end of the time series (Doncaster et al. 2016).

CONCLUSION
The findings presented here suggest that the mecha nisms driving long-term ecological change in phytoplankton community composition within the North-East Atlantic are highly complex and perhaps more uncertain than previously described. The temporal and spatial extent of the CPR survey allows for a unique insight into the long-term variability of phyto plankton communities in the North-East At lantic due to its unchanged methodology and consistent taxonomic detail. We found little or no correlation between diatom and dinoflagellate community composition and SST, MLD or the NAO index.
There was a change in the community composition be tween 1985 and 1995 from its initial state due to an in creased dominance of 2 diatom taxa; however, after this period the community returned to its previous composition. As such, the diversity of the community remained relatively unchanged over the time series despite a shift to increased SST in the mid-1990s. This stability was confirmed by tests of ecological resilience, which revealed high resilience of community composition over the time series. The re silience may be explained by the fact that the majority of taxa identified by the CPR have thermal tolerances greater than the environmental variability ex perienced. Disorder analysis, however, suggested a recent increase in structural rigidity in some communities, which may be an early warning signal to an eventual community collapse, and should be investigated in a future study. This study examined changes in the diatom and dinoflagellate fraction of the phyto plankton community captured by the CPR survey. We recognise that the community assessed here does not contain the smaller size fraction of the in situ phytoplankton community, which may manifest different patterns of variability than the larger size fraction captured by the CPR survey.
To conclude, our analysis demonstrates the ability of the phytoplankton community to return quickly to equilibrium composition after environmental perturbation, with a relatively stable community composition between 1969 and 2013 in the North-East Atlan tic. With no clear relationship with the physical environment, we suggest that interannual variability in phytoplankton community composition may be driven by other factors, including internal community forcing such as chaos (Dakos et al. 2009), and it is therefore important that non-linear relationships are considered in future studies.