Study of Temporal Variations in Species–Environment Association through an Innovative Multivariate Method: MixSTATICO

Abstract

The study of biotic and abiotic factors and their interrelationships is essential in the preservation of sustainable marine ecosystems and for understanding the impact that climate change can have on different species. For instance, phytoplankton are extremely vulnerable to environmental changes and thus studying the factors involved is important for the species’ conservation. This work examines the relationship between phytoplankton and environmental parameters of the eastern equatorial Pacific, known as one of the most biologically rich regions in the world. For this purpose, a new multivariate method called MixSTATICO has been developed, allowing mixed-type data structured in two different groups (environment and species) to be related and measured on a space–time scale. The results obtained show how seasons have an impact on species–environment relations, with the most significant association occurring in November and the weakest during the month of May (change of season). The species Lauderia borealis, Chaetoceros didymus and Gyrodinium sp. were not observed in the coastal profiles during the dry season at most stations, while during the rainy season, the species Dactyliosolen antarcticus, Proboscia alata and Skeletonema costatum were not detected. Using MixSTATICO, species vulnerable to specific geographical locations and environmental variations were identified, making it possible to establish biological indicators for this region.

1. Introduction

Water sustainability has often been treated as sustainable freshwater management for human consumption. However, from a more holistic point of view, water, both freshwater and saltwater, is an important resource and essential for sustaining ecosystems that generate or maintain the environmental conditions necessary for sustaining life in general [1]. The demand for water in rural agriculture, industrial operations, human consumption and other purposes is high, but in recent years, the water contamination caused by chemicals, pharmaceuticals and cleaning processes [2] have polluted rivers, creating serious problems in the marine ecosystems once this contamination reaches the sea.

For the United Nations (UN), water is an essential element in the adaptation to climate change [1], and is a significant problem that society and the environment are currently experiencing. Climate change has led to alterations in the natural variability of the environment, which is reflected worryingly in extreme temperatures and precipitation, affecting pollution in the hydrological component due to the high nutrient values generated [3]. Environmental changes are considered one of the main elements that cause alterations in the environment and ecology, leading to vulnerability [4]. To assess ecological vulnerability, it is necessary to have quantitative and qualitative knowledge about ecosystems and to thus understand the factors that cause their vulnerability [5].

In the study of the hydrological component, mathematical models [6] have increasingly been applied through predictive and descriptive methods that are classified into deterministic and stochastic. Stochastics are characterized by using multivariate statistical methods [7]. These methods are considered helpful for studying water quality [8] using evaluation factors such as correlation, covariation, similarities and distances [7]. In addition, in order to analyse temporal hydrological behaviour from a stochastic perspective, it is necessary to study the temporal–spatial correlation of hydrological parameters [7].

In this way, it is of vital importance to study the sustainability of water associated with climate change, not necessarily in its most common environment, but rather from a holistic perspective.

It is known that climate change affects the sustainability of water and its ecosystem, whether inland waters, seawater or the open ocean, as well as the environment and society at large.

Multivariate statistical methods have been developed to study data matrices with multiple variables in order to reduce the dimensionality of a matrix. These methods have a long history in science. The first was Principal Component Analysis (PCA) [9,10,11], a technique that projects the observation matrix onto a bidimensional space built from several inter-correlated orthogonal quantitative variables called principal components. The components are linear combinations of the original variables. The first component represents the largest variance retained by the method, the second component represents the largest variance possible while orthogonal to the first, and so on for the other components.

Subsequent developments include Factorial Analysis (FA) [12,13], which projects the correlations between variables and among variables and factors; Correspondence Analysis (CA) [14,15], which studies the association between the categories of two qualitative variables; Multiple Correspondence Analysis (MCA) [16], which extends CA to several qualitative variables; and Canonical Correspondence Analysis (CCA) [17]. The latter is an extension of CA widely used in biology which relates community composition to known variation in the environment (two matrices). The ordination axes are linear combinations of the environmental variables.

These methods are useful in analysing data arrays with a two-way dimensional structure, as they have the constraint that the number of variables must be less than the number of individuals (PCA). However, this condition cannot be met in some areas of research that analyse several variables to describe communities and environments, and for ecological data sets that are for the most part multivariate [18]. At present, with the advancement of multivariate data analysis techniques, it is also possible to analyse data that consider multiple conditions, which give rise to arrays that have a k-way structure and multi-way data arrays. For the analysis of this new array set, multi-way multivariate methods are used that allow multiple conditions to be studied in combination, identifying underlying patterns that other methods are unable to identify. There are valuable contributions to these multivariate methods cited below.

1.1. Contributions of the French School

Multiple Factor Analysis (FMA) [19,20,21] analyses information from different variable sets (matrices) defined on the same set of individuals (observations). This method performs PCAs in two phases: first it normalizes each matrix (dividing its elements by the first eigenvalue obtained from applying a PCA to each matrix) and then merges the matrices into one that represents the overall structure, then a PCA is applied to that global matrix.

STATIS [22] is a technique used to study a set of K-way matrices and is also known as ACT-STATIS [23]. It analyses data matrices by calculating Euclidean distances between configurations of the same individuals measured by different variables and K conditions (times or situations). This method bases its analysis on operator matrices (covariance matrices); it is a generalization of PCA and begins by analysing the similarities between the K-matrices. From these, it obtains weights (first eigenvector) used to find the stable part of the common structure (compromise or consensus) of the set of matrices, which is calculated as the weighted linear combination of the K-matrices and represents it in a low dimension space. It also allows the original individuals of each k-matrix (k = 1, …, K) to be projected over the consensus dimensions obtained.

STATIS developments have had important and numerous contributions and can be classified according to the input data (same variables or individuals in the k-matrices), according to the weights assigned to the compromised matrix, when considering external information and if pairs of matrices are available in different situations or times. These methods started with the STATIS and STATIS-DUAL [22], to be continued with X-STATIS or Triadic Partial Analysis (PTA) [24], STATIS-CoA [25], STATICO (STATIS and COINERCIA) [26,27], DOUBLE STATIS (DO-ACT) [28], K + 1 STATIS [29], DISTATIS [30], STATIS-4 [31], Kernel-STATIS [32], COVSTATIS and COSTATIS (COINERCIA and STATIS) [33], POWER-STATIS [34], Anisotropic STATIS or ANISOSTATIS [35], STATIS-LDA [36], INTERSTATIS [37], Sir-STATIS [38], HiDISTATIS and DiDISTATIS [39], CATATIS [40], STATICO-CoA [41] and recently CLUSTATIS [42].

However, despite all these contributions, none of these methods analyse matrices with mixed data.

The methodological basis of this study is based on the methods of the French school, but it is also necessary to present the methods developed by other schools.

1.2. Contributions of the Anglo-Saxon School

There are the methods Tucker [43] and its variants Tucker 1 [44,45], Tucker 2 [45,46,47], Tucker 3 [43,47], PARAFAC/CANDECOMP [48,49], INDSCAL [49], Three Mode Scaling [44], IDIOSCAL [50], DEDICOM [51], TUCKALS2 and TUCKALS3 [47], CANDELINC [52], PARATUCK2 [53], PARARLIND [54], CONFAC [55], CP-RIDGE [56] and CP-LASSO [57].

1.3. Contributions of the Salmantina School

There are the methods METABLIPOT [58,59], Multiple Biplot [60], CANOSTATIS [61], ACPR Three ways [62], Triadic Biplot [63], STATIS Dual Canónico [64], Biplot Consensus [65], Co-Tucker [66,67], Coupled Data AnalysisT3-PCA [68], Dynamic Biplot [69,70], Co-Tucker 3 [71,72] and C_enet Tucker [73].

STATIS methods are applied in different fields research as medicine [74,75], policy [76], quality control [77,78,79,80], sensory profiles [42,81,82,83], economy [84,85], customer research [86], education [87], quality of life [88], molecular biology [89,90,91] and other fields.

Among the numerous research works developed in similar fields to this study, these methods can be employed in environmental studies [41,66,92,93,94,95,96,97,98], ecology [99,100], sustainability [101,102,103] and hydrology [104,105].

This work proposes a multivariate alternative called MixSTATICO in order to analyse together the information provided by quantitative and qualitative variables, which allows us to identify, through a space–time analysis, the effect of variations of the physical and chemical variables of the ocean on certain phytoplanktonic species from a perspective similar to that of the STATICO method.

2. Materials and Methods

2.1. Materials

The study area consists of four sampling stations, located spatially from the north to the south of the coastal profile of Continental Ecuador, located 16 km from the coastline. The samples sites of Emeralds and Puerto Bolívar (Pto. Bolivar) are influenced by an estuarine system, while those of Manta and La Libertad maintain a direct influence of the Pacific Ocean. The area has two climatic epochs: rainy (December–May) and dry (June–November).

Data were collected in 2013. At each station, multiple environmental and biological variables were collected monthly at seven depth levels of about 0, 10, 20, 30, 40, 50 and 75 m in the morning (8:00 a.m. to 12:00 p.m.) from February to December in 2013. For a further description of the data, consult González-Narváez [106].

2.1.1. Hydrography and Sampling

A vertical profile of conductivity–temperature–depth (CTD), SeaBird 9/11-plus was measured daily at a fixed position, providing temperature (°C) and salinity (psu) profiles. Water samples were collected at the four sampling sites, using a cast of 3 L nontoxic sampling Van Dorn bottles to obtain water samples for dissolved oxygen concentration (mg L⁻¹) and analysis of nutrients including nitrite, nitrate, phosphate and silicate (µgat L⁻¹), as well as phytoplankton examination. To measure dissolved oxygen, the volumetric method was applied (PEE/LAB-DOQ/01) based on Eaton et al. [107].

2.1.2. Determination of Inorganic Nutrients and Ratios

Samples for nutrient analysis were filtered through washed glass-fibre filters with 0.45 µm (Whatman, GF/C). They were frozen for analysis 1–3 weeks later onshore. Ammonium, nitrite, nitrate, phosphate and silicic acid amounts were determined by the method of Strickland and Parsons [108].

2.1.3. Phytoplankton Analysis

Phytoplankton analysis was performed with approximately 250 cm³ water samples in glass bottles fixed with Lugol’s solution. >A 25 cm³ composite chamber was subsequently filled with the sample water and its contents allowed to settle for 24 h. At least two transects of the chamber bottom were observed with an inverted microscope [109] at 400× magnification to count the small, frequently occurring phytoplankton forms. Additionally, the whole chamber bottom was examined at 125× magnification to count larger, less frequent cells [110,111]. Data are expressed as cells per liter. Classification was performed at the genus or species level when possible, but many taxa could not be identified and were pooled into categories such as “small flagellates” or “small dinoflagellates”. Abbreviations were added for statistical reporting (see Table 1 for species lists); references used can be found in Jiménez [112], Pesantes [113], Balech [114], Tomas [115]), Taylor et al. [116], Tomas [117,118], Young et al. [119] and Jiménez [120].

Appendix B and Appendix C show the average values at 20 m (m) the depth that was initially analysed. The data are of mixed type and thus suitable for testing this method. There are two series of three-way matrices, ${\mathbb{X}}_{\left[\mathit{n}\mathit{p}\mathit{K}\right]}$(environment) and ${\mathbb{Y}}_{\left[\mathit{n}\mathit{q}\mathit{K}\right]}$ (species), each with k = 3 (time—months). The individuals (rows) are the average values obtained within the first 20 m at four sampling stations, while the variables (columns) are the physical and chemical parameters and the abundances of phytoplankton species.

The species matrices include species abundance for Stauroneis membranacea (ordinal-qualitative: 0 if $y_{\left[ij\right]}<Q_1$; 1 if $Q_1\le y_{\left[ij\right]}<Q_2$; 2 if $Q_2\le y_{\left[ij\right]}<Q_3$ and 3 if $y_{\left[ij\right]}\ge Q_3$ for quartiles calculated using all data collected throughout the entire study period). As the last variable (nominal type), the location is used to indicate the sector on which the sampling station is situated: north for Esmeraldas, centre north for Manta, centre south for La Libertad and south for Pto. Bolívar. The seasons were indicated as “rainy from December to May” and “dry from June to November”.

2.2. Methodological Description of the MixSTATICO Method

Boldface lowercase letters (e.g., ${\stackrel{\rightharpoonup }{\mathit{v}}}_{\left[\mathit{n},\mathit{j},\mathit{k}\right]}$) will be used to represent vectors of dimensionality n corresponding to j-th column and k-th condition. Boldface uppercase letters ${\mathit{X}}_{\left[\mathit{n}\mathit{p},\mathit{k}\right]}$ will denote $\left(n\times p\right)$ matrices corresponding to the k-th condition. A set of K such matrices (a “data cube”) will be denoted by a boldface uppercase double-stroke letter such as ${\mathbb{X}}_{\left[\mathit{n}\mathit{p}\mathit{K}\right]}$. The positions of matrix elements will be represented by boldface lowercase letters $\{{\mathit{x}}_{\mathit{i}\mathit{j}\mathit{k}}|i=1,\dots ,n;j=1,\dots p,k=1,\dots ,K\}$.

The transpose of a matrix is written ${\mathit{X}}_{\left[\mathit{n}\mathit{p},\mathit{k}\right]}^{\mathit{T}}$, its rank is $Rank({\mathit{X}}_{\left[\mathit{n}\mathit{p},\mathit{k}\right]}$), its diagonalization matrix is $Diag\left({\mathit{X}}_{\left[\mathit{n}\mathit{p},\mathit{k}\right]}\right)$ and its vectorised form is $Vec\left({\mathit{X}}_{\left[\mathit{n}\mathit{p},\mathit{k}\right]}\right)$. The merger of two matrices is written as $\mathit{X}{\mathit{Y}}_{\left[\mathit{n}{\mathit{r}}^{\prime },\mathit{k}\right]}$, where there are n individuals and $r^{\prime }=p+q$ variables. The identity matrix is ${\mathit{I}}_{\left[\mathit{n}\mathit{n}\right]}$ and the unit vector is $1_{\left[\mathit{n}\right]}$. The multiplication of a matrix by a scalar ${\mathit{\lambda }}_{\left[\mathbf{11},\mathit{k}\right]}$ is represented by ${\mathit{\lambda }}_{\left[\mathbf{11},\mathit{k}\right]}{\mathit{X}}_{\left[\mathit{n}\mathit{p},\mathit{k}\right]}$, while the product of two matrices is symbolised by ${\mathit{X}}_{\left[\mathit{n}\mathit{p},\mathit{k}\right]}\times X_{\left[\mathit{p}\mathit{n},\mathit{k}\right]}$.

Consider two sets of three-way matrices (data cubes) ${\mathbb{X}}_{\left[\mathit{n}\mathit{p}\mathit{K}\right]}$ and ${\mathbb{Y}}_{\left[\mathit{n}\mathit{q}\mathit{K}\right]}$, structured by mixed data that may be quantitative (discrete, continuous) and/or qualitative (binary, ordinal, nominal). The data cubes have the characteristic that their matrices contain the same variables and that the individuals are the same between pairs of matrices for the k-th condition, that is:

(1)${\mathbb{X}}_{\left[\mathit{n}\mathit{p}\mathit{K}\right]}=\left[{\mathit{X}}_{\left[\mathit{n}\mathit{p},\mathbf{1}\right]} | {\mathit{X}}_{\left[\mathit{n}\mathit{p},\mathbf{2}\right]} \left|,\dots ,\right| {\mathit{X}}_{\left[\mathit{n}\mathit{p},\mathit{K}\right]}\right],$

(2)${\mathbb{Y}}_{\left[\mathit{n}\mathit{q}\mathit{K}\right]}=\left[{\mathit{Y}}_{\left[\mathit{n}\mathit{q},\mathbf{1}\right]} | {\mathit{Y}}_{\left[\mathit{n}\mathit{q},\mathbf{2}\right]} \left|,\dots ,\right| {\mathit{Y}}_{\left[\mathit{n}\mathit{q},\mathit{K}\right]}\right].$

Subsequently, it is suggested to pre-treat ${\mathbb{X}}_{\left[\mathit{n}\mathit{p}\mathit{K}\right]}$ and ${\mathbb{Y}}_{\left[\mathit{n}\mathit{q}\mathit{K}\right]}$ as follows:

• Quantitative variables must be applied a mathematical transformation as the suggest principally Kroonenberg [121] and Legendre and Legendre [122], such as: normalization by column, centre by column, logarithm;

• Binary variables must be coded with 0 and 1;

• Ordinal variables must be coded for each scale in ascending numerical order;

• Nominal variables must be transformed to a disjunctive matrix (each category is a new dichotomous variable).

With these modifications, the dimensions of the two data cubes are temporarily modified. The dimensions p and q belonging to the variables will increase by the creation of the disjunctive matrix for each nominal variable, resulting in:

(3)${\mathbb{X}}_{\left[\mathit{n}{\mathit{p}}^{\prime }\mathit{K}\right]}^{\prime }=\left[{\mathit{X}}_{\left[\mathit{n}{\mathit{p}}^{\prime },\mathbf{1}\right]} | {\mathit{X}}_{\left[\mathit{n}{\mathit{p}}^{\prime },\mathbf{2}\right]} \left|,\dots ,\right| {\mathit{X}}_{\left[\mathit{n}{\mathit{p}}^{\prime },\mathit{K}\right]}\right],$

(4)${\mathbb{Y}}_{\left[\mathit{n}{\mathit{q}}^{\prime }\mathit{K}\right]}^{\prime }=\left[{\mathit{Y}}_{\left[\mathit{n}{\mathit{q}}^{\prime },\mathbf{1}\right]} | {\mathit{Y}}_{\left[\mathit{n}{\mathit{q}}^{\prime },\mathbf{2}\right]} \left|,\dots ,\right| {\mathit{Y}}_{\left[\mathit{n}{\mathit{q}}^{\prime },\mathit{K}\right]}\right].$

Among the methods commonly used to analyse the relationships between two sets of variables measured on the same individuals are CCA, Redundancy Analysis (RDA) [123], Canonical Correlation Analysis (CANCOR) [11] and Co-Inertia Analysis (Co-IA) [124]. Co-Inertia differs from the others in that it focuses on finding axes maximizing the covariance between the variables of the two tables.

To implement the methodology of analysing the co-structure between two sets of variables measured on the same individuals, a suggestion made by Abdi et al. [30] for the distance matrix based DISTATIS method is adopted. For mixed data, the distance metric of Gower [125] is used, between the variables of the paired matrices ${\mathit{X}}_{\left[\mathit{n}\mathit{p},\mathit{k}\right]}$ and ${\mathit{Y}}_{\left[\mathit{n}\mathit{q},\mathit{k}\right]}$, having previously merged these matrices in each K-condition to obtain the following set of K-distance matrices of dimension $\left(q^{\prime }\times p^{\prime }\times K\right)$:

(5)$\mathbb{X}{\mathbb{Y}}_{\left[\mathit{n}{\mathit{r}}^{\prime }\mathit{K}\right]}=\left[\mathit{X}{\mathit{Y}}_{\left[\mathit{n}{\mathit{r}}^{\prime },\mathbf{1}\right]} | \mathit{X}{\mathit{Y}}_{\left[\mathit{n}{\mathit{r}}^{\prime },\mathbf{2}\right]} \left| ,\dots , \right| \mathit{X}{\mathit{Y}}_{\left[\mathit{n}{\mathit{r}}^{\prime },\mathit{K}\right]}\right] , \mathrm{where} r^{\prime }=p+q.$

The Gower distance d_ij to each K-matrix is calculated independently from:

(6)$d_{ij}^2=1-s_{ij},$

where

(7)$s_{ij}=\frac{{{\displaystyle \sum }}_{h=1}^{p_1}\left(1-\frac{\left|x{y}_{ih}-x{y}_{jh}\right|}{G_h}\right)+a+\alpha }{p_1+\left(p_2-d\right)+p_3}$

Here,

• $s_{ij}$ is the Gower similarity coefficient;

• $p_1$ is the number of continuous quantitative variables;

• $p_2$ is the number of binary variables;

• $p_3$ is the number of qualitative non-binary variables;

• $a$ and $d$ are the number of matches in binary variables (1,1) and (0,0), respectively;

• $\alpha$ is the number of matches in qualitative non-binary variables; and

• $G_h$ is the rank of the h-th quantitative variable.

The analysis is based only on the distances between the variables p and q in ${\mathit{X}}_{\left[\mathit{n}\mathit{p},\mathit{k}\right]}$ and ${\mathit{Y}}_{\left[\mathit{n}\mathit{q},\mathit{k}\right]}$. Therefore, when at least one of these is nominal, the average Gower distance between the distances of the non-nominal variable and each of the categories of the nominal variable that were previously transformed into dichotomous variables can be obtained from:

(8)${\overline{d}}_{v,w}=\frac{{{\displaystyle \sum }}_{i=1}^{nc}{d}_{v,w_i}}{nc},$

where

• $d$ is the Gower distance;

• $v$ represents the non-nominal variable;

• $w$ represents the nominal variable; and

• $nc$ is the category number of the nominal variable.

If both variables are nominal, then the centroid method is applied, using the generalised average of all distances.

Finally, a set of K-distance matrices of dimension $\left(q\times p\times K\right)$ is obtained that expresses the common structure of the two initial data cubes:

(9)${\mathbb{Z}}_{\left[\mathit{q}\mathit{p}\mathit{K}\right]}=\left[{\mathit{Z}}_{\left[\mathit{q}\mathit{p},\mathbf{1}\right]} | {\mathit{Z}}_{\left[\mathit{q}\mathit{p},\mathbf{2}\right]} \left| ,\dots , \right| {\mathit{Z}}_{\left[\mathit{q}\mathit{p},\mathit{K}\right]}\right].$

Continuing to parallel the approach taken by Abdi et al. [30] to DISTATIS, a normalisation of the distance matrices, using a similar method to that applied in Factorial Multiple Analysis (FMA) [21,126], is developed to obtain new weighted co-inertia matrices ${\mathit{Z}}_{\left[\mathit{q}\mathit{p},\mathit{k}\right]}^{\prime }$ comparable in importance to the original distance matrices when comparing the K-matrices in the factorial plane.

To perform the normalisation, the Singular Value Decomposition (SVD) is applied to each K-distance matrix with rank L. This technique allows us to decompose rectangular matrices ${\mathit{Z}}_{\left[\mathit{q}\mathit{p},\mathit{k}\right]}$ into factors ${\mathit{U}}_{\left[\mathit{q}\mathit{L}\right]}$ and ${\mathit{V}}_{\left[\mathit{p}\mathit{L}\right]}$, matrices containing as their columns orthonormal singular vectors associated with the singular values ${\gamma }_{\mathit{i}\mathit{i}}$ contained in the principal diagonal of the ${\mathsf{\Gamma }}_{\left[\mathit{L}\mathit{L}\right]}$ matrix. ${\mathit{U}}_{\left[\mathit{q}\mathit{L}\right]}$ and ${\mathit{V}}_{\left[\mathit{p}\mathit{L}\right]}$ have the property $({\mathit{U}}_{\left[\mathit{q}\mathit{L},\mathit{k}\right]}^{\mathit{T}}\times {\mathit{U}}_{\left[\mathit{q}\mathit{L},\mathit{k}\right]})=({\mathit{V}}_{\left[\mathit{p}\mathit{L},\mathit{k}\right]}^{\mathit{T}}\times {\mathit{V}}_{\left[\mathit{p}\mathit{L},\mathit{k}\right]})=\mathit{I}$. The singular vectors of ${\mathit{U}}_{\left[\mathit{q}\mathit{L}\right]}$ and ${\mathit{V}}_{\left[\mathit{p}\mathit{L}\right]}$ are eigenvectors of the matrices $({\mathit{Z}}_{\left[\mathit{q}\mathit{p},\mathit{k}\right]}\times {\mathit{Z}}_{\left[\mathit{q}\mathit{p},\mathit{k}\right]}^{\mathit{T}})$ and $({\mathit{Z}}_{\left[\mathit{q}\mathit{p},\mathit{k}\right]}^{\mathit{T}}\times {\mathit{Z}}_{\left[\mathit{q}\mathit{p},\mathit{k}\right]}),$ respectively, with eigenvalues ${\mathit{\lambda }}_{\left[\mathbf{11},\mathit{k}\right]}={\mathit{\gamma }}_{\left[\mathit{i}\mathit{i},\mathit{k}\right]}^{\mathbf{2}}$. The matrices are normalised by multiplying their entries by the scalar ${\mathit{\lambda }}_{\left[\mathbf{11},\mathit{k}\right]}^{-\mathbf{1}}$.

Formally,

(10)$\mathrm{SVD}({\mathit{Z}}_{\left[\mathit{q}\mathit{p},\mathit{k}\right]})=U\mathsf{\Gamma }{V}^T, \mathrm{where} U\in {ℝ}^{qxL}, V\in {ℝ}^{pxL}, \mathsf{\Gamma }\in {ℝ}^{LxL}, \mathrm{and} L=Rank({\mathit{Z}}_{\left[\mathit{q}\mathit{p},\mathit{K}\right]});$

and the normalised K-matrices are:

(11)${\mathit{Z}}_{\left[\mathit{q}\mathit{p},\mathit{k}\right]}^{\prime }={\mathit{\lambda }}_{\left[\mathbf{11},\mathit{k}\right]}^{-\mathbf{1}}{\mathit{Z}}_{\left[\mathit{q}\mathit{p},\mathit{k}\right]}.$

The ${\mathit{Z}}_{\left[\mathit{q}\mathit{p},\mathit{k}\right]}^{\prime }$ form a set of weighted K-distance matrices:

(12)${\mathbb{Z}}_{\left[\mathit{q}\mathit{p}\mathit{K}\right]}^{\prime }=\left[{\mathit{Z}}_{\left[\mathit{q}\mathit{p},\mathbf{1}\right]}^{\prime } | {\mathit{Z}}_{\left[\mathit{q}\mathit{p},\mathbf{2}\right]}^{\prime } \left|,\dots ,\right| {\mathit{Z}}_{\left[\mathit{q}\mathit{p},\mathit{K}\right]}^{\prime }\right].$

The next step is to calculate the cross-product matrix ${\mathit{C}}_{\left[\mathit{K}\mathit{K}\right]}$, known as the variance-and-covariance matrix. It is first necessary to obtain the vectorisation of the matrices ${\mathit{Z}}_{\left[\mathit{q}\mathit{p},\mathit{k}\right]}^{\prime }$:

(13)${\mathbb{Z}}_{\left[{\mathit{s}}^{\prime }\mathit{K}\right]}^{″}={\left[Vec\left( Z_{\left[\mathit{q}\mathit{p},\mathbf{1}\right]}^{\prime }\right),Vec\left( {\mathit{Z}}_{\left[\mathit{q}\mathit{p},\mathbf{2}\right]}^{\prime }\right),\dots ,Vec\left( {\mathit{Z}}_{\left[\mathit{q}\mathit{p},\mathit{K}\right]}^{\prime }\right) \right]}^T ; \mathrm{where} s^{\prime }=q\times p.$

Then,

(14)${\mathit{C}}_{\left[\mathit{K}\mathit{K}\right]}={\mathbb{Z}}_{\left[{\mathit{s}}^{\prime }\mathit{K}\right]}^{″}\times {\mathbb{Z}}_{\left[{\mathit{s}}^{\prime }\mathit{K}\right]}^{{}^{″}\mathit{T}}.$

To compute the inter-structure, Equation (14) is used as well as to calculate a scalar product matrix that contains the Rayleigh coefficients RV [127], the vectorial correlation matrix:

(15)$\mathit{R}\mathit{V}\left(k,k^{\prime }\right)=\frac{{\mathit{C}}_{\left[\mathit{k}{\mathit{k}}^{\prime }\right]}}{{\left({\mathit{C}}_{\left[\mathit{k}\mathit{k}\right]}\right)}^{1/2} {\left({\mathit{C}}_{\left[{\mathit{k}}^{\prime }{\mathit{k}}^{\prime }\right]}\right)}^{1/2}} ; \mathrm{with} k=1,\dots ,K \mathrm{and} k^{\prime }=1,\dots ,K$

This matrix shows the proximity between K-matrices; its interpretation is similar to that of the Pearson linear correlation coefficient.

Abdi et al. [35] apply SVD to the ${\mathit{C}}_{\left[\mathit{K}\mathit{K}\right]}$ matrix to analyse the similarities between K-matrices because it meets the symmetry requirements and is positive defined. Nevertheless, Abdi et al. [30] suggest calculating the eigenvectors and eigenvalues of the $\mathit{R}{\mathit{V}}_{\left[\mathit{K}\mathit{K}\right]}$ matrix (the method adopted in this investigation), such that:

(16)$\mathit{R}{\mathit{V}}_{\left[\mathit{K}\mathit{K}\right]}=U\mathsf{\Xi }{V}^T, \mathrm{with} U\in {ℝ}^{KxL}, V\in {ℝ}^{KxL}, \mathsf{\Xi }\in {ℝ}^{LxL} \mathrm{and} L=Rank(\mathit{R}{\mathit{V}}_{\left[\mathit{K}\mathit{K}\right]})$

After the decomposition of the $\mathit{R}{\mathit{V}}_{\left[\mathit{K}\mathit{K}\right]}$ matrix and the computation of its eigenvectors and eigenvalues (the PCA method), the K-matrices can be projected onto the factorial plane as points, allowing the similarities between them to be studied. The elements of the first eigenvector have the same sign for being a positive semi-definite matrix.

The first eigenvector of ${\mathit{U}}_{\left[\mathit{K}\mathit{K}\right]}$ supplies the weight vector ${\stackrel{\rightharpoonup }{\mathit{\alpha }}}_{\left[\mathit{K}\right]}$, subject to ${{\displaystyle \sum }}_{\mathit{k}=\mathbf{1}}^{\mathit{K}}{\stackrel{\rightharpoonup }{\mathit{\alpha }}}_{\left[\mathit{k}\right]}=\mathbf{1}$:

(17)${\mathit{U}}_{\left[\mathit{k}\mathit{k}\right]}=\left[\begin{array}{ccc}{u}_{11}& \dots & u_{1K}\\ ⋮& \ddots & ⋮\\ u_{K1}& \dots & u_{KK}\end{array}\right]=>\left(\begin{array}{c}{u}_{11}\\ ⋮\\ u_{K1}\end{array}\right)={\stackrel{\rightharpoonup }{\mathit{u}}}_{\left[\mathit{K},\mathbf{1}\right]}$

(18)${\stackrel{\rightharpoonup }{\mathit{\alpha }}}_{\left[\mathit{K}\right]}={\left({\stackrel{\rightharpoonup }{\mathit{u}}}_{\left[\mathit{K},\mathbf{1}\right]}^{\mathit{T}}\times 1_{\left[\mathit{K}\right]}\right)}^{-1} {\stackrel{\rightharpoonup }{\mathit{u}}}_{\left[\mathit{K},\mathbf{1}\right]}$

(19)${\stackrel{\rightharpoonup }{\mathit{\alpha }}}_{\left[\mathit{K}\right]}=\left(\begin{array}{c}{\alpha }_1\\ ⋮\\ {\alpha }_K\end{array}\right)=>{\displaystyle \sum }_{k=1}^K{\alpha }_k=1$

Abdi et al. [35] mention that K-matrix similarities are represented in the first eigenvector ${\stackrel{\rightharpoonup }{\mathit{u}}}_{\left[\mathit{K},\mathbf{1}\right]}$ which provides greater importance at matrices that better represent the co-inertia.

To show the inter-structure graphically, one calculates the coordinates of the K-matrices using the first two vectors of ${\mathit{U}}_{\left[\mathit{K}\mathit{K}\right]}$ (the first columns of which must be positive; see Abdi et al. [35]). From these, one obtains:

(20)$G={\mathit{U}}_{\left[\mathit{K}\mathbf{2}\right]}\times {\mathit{D}}_{\left[\mathbf{22}\right]}{}^{1/2}, \mathrm{where} {\mathit{D}}_{\left[\mathit{L}\mathit{L}\right]}=Diag\left(\mathsf{\Xi }\right)$

With the ${\stackrel{\rightharpoonup }{\mathit{\alpha }}}_{\left[\mathit{k}\right]}$ values, the consensus matrix is calculated, which represents the stable part or average image of the relations between the paired matrices (co-inertia).

(21)${\mathit{W}}_{\left[\mathit{q}\mathit{p}\right]}={\displaystyle \sum }_{k=1}^K{\stackrel{\rightharpoonup }{\mathit{\alpha }}}_{\left[\mathit{k}\right]}{{\mathit{Z}}^{\prime }}_{\left[\mathit{q}\mathit{p},\mathit{k}\right]}$

The consensus matrix ${\mathit{W}}_{\left[\mathit{q}\mathit{p}\right]}$ is analysed using Principal Coordinates Analysis (PcoA) [128] with a double centre. To represent its configuration graphically in a low dimension space, it is necessary to apply SVD:

(22)${\mathit{W}}_{\left[\mathit{q}\mathit{p}\right]}=U^{\prime }{\mathsf{\Xi }}^{\prime }{V}^{{\prime }^T}, \mathrm{where} U^{\prime }\in {ℝ}^{q\times L}, V^{\prime }\in {ℝ}^{p\times L}, {\mathsf{\Xi }}^{\prime }\in {ℝ}^{L\times L} \mathrm{and} L=Rank\left({\mathit{W}}_{\left[\mathit{q}\mathit{p}\right]}\right)$

To generate the graph of the consensus matrix, the scores ${\mathit{F}}_{\left[\mathit{q}\mathit{L}\right]}$ for the row variables and the loadings ${\mathit{Q}}_{\left[\mathit{p}\mathit{L}\right]}$ for the column variables are calculated:

(23)${\mathit{F}}_{\left[\mathit{q}\mathit{L}\right]}={{\mathit{U}}^{\prime }}_{\left[\mathit{q}\mathit{L}\right]}\times {\mathit{D}}_{\left[\mathit{L}\mathit{L}\right]}^{\mathbf{1}/\mathbf{2}} \mathrm{and} {\mathit{Q}}_{\left[\mathit{p}\mathit{L}\right]}={{\mathit{V}}^{\prime }}_{\left[\mathit{p}\mathit{L}\right]}\times {\mathit{D}}_{\left[\mathit{L}\mathit{L}\right]}^{\mathbf{1}/\mathbf{2}} ; \mathrm{with} {{\mathit{D}}^{\prime }}_{\left[\mathit{L}\mathit{L}\right]}=Diag\left({\mathsf{\Xi }}^{\prime }\right)$

In the same way, the intra-structure is studied, which shows the Euclidean image of the consensus and consists of projecting the co-inertia of the K-matrices onto the consensus axes. A double centring is carried out on each K-matrix of ${{\mathit{Z}}^{\prime }}_{\left[\mathit{q}\mathit{p},\mathit{k}\right]}$, and the FK (scores for the row variables) and QK (loadings for the column variables) are obtained:

(24)$FK={{\mathit{Z}}^{\prime }}_{\left[\mathit{q}\mathit{p},\mathit{k}\right]}\times {\mathit{V}}_{\left[\mathit{p}\mathit{L}\right]}^{\prime } \mathrm{and} QK={{\mathit{Z}}^{\prime }}_{\left[\mathit{q}\mathit{p},\mathit{k}\right]}^{\mathit{T}}\times {\mathit{U}}_{\left[\mathit{q}\mathit{L}\right]}^{\prime }$

The procedure of the method is summarized in Figure 1 and the algorithm is shown in Appendix A.

The algorithm that summarizes the MixSTATICO method steps is in Appendix A.

2.3. Case Study

In Ecuador, there is limited information regarding previous studies carried out on this subject. Therefore, it is necessary to do a space–time analysis to identify the effect that typical environmental variations have on this region. Increases in temperature and nutrient concentrations can occur, owing to the effect of the geographical location and the climatic season of this area on the abundance of phytoplanktonic, which is predominately the diatoms group.

The coast of Continental Ecuador is located in the eastern equatorial Pacific (Figure 2), with the north influenced by the warm current of Panama (heading south and producing an increase in SST and decrease in nutrients; begins in December and intensifies in February and April) and the south by the cold stream of Humboldt (heading north, its waters are coldest, high in salinity and rich in nutrients; it mainly presents in July and November, weakening in December), and there are two climatic seasons (dry and rainy) [129]. The southern part is considered to be the most significant tropical estuary on the west coast of South America [120]. However, much of its coastline is contained in the Niño Region 1 + 2 (0°–10° S, 90°–80° W), an area sensitive to the variations produced by the surface sea temperature from the Central Pacific due to the presence of ENSO events (El Niño—La Niña), which are events that alter the physical and chemical parameters of the sea. On the west coast of South America, the most evident signs of ENSO are manifested in ocean–atmospheric systems and by the impact it has on natural ecosystems [130].

Due to the importance of this area, constant monitoring of the chemical and biological conditions of the sea is carried out in space and in time.

3. Results

All data were analysed using the software R Studio version 1.2.1335 and the functions daisy (Gower distance) and svd (Singular Value Decomposition).

Below is a summary of the spatial–temporal data describing the environmental and biological characteristics of the coastal region of Ecuador (see Appendix B and Appendix C for all the environmental and species data).

Table 1 and Table 2 show the average value and standard deviation (SD) of the environmental and biological variables, respectively, according to the location of the sampling stations and the climatic period.

The stations in the north (Esmeraldas and Manta) presented the highest average values with regard to temperature and the lowest average values for nutrients. In addition, the environmental variables (except salinity) showed the highest average values during the rainy season and the Esmeraldas estuary had the warmest waters of the entire coast (Table 1).

The phytoplankton species had different average abundance values that varied according to the geographical location of each station (which had different environmental characteristics) and the influence of the season (Table 2).

The inter-structure graph, the result of the first step of the MixSTATICO analysis, showed that the environment–species association had two types of behaviour, which varied according to the season (Figure 3A). The members of the first type presented a similar behaviour during the dry season (red in Figure 3A), while those of the second type behaved similarly during the months of the rainy season (blue), except for October and December. The two first components explained 93.07% of the inertia (Figure 3B).

Moreover, there was a solid typical structure between the environment and species, as it was observed that all K-matrices (months) reached values close to the unitary cycle. They also had high values for weights. Therefore, all matrices substantially contributed toward building the consensus matrix (Figure 3A,C). November contributed the most to the consensus (weight × 0.094), followed by August and September (dry season), but the consensus of the environment–species association was more notably represented during November (highest value of cos²) (Figure 3C).

April, February and May (rainy season) made the lowest contributions to the consensus, while May was considered a month of seasonal change and the end of the rainy season.

In relation to the consensus analysis, which showed the stable part of the species–environment association, the two first components explained 99.92% of the inertia (Figure 4C). The consensus for the environment variables showed a strong positive correlation among temperature, dissolved oxygen and salinity, and a strong negative correlation with phosphate, silicate and nitrate (Figure 4A).

In 2013, the average temperature value was 23.6 °C, with the highest values recorded at the stations located to the north, while the highest average nutrient values were recorded at the stations south of the coastal profile (Table 1).

The ordination of species in the consensus (Figure 4B) established a high level of association between S. membranacea (e23—continuous variable) and S. membranacea (e24ordinal variable), where both variables contain the same information, showing the suitability of the method when different variable types are presented. The abundance of centric diatoms Dactyliosolen antarcticus (e4), Proboscia alata (e5), Skeletonema costatum (e6), Lauderia borealis (e7), Chaetoceros didymus (e15) and Ditylum brightwellii (e21) and dinoflagellate Gyrodinium sp. (e19) were associated (same direction) with the location and season variables. This indicated that the vulnerability of the species varied according to the geographic location of the sampling stations and the season, with different levels of average abundance (Table 2). In addition, these species had the lowest average abundance throughout the whole coastline (Total column of Table 2), being absent during some of the months and stations (Table A2). In addition, the species with stable abundance levels (ordered in the opposite direction to the location and season variables), during both periods and in the entire coastal profile, were mainly P. pungens (e10), Gymnodinium sp. (e18) and Dactyliosolen fragilissimus (e1), S. membranacea (e23) and S. membranacea (e24).

The vulnerable species indicated in Figure 4B (located on the factorial environmental plane Figure 4A) show preferential environmental conditions with below-average values in relation to temperature, dissolved oxygen and salinity (opposite direction to these variables).

The following species where those with the highest abundance (longer vector) over the entire coastline: the pennate diatoms Pseudo-nitzschia pungens (e10) and N. longissima (e8) and the dinoflagellate Gymnodinium sp. (e18). These species preferred above-average values in relation to temperature, dissolved oxygen and salinity. Dactyliosolen antarcticus (e4), on the other hand, was detected during the dry season at the stations located in the north (Esmeraldas and Manta) (Figure 3A,B).

The trajectories of both variables and species for the rainy and the dry seasons are presented separately for the purpose of comparison.

The same interpretation is applied to the intra-structure.

During the rainy months (Figure 5A,C,E,G,I), it was observed that the environmental and species variables showed a non-stable management pattern. For February and March, temperature, dissolved oxygen and salinity showed a direct and inverse association with phosphate and nitrate, while in April and December, that association became weaker, mainly being between temperature and dissolved oxygen. The association pattern in May (month of change of season) differs from the rest of the months since it was observed that temperature does not directly associate with any other variable. At this time, the highest average values for temperature and nutrients were recorded (Table 1).

The species with the highest abundance during the rainy season (Figure 5B,D,F,H,J) were mainly the species Leptocylindrus danicus (e11), Chaetoceros affinis (e13) and Gymnodinium sp. (e18). In February to April (they showed a longer vector), their high abundance was associated with above-average temperature and salinity values and low abundance was associated with nitrate values, among other nutrients.

Each sampling month showed the species that had a decrease in their abundance or disappeared due to the influence of the sampling station location and the variability of the corresponding climatic feature in that month. At this time of year (except in December), the main vulnerable species were the centric diatoms D. antarcticus (e4), P. alata (e5) and S. costatum (e6), which were found mainly at the southern station, and D. brightwellii (e21) at the stations in the centre of the coastal profile that were influenced directly by the ocean. These species seemed to thrive more in waters with below-average temperature and salinity values (Table 2).

During the dry season, like during the months of the rainy season, a non-stable environmental ratio pattern (Figure 6A,C,E,G,I,K) was displayed, changing according to the month the samples were taken. In addition, lower values were obtained for nutrients and temperature.

The species that stood out as being more abundant during these months were centric diatoms D. antarcticus (e4), S. costatum (e6), P. pungens (e10) and D. brightwellii (e21), among others, and their preferred environmental conditions varied according to the month. The species that were vulnerable to environmental variations and geographical location were mainly Lauderia Borealis (e7), Chaetoceros curvisetus (e14), C. didymus (e15), Gyrodinium sp. (e19) and Navicula sp. (e22), as they had the lowest average abundance values (see Table 2). However, the species that had a steady abundance were Guinardia striata (e2), Rhizosolenia imbricata (e3), D. antarcticus (e4), P. alata (e5), S. membranacea (e24) and S. membranacea (e23) (Figure 6B,D,F,H,J,L).

Moreover, there are more apparent patterns of association between species in the dry season (except in June, the start of the dry season) as compared to the rainy season.

4. Discussion

Water sustainability is an overriding issue for ecosystem support [1].

Human actions and natural events are factors that cause alterations in the environment, affecting ecosystem biotic and abiotic components [132,133,134]. Therefore, planktonic components are of great importance in assessing the water quality of aquatic systems, as they are subject to environmental variations [135] in the area and the geographical locations where they grow [95]. This was established by the analysis of the consensus and the intra-structure, where phytoplankton species vulnerable to climatic variability and environmental conditions of the estuarine and oceanic stations were identified. These stations are seen to have a positive effect on the ecosystem, which we believe gives support to the importance of studying the hydrological variations and the effects they have on biota and biotic processes. The estuarine and oceanic stations are known as dynamic environments that exhibit spatial and temporal variations in their biotic and abiotic parameters [136]. This dynamic behaviour was classified into two groups based on the months of the year (inter-structure), where the rainy season was established as February–May and the dry season as June–December [129], except for October and December. December was the start of the rainy season of 2014, the year in which a weak El Niño event occurred [137] with aforementioned anomalies >1.5 (Oceanic Niño Index—ONI) that decreased as the month passed.

In addition to the intra-structure, differences in species population patterns were observed between both seasons, with the dry season showing a more pronounced pattern of association among the species. The reason for this is that during the rainy season, the species–environment co-structure was weaker and the environmental variables presented unstable behaviour due to the presence of high outlier values for phosphate, nitrate, nitrite and salinity during February, March and April (in the different sample stations). This could be due to an increased amount of runoff originating from sedimentary basins (agriculture fields wash off). However, regarding the temperature, it only increases in August, a month that tends to experience upwellings along the Ecuador coast at the northern station (Esmeraldas) [106]. The highest average temperature and the lowest nutrient values were recorded in the stations located to the north (Esmeraldas and Manta). In contrast, an opposite effect occurred in the southern stations (La Libertad and Pto. Bolivar) (see Table 2). This behaviour coincides with the presence of warm waters coming from the north caused by the influence of the Panama current. However, in the south, oceanic waters are influenced by the cold Humboldt current, as well as upwellings that favour planktonic productivity [111].

During most of the months of the dry season (intra-structure), the species D. antarcticus (e4) showed a high abundance, which was associated with above-average nitrate values and low temperatures [138], which occurred mainly in June and July (Table A2). This is due to the presence of a cold Humboldt current, rich in nutrients, which is consistent with optimal environmental conditions to support the high productivity of diatoms. In the intra-structure chart for February, March, May and December, the dinoflagellate Gyrodinium sp. (e19) was found to be the species that associated its above-average abundance with geographical location (north, Esmeraldas and Manta—Table 2) and climate (rainy—Table 2). In 2013, other authors [139] stated that nutrient concentrations are the drivers of dinoflagellate productivity, drawing attention to species such as Gymnodinium cf catenatum, Oxytoxum turbo and Prorocentrum micans, which we found were associated with warm waters (typical of the Panama current) with low salinity and nutrients during the rainy season in the northern station. This may have occurred because dinoflagellates trophic behaviour is not dependent on nutrient levels, and as such is considered an indicator of species of warm ocean waters with low nutrient levels [140,141].

In the intra-structure, it was possible to identify the high abundance species, such as P. alata (e5) in November, which was associated with above-average nitrate and temperature values in the north station, Esmeraldas (Table 2). In August, the species D. fragilissimus (e1) was abundant with above-average values in temperature and oxygen association values in the north centre station, Manta (Table 2), but whose presence was found to be associated with nutrient availability [137]. In February, March, April and May, the abundance of P. pungens (e10) was found to be associated with above-average temperature and oxygen values in the south centre station, La Libertad (Table 2), which is also consistent with the warmer temperatures associated with conditions during the rainy season. Finally, throughout the rainy season (Feb–Dec), the species Hemiaulus sinensis (e16) associated its high abundance to above-average nitrate, nitrite and salinity values in the south station, Pto. Bolívar. These results reaffirm the reports of other studies [106] that describe how P. alata (e5) is the most abundant species in Esmeraldas and D. fragilissimus (e1) in Manta during the dry season, and P. pungens (e10) in La Libertad and H. sinensis (e16) in Pto. Bolívar during the rainy season. Nevertheless, the findings of this study highlight the association of phytoplankton growth with location and seasons, where P. alata (e5) was found to be associated with these variables. However, in the case of D. fragilissimus (e1), P. pungens (e10) and H. sinensis (e16), these species were found to have had a stable level of abundance throughout all the seasons and along the entire coastal profile.

This study has shown the utility of this new multivariate method in analysing associations between biota abundance and environmental variables such as location and season within ENSO neutral conditions. In addition, this work highlights the seasonal patterns and possible association with freshwater run-off, which proves that the method could be useful for the assessment of other environmental contexts.

5. Conclusions

In this work, using the MixSTATICO multi-way mixed multivariate method, the need for additional analysis has been reduced. In this sense, descriptive statistics techniques such as PCA and CCA were employed and applied independently to the data collected at each sampling station. In addition, this work introduces geographical location and climatic periods as qualitative variables and analyses them in conjunction with the other quantitative variables, namely temperature, dissolved oxygen, salinity, nutrients and phytoplankton abundance.

Our findings support the effectiveness of the proposed method, which reduces the amount of analysis required to obtain the same results. This is conducted through the analysis of data cubes that contain sets of quantitative and qualitative variables. In this sense, the STATICO method calculates the stable part of the common structure between two data cubes.

This innovative method has identified the phytoplankton species showing poor, stable and high levels of abundance according to climatic variability and the environmental conditions typical of these geographical locations and seasons. In addition, this work has identified the more vulnerable species in the coast profile as D. antarcticus, P. alata, S. costatum, L. borealis, C. didymus and D. brightwellii (centric diatoms) and Gyrodinium sp. (dinoflagellate). During the dry season, it was possible to observe more apparent patterns of association between species.

These results provide proof of the validity of this method for analysing mixed data, and even offer the possibility of introducing qualitative characteristics to the study. In this sense, the results obtained are further enriched for the knowledge of the ecosystem. This is performed from a quantitative and qualitative analysis perspective. Furthermore, this method, supported by statistical data, efficiently introduces the insight of environmental assessment in terms of environmental parameters’ relationships.

Author Contributions

Conceptualization, M.J.F.-G. and O.R.-B.; methodology, M.G.-N. and O.R.-B.; software, M.G.-N.; validation, S.M.; formal analysis, M.J.F.-G.; investigation, S.M. and M.J.F.-G.; resources, S.M., M.J.F.-G., P.G.-V. and J.-L.M.; data curation, M.G.-N.; writing—original draft preparation, M.G.-N.; writing—review and editing, P.G.-V., J.-L.M., S.M. and M.J.F.-G.; visualization, M.G.-N.; supervision, P.G.-V., M.J.F.-G., S.M. All authors have read and agreed to the published version of the manuscript.

Funding

This research received no external funding.

Institutional Review Board Statement

Not applicable.

Informed Consent Statement

Not applicable.

Data Availability Statement

The data presented in this study are available in Appendix B and Appendix C.

Acknowledgments

We acknowledge the invaluable collaboration of INOCAR (Oceanographic Institute of the Navy, Ecuador), CPFG-EM Edwin Pinto Uscocovich, and Gladys Torres, who have provided us with oceanographic data from a project that studied El Niño under the regional project framework El Niño-ERFEN-ECUADOR. We have no business and/or financial interest in the project; we have disclosed our interests fully and have in place an approved plan for managing any potential conflicts arising from this arrangement.

Conflicts of Interest

The authors declare no conflict of interest.

Appendix A. MixSTATICO Algorithm

• Input two pre-processed data-cubes ${\mathbb{X}}_{\left[\mathit{n}\mathit{p}\mathit{K}\right]}$ and ${\mathbb{Y}}_{\left[\mathit{n}\mathit{q}\mathit{K}\right]}$.

• Merge all ${\mathit{X}}_{\left[\mathit{n}\mathit{p},\mathit{k}\right]}$ and ${\mathit{Y}}_{\left[\mathit{n}\mathit{q},\mathit{k}\right]}$ matrices.

  For k = 1 to K

  $\left\{\mathit{X}{\mathit{Y}}_{\left[{\mathit{n}\mathit{r}}^{\prime },\mathit{k}\right]}=merge\left({\mathit{X}}_{\left[\mathit{n}\mathit{p},\mathit{k}\right]},{\mathit{Y}}_{\left[\mathit{n}\mathit{q},\mathit{k}\right]}\right) \right\} ; \mathrm{where} r^{\prime }=p+q$

• Compute the Gower distance, considering Equations (6)–(8).

  For k = 1 to K

  $\{ {\mathbb{Z}}_{\left[\mathit{q}\mathit{p}\mathit{K}\right]}=gower\left(\mathit{X}{\mathit{Y}}_{\left[{\mathit{n}\mathit{r}}^{\prime },\mathit{k}\right]}\right) \}$

• Apply FMA to normalize

  For k = 1 to K

  $\mathrm{SVD} ({\mathit{Z}}_{\left[\mathit{q}\mathit{p},\mathit{k}\right]}) \mathrm{to} \mathrm{obtain} {\mathit{\lambda }}_{\left[\mathbf{11},\mathit{k}\right]} \mathrm{and} \mathrm{multiply} \mathrm{each} \mathrm{matrix} \mathrm{element} \mathrm{by} {\mathit{\lambda }}_{\left[\mathbf{11},\mathit{k}\right]}^{-1}$

  $\{ {\mathbb{Z}}_{\left[\mathit{q}\mathit{p}\mathit{K}\right]}^{\prime }={\mathit{\lambda }}_{\left[\mathbf{11},\mathit{k}\right]}^{-\mathbf{1}}{\mathit{Z}}_{\left[\mathit{q}\mathit{p},\mathit{k}\right]}\}$

• Vectorize ${\mathbb{Z}}_{\left[\mathit{q}\mathit{p}\mathit{K}\right]}^{\prime }$

  ${\mathbb{Z}}_{\left[{\mathit{s}}^{\prime }\mathit{K}\right]}^{″}={\left[Vec\left({{\mathit{Z}}^{\prime }}_{\left[\mathit{q}\mathit{p},\mathbf{1}\right]}\right), Vec\left({{\mathit{Z}}^{\prime }}_{\left[\mathit{q}\mathit{p},\mathbf{2}\right]}\right),\dots ,Vec\left({{\mathit{Z}}^{\prime }}_{\left[\mathit{q}\mathit{p},\mathit{K}\right]}\right)\right]}^T ; \mathrm{where} s^{\prime }=q\times p$

• Calculate ${\mathit{C}}_{\left[\mathit{K}\mathit{K}\right]}$

  ${\mathit{C}}_{\left[\mathit{K}\mathit{K}\right]}={\mathbb{Z}}_{\left[{\mathit{s}}^{\prime }\mathit{K}\right]}^{″}\times {\mathbb{Z}}_{\left[s^{\prime }\mathit{K}\right]}^{″\mathit{T}}$

• Calculate $\mathit{R}{\mathit{V}}_{\left[\mathit{K}\mathit{K}\right]}$.

  For k = 1 to K

  $\{\mathit{R}{\mathit{V}}_{\left[\mathit{K}\mathit{K}\right]}=\frac{{\mathit{C}}_{\left[\mathit{k}{\mathit{k}}^{\prime }\right]}}{\sqrt{{\mathit{C}}_{\left[\mathit{k}\mathit{k}\right]}}\times \sqrt{{\mathit{C}}_{\left[{\mathit{k}}^{\prime }{\mathit{k}}^{\prime }\right]}}} \}; \mathrm{with} k=1,\dots ,K \mathrm{and} k^{\prime }=1,\dots ,K$

• SVD $\left(\mathit{R}{\mathit{V}}_{\left[\mathit{K}\mathit{K}\right]}\right)$ to obtain ${\stackrel{\rightharpoonup }{\mathit{\alpha }}}_{\left[\mathit{K}\right]}$

• Calculate consensus matrix

  ${\mathit{W}}_{\left[\mathit{q}\mathit{p}\right]}={\displaystyle \sum }_{k=1}^K{\stackrel{\rightharpoonup }{\mathit{\alpha }}}_{\left[\mathit{k}\right]}{{\mathit{Z}}^{\prime }}_{\left[\mathit{q}\mathit{p},\mathit{k}\right]}$

• SVD $({\mathit{W}}_{\left[\mathit{q}\mathit{p}\right]}$).

• Generate consensus graph using PCoA.

• Generate intra-structure graph using PCoA.

Appendix B. Table Environmental Data

Table A1

Text Values for physical and chemical variables.

Space	Time	Temperature (T—°C)	Salinity (S—psu)	Dissolved Oxygen (DO—mg.L−1)	Nitrate(NO3—µg−at.L−1)	Nitrite(NO2—µg−at.L−1)	Phosphate($\mathbf{P}{\mathbf{O}}_4^{3^{-}}$—µg−at.L−1)	Silicate($\mathbf{S}\mathbf{i}{\mathbf{O}}_4^{4^{-}}$—µg−at.L−1)
Esmeraldas	Feb	26.74	32.65	5.29	0.77	0.08	0.71	10.42
Manta	Feb	24.02	33.24	4.35	6.91	0.16	1.62	5.12
La Libertad	Feb	25.25	33.11	5.38	0.73	0.28	0.58	11.88
Pto. Bolivar	Feb	22.23	33.47	4.35	8.76	0.22	1.31	7.52
Esmeraldas	Mar	21.92	34.11	4.28	5.95	0.39	2.00	14.36
Manta	Mar	23.04	33.69	4.37	2.73	0.43	1.75	13.73
La Libertad	Mar	23.58	34.29	4.62	3.64	0.30	1.33	5.24
Pto. Bolivar	Mar	21.17	31.56	3.35	9.36	0.29	0.97	12.87
Esmeraldas	Apr	27.45	33.24	4.45	0.50	0.06	0.71	11.06
Manta	Apr	24.47	34.15	4.36	1.55	0.28	0.38	8.07
La Libertad	Apr	22.79	34.47	4.46	4.96	0.09	1.26	10.08
Pto. Bolivar	Apr	23.12	33.56	4.16	3.86	0.73	1.02	8.14
Esmeraldas	May	26.11	33.21	4.60	1.80	0.19	0.70	3.22
Manta	May	26.07	33.40	4.35	3.85	0.15	0.75	2.10
La Libertad	May	24.65	33.64	4.12	1.70	0.56	1.17	19.10
Pto. Bolivar	May	24.95	33.81	4.53	4.92	0.02	0.51	10.66
Esmeraldas	Jun	26.33	33.17	4.60	2.02	0.20	0.38	10.11
Manta	Jun	24.31	33.95	4.32	0.65	0.34	0.31	1.63
La Libertad	Jun	20.45	34.41	3.20	7.75	0.41	0.95	9.55
Pto. Bolivar	Jun	22.18	34.22	5.12	0.36	0.11	0.33	2.65
Esmeraldas	Jul	25.50	32.91	4.73	1.13	0.03	1.02	3.66
Manta	Jul	23.90	33.56	4.62	0.99	0.10	0.84	3.16
La Libertad	Jul	22.68	33.78	4.82	0.39	0.04	0.83	9.34
Pto. Bolivar	Jul	22.80	33.99	5.01	1.47	0.24	1.25	8.87
Esmeraldas	Aug	25.49	32.89	4.62	0.24	0.01	0.25	1.92
Manta	Aug	22.73	33.59	4.47	2.27	0.05	0.71	6.73
La Libertad	Aug	19.09	34.38	3.89	4.04	0.23	1.00	12.21
Pto. Bolivar	Aug	21.60	33.93	3.90	3.19	0.04	0.56	12.46
Esmeraldas	Sep	26.33	32.68	4.32	0.72	0.04	0.26	1.16
Manta	Sep	21.35	34.20	4.53	1.09	0.01	0.38	8.84
La Libertad	Sep	21.85	33.96	4.00	2.43	0.09	0.58	5.01
Pto. Bolivar	Sep	22.54	33.84	4.41	3.36	0.25	0.50	4.66
Esmeraldas	Oct	26.09	32.20	4.23	0.58	0.18	0.10	3.73
Manta	Oct	24.31	33.43	3.81	2.59	0.23	0.18	5.11
La Libertad	Oct	23.18	33.72	4.51	0.29	0.13	0.33	6.87
Pto. Bolivar	Oct	22.60	33.81	4.17	1.54	0.38	0.46	11.88
Esmeraldas	Nov	26.06	32.41	4.39	0.12	0.08	0.11	2.23
Manta	Nov	24.93	32.99	4.56	0.30	0.05	0.13	2.82
La Libertad	Nov	20.30	34.38	3.27	5.03	0.26	0.61	9.11
Pto. Bolivar	Nov	21.21	34.24	3.58	6.08	0.31	0.63	10.89
Esmeraldas	Dec	25.45	32.65	3.96	0.26	0.03	0.21	6.24
Manta	Dec	24.38	33.09	4.48	1.03	0.03	0.37	1.27
La Libertad	Dec	21.64	34.23	4.25	3.52	0.29	0.69	7.41
Pto. Bolivar	Dec	21.56	34.23	3.13	3.81	0.37	0.09	8.48

Appendix C. Table Species Data

Table A2

Text The abundance of the phytoplankton species.

Space	Time	Phytoplankton Species 1—Abundance Cells L−1																								Localization	Sea-son
		e1	e2	e3	e4	e5	e6	e7	e8	e9	e10	e11	e12	e13	e14	e15	e16	e17	e18	e19	e20	e21	e22	e23	e24
Esmeraldas	Feb	0.00	4.62	3.74	0.00	3.80	0.00	0.00	4.36	0.00	2.90	4.88	0.00	0.00	0.00	0.00	0.00	0.00	3.94	4.44	3.74	0.00	0.00	0.00	0	North	Rainy
Manta	Feb	0.00	3.37	3.85	0.00	0.00	0.00	0.00	3.59	0.00	4.04	3.67	0.00	3.67	0.00	0.00	0.00	0.00	2.90	0.00	0.00	0.00	0.00	0.00	0	Centre North	Rainy
Libertad	Feb	0.00	4.46	0.00	0.00	0.00	0.00	0.00	0.00	0.00	0.00	4.37	0.00	3.59	0.00	3.74	3.67	0.00	0.00	0.00	3.37	0.00	0.00	0.00	0	Centre South	Rainy
Pto. Bolivar	Feb	0.00	3.37	3.50	0.00	0.00	0.00	0.00	4.84	0.00	5.09	3.59	0.00	5.18	0.00	0.00	0.00	2.90	3.80	0.00	4.55	0.00	4.07	0.00	0	South	Rainy
Esmeraldas	Mar	4.33	4.17	4.66	0.00	2.90	0.00	0.00	4.62	0.00	4.69	4.36	3.20	3.74	3.80	0.00	3.37	3.20	3.80	0.00	3.50	0.00	0.00	4.13	3	North	Rainy
Manta	Mar	4.01	4.28	4.15	3.20	3.20	0.00	0.00	5.58	2.90	4.54	3.74	4.34	4.90	4.55	3.50	4.10	3.37	4.22	4.24	3.80	2.90	2.90	4.10	3	Centre North	Rainy
Libertad	Mar	5.35	4.52	3.80	0.00	2.90	2.90	0.00	4.52	0.00	0.00	4.28	4.41	2.90	0.00	0.00	3.20	0.00	4.22	3.50	3.67	0.00	0.00	3.50	2	Centre South	Rainy
Pto. Bolivar	Mar	3.80	4.26	0.00	0.00	0.00	6.30	0.00	5.06	3.20	5.45	4.56	0.00	5.55	5.03	4.04	3.67	3.80	4.46	0.00	2.90	3.20	3.37	3.74	2	South	Rainy
Esmeraldas	Apr	0.00	0.00	0.00	3.94	0.00	0.00	0.00	5.10	0.00	3.37	3.67	0.00	3.50	0.00	0.00	0.00	4.13	3.50	3.37	4.01	0.00	2.90	0.00	0	North	Rainy
Manta	Apr	4.69	4.52	0.00	4.31	0.00	0.00	0.00	5.51	3.37	4.01	5.22	3.74	5.25	0.00	4.29	4.10	4.04	3.80	2.90	3.37	0.00	3.50	4.31	3	Centre North	Rainy
Libertad	Apr	0.00	0.00	0.00	0.00	0.00	0.00	0.00	4.49	0.00	0.00	0.00	0.00	0.00	0.00	0.00	0.00	0.00	4.46	0.00	4.67	0.00	3.74	0.00	0	Centre South	Rainy
Pto. Bolivar	Apr	0.00	0.00	0.00	0.00	0.00	0.00	0.00	4.17	4.24	3.97	0.00	2.90	0.00	0.00	0.00	0.00	0.00	0.00	2.90	3.20	0.00	0.00	3.37	2	South	Rainy
Esmeraldas	May	0.00	4.51	2.90	0.00	2.90	4.78	0.00	4.45	3.37	4.59	4.07	2.90	5.46	4.58	0.00	4.73	3.85	0.00	3.20	3.67	3.67	0.00	0.00	0	North	Rainy
Manta	May	4.01	4.22	3.59	0.00	0.00	0.00	0.00	3.94	4.41	0.00	2.90	0.00	4.47	3.67	0.00	0.00	0.00	0.00	2.90	3.59	0.00	2.90	0.00	0	Centre North	Rainy
Libertad	May	3.80	3.80	0.00	0.00	0.00	0.00	0.00	0.00	0.00	2.90	0.00	3.80	0.00	0.00	0.00	2.90	0.00	2.90	0.00	2.90	0.00	0.00	3.20	1	Centre South	Rainy
Pto. Bolivar	May	0.00	0.00	0.00	0.00	0.00	0.00	0.00	4.13	0.00	0.00	0.00	2.90	0.00	0.00	0.00	0.00	0.00	0.00	3.20	4.53	0.00	0.00	0.00	0	South	Rainy
Esmeraldas	Jun	4.94	5.25	4.34	4.74	3.74	4.72	0.00	3.74	0.00	4.15	4.87	0.00	4.34	0.00	0.00	5.88	3.50	0.00	0.00	3.37	3.20	0.00	4.04	3	North	Dry
Manta	Jun	4.76	5.03	5.12	3.90	4.89	0.00	0.00	4.47	2.90	0.00	4.53	3.20	4.73	0.00	0.00	4.75	3.37	2.90	0.00	3.20	0.00	0.00	0.00	0	Centre North	Dry
Libertad	Jun	4.67	5.22	4.41	0.00	4.62	0.00	0.00	4.61	2.90	0.00	4.40	4.01	4.58	0.00	0.00	4.93	3.94	3.80	0.00	2.90	3.74	0.00	4.01	3	Centre South	Dry
Pto. Bolivar	Jun	5.01	6.16	4.01	5.06	4.99	0.00	4.34	4.61	3.59	4.15	4.95	3.20	4.10	0.00	0.00	5.87	0.00	3.67	0.00	4.01	3.20	0.00	4.10	3	South	Dry
Esmeraldas	Jul	4.86	5.14	4.99	4.49	4.41	4.28	3.67	4.40	3.74	3.74	4.68	0.00	4.49	5.27	0.00	4.44	0.00	3.90	0.00	0.00	3.37	0.00	3.20	1	North	Dry
Manta	Jul	2.90	0.00	4.29	4.78	0.00	0.00	0.00	4.07	3.37	3.20	4.34	3.59	4.20	0.00	0.00	0.00	3.20	4.04	0.00	3.20	2.90	2.90	4.13	3	Centre North	Dry
La Libertad	Jul	4.53	5.33	4.71	5.79	3.20	3.59	0.00	3.67	2.90	0.00	3.50	4.24	4.41	0.00	0.00	0.00	0.00	3.59	0.00	2.90	0.00	0.00	3.67	2	Centre South	Dry
Pto. Bolivar	Jul	4.34	5.16	4.01	5.15	0.00	4.26	4.52	3.80	3.74	5.39	4.53	4.26	4.26	4.53	0.00	0.00	4.41	3.20	0.00	3.20	3.37	3.20	0.00	0	South	Dry
Esmeraldas	Aug	4.39	5.04	4.66	4.37	4.13	3.94	0.00	4.28	3.97	4.46	4.51	3.50	4.56	3.94	0.00	3.67	3.67	4.01	0.00	3.20	2.90	0.00	3.67	2	North	Dry
Manta	Aug	4.86	5.22	5.84	5.57	4.24	4.40	0.00	4.64	3.85	4.39	5.18	3.80	4.67	3.80	0.00	4.45	4.36	3.97	0.00	0.00	0.00	3.50	4.45	3	Centre North	Dry
La Libertad	Aug	4.07	4.94	5.83	5.71	4.71	0.00	4.13	3.94	3.20	3.67	5.01	4.51	4.31	0.00	0.00	4.68	4.65	4.37	3.20	3.80	0.00	2.90	0.00	0	Centre South	Dry
Pto. Bolivar	Aug	0.00	3.50	4.40	0.00	0.00	0.00	0.00	3.85	3.20	0.00	3.59	3.20	0.00	0.00	0.00	0.00	0.00	4.04	0.00	4.22	0.00	0.00	0.00	0	South	Dry
Esmeraldas	Sep	3.37	4.93	4.66	3.50	4.26	0.00	0.00	4.07	3.50	0.00	4.39	0.00	4.44	0.00	0.00	3.20	3.20	4.04	0.00	3.37	0.00	3.90	3.80	2	North	Dry
Manta	Sep	4.53	5.02	4.95	3.94	3.94	0.00	0.00	4.92	3.94	3.74	5.12	0.00	4.41	0.00	0.00	3.37	0.00	3.90	0.00	3.90	2.90	3.59	3.67	2	Centre North	Dry
La Libertad	Sep	0.00	4.33	4.82	2.90	0.00	0.00	0.00	4.41	3.20	4.10	4.34	0.00	4.01	0.00	0.00	0.00	0.00	4.17	0.00	3.37	0.00	2.90	0.00	0	Centre South	Dry
Pto. Bolivar	Sep	0.00	3.37	3.20	0.00	0.00	0.00	0.00	4.59	3.67	0.00	0.00	2.90	3.80	0.00	0.00	0.00	3.37	4.46	0.00	3.90	0.00	0.00	4.40	3	South	Dry
Esmeraldas	Oct	0.00	3.85	3.97	0.00	3.20	0.00	0.00	3.37	3.59	0.00	3.37	3.20	0.00	0.00	0.00	0.00	0.00	4.04	0.00	3.59	0.00	3.59	0.00	0	North	Dry
Manta	Oct	4.10	4.69	3.80	4.10	3.59	4.29	0.00	3.74	3.85	4.22	4.39	3.59	4.40	0.00	0.00	3.67	3.50	4.51	0.00	3.74	0.00	3.20	3.50	2	Centre North	Dry
La Libertad	Oct	0.00	0.00	0.00	0.00	0.00	0.00	0.00	5.13	3.85	3.37	0.00	4.01	0.00	0.00	0.00	3.67	3.50	4.44	0.00	3.50	0.00	3.80	0.00	0	Centre South	Dry
Pto. Bolivar	Oct	0.00	0.00	0.00	0.00	0.00	0.00	0.00	3.85	2.90	0.00	0.00	3.67	0.00	0.00	0.00	0.00	0.00	4.51	0.00	3.37	0.00	2.90	0.00	0	South	Dry
Esmeraldas	Nov	4.46	5.03	4.20	4.55	4.50	0.00	0.00	3.80	3.97	4.45	4.62	4.07	3.80	0.00	0.00	4.43	3.90	4.33	0.00	3.90	0.00	4.36	4.31	3	North	Dry
Manta	Nov	4.22	5.08	4.07	3.80	4.40	5.94	0.00	4.67	4.39	4.07	5.62	3.97	5.09	3.80	0.00	4.69	4.40	4.15	0.00	3.74	0.00	3.90	3.85	2	Centre North	Dry
La Libertad	Nov	4.45	5.31	6.14	4.07	3.50	0.00	4.37	5.12	4.24	4.15	5.12	3.50	4.57	0.00	0.00	3.80	4.26	4.31	0.00	3.94	0.00	3.37	0.00	0	Centre South	Dry
Pto. Bolivar	Nov	3.80	0.00	0.00	0.00	0.00	0.00	0.00	4.87	3.85	3.20	4.13	3.50	3.74	0.00	0.00	3.20	4.31	3.37	0.00	3.20	0.00	2.90	0.00	0	South	Dry
Esmeraldas	Dec	4.24	4.64	3.20	4.75	4.20	3.97	3.20	0.00	3.85	3.37	4.78	3.94	4.26	0.00	0.00	3.85	0.00	3.97	0.00	3.85	2.90	3.80	4.26	3	North	Rainy
Manta	Dec	4.65	5.07	4.26	4.68	3.80	0.00	4.44	4.15	4.50	0.00	4.59	3.74	4.57	0.00	0.00	4.33	3.97	4.13	0.00	3.59	0.00	4.15	0.00	0	Centre North	Rainy
La Libertad	Dec	4.91	5.09	5.21	4.28	4.24	4.52	4.43	5.35	5.35	4.31	4.75	4.26	4.59	0.00	0.00	4.77	3.74	3.74	0.00	0.00	3.20	4.13	3.85	2	Centre South	Rainy
Pto. Bolivar	Dec	4.40	5.43	4.61	0.00	4.22	0.00	4.01	5.25	4.29	3.74	4.75	4.01	4.15	0.00	0.00	3.74	3.85	3.94	0.00	3.90	0.00	4.49	0.00	0	South	Rainy

¹Dactyliosolen fragilissimus (e1), Guinardia striata (e2), Rhizosolenia imbricata (e3), Dactyliosolen antarcticus (e4), Proboscia alata (e5), Skeletonema costatum (e6), Lauderia borealis (e7), Nitzschia longissima (e8), Nitzschia sp. (e9), Pseudo-nitzschia pungens (e10), Leptocylindrus danicus (e11), Thalassiosira sp. (e12), Chaetoceros affinis (e13), Chaetoceros curvisetus (e14), Chaetoceros didymus (e15), Hemiaulus sinensis (e16), Thalassionema nitzschioides (e17), Gymnodinium sp. (e18), Gyrodinium sp. (e19), Mesodinium rubrum (e20), Ditylum brightwellii (e21), Navicula sp. (e22), Stauroneis membranacea (e23—continuous variable), Stauroneis membranacea (e24—ordinal variable).

Publisher’s Note: MDPI stays neutral with regard to jurisdictional claims in published maps and institutional affiliations.

Figures and Tables

Figure 1. Scheme of the MixSTATICO method. For this application, XnpK is the environment and YnpK is the species.

Figure 2. Localization of the sample stations—1. Esmeraldas (0.92° N, 80.11° W), 2. Manta (0.88° S, 80.83° W), 3. La Libertad (2.08° S, 81.12° W), 4. Pto. Bolívar (3.11 ° S, 80.49° W). The stations 2, 3 and 4 are contained in Region Niño 1+2. Self-authorship using Ocean Data View [131].

Figure 3. (A) Inter-structure graph (dry season: red colour; rainy season: blue colour). (B) Eigenvalues and (C) weights (α) and Cos2 for each K matrix.

Figure 4. Consensus: factorial plane to (A) environment and (B) species; (C) eigenvalues. Species names: Dactyliosolen fragilissimus (e1), Guinardia striata (e2), Rhizosolenia imbricata (e3), Dactyliosolen antarcticus (e4), Proboscia alata (e5), Skeletonema costatum (e6), Lauderia borealis (e7), Nitzschia longissima (e8), Nitzschia sp. (e9), Pseudo-nitzschia pungens (e10), Leptocylindrus danicus (e11), Thalassiosira sp. (e12), Chaetoceros affinis (e13), Chaetoceros curvisetus (e14), Chaetoceros didymus (e15), Hemiaulus sinensis (e16), Thalassionema nitzschioides (e17), Gymnodinium sp. (e18), Gyrodinium sp. (e19), Mesodinium rubrum (e20), Ditylum brightwellii (e21), Navicula sp. (e22), Stauroneis membranacea (e23—continuous variable), Stauroneis membranacea (e24—ordinal variable).

Figure 5. Intra-structure—rainy season: projections onto factorial plane of the consensus to environment and species corresponding each matrix (time). In the months of the rainy season, the environmental and species variables showed a non-stable management pattern. Distinct environmental correlations and weak patterns of association between species are observed. Species names: Dactyliosolen fragilissimus (e1), Guinardia striata (e2), Rhizosolenia imbricata (e3), Dactyliosolen antarcticus (e4), Proboscia alata (e5), Skeletonema costatum (e6), Lauderia borealis (e7), Nitzschia longissima (e8), Nitzschia sp. (e9), Pseudo-nitzschia pungens (e10), Leptocylindrus danicus (e11), Thalassiosira sp. (e12), Chaetoceros affinis (e13), Chaetoceros curvisetus (e14), Chaetoceros didymus (e15), Hemiaulus sinensis (e16), Thalassionema nitzschioides (e17), Gymnodinium sp. (e18), Gyrodinium sp. (e19), Mesodinium rubrum (e20), Ditylum brightwellii (e21), Navicula sp. (e22), Stauroneis membranacea (e23—continuous variable), Stauroneis membranacea (e24—ordinal variable).

Figure 6. Intra-structure—dry season: projections onto factorial plane of the consensus to environment and species corresponding each matrix (time). In the months of the dry season, there are more apparent patterns of association between species with a non-stable environmental ratio pattern. Species names: Dactyliosolen fragilissimus (e1), Guinardia striata (e2), Rhizosolenia imbricata (e3), Dactyliosolen antarcticus (e4), Proboscia alata (e5), Skeletonema costatum (e6), Lauderia borealis (e7), Nitzschia longissima (e8), Nitzschia sp. (e9), Pseudo-nitzschia pungens (e10), Leptocylindrus danicus (e11), Thalassiosira sp. (e12), Chaetoceros affinis (e13), Chaetoceros curvisetus (e14), Chaetoceros didymus (e15), Hemiaulus sinensis (e16), Thalassionema nitzschioides (e17), Gymnodinium sp. (e18), Gyrodinium sp. (e19), Mesodinium rubrum (e20), Ditylum brightwellii (e21), Navicula sp. (e22), Stauroneis membranacea (e23—continuous variable), Stauroneis membranacea (e24—ordinal variable).

Table 1

Environment variables (average ± SD).

	Location (Both Seasons)				Season (Entire Coastal Profile)		Total
	North	Centre North	Centre South	South	Rainy	Dry
Temperature (T—°C)	25.77 ± 1.34	23.96 ± 1.18	22.31 ± 1.79	22.36 ± 1.03	24.03 ± 1.79	23.24 ± 2.04	23.6 ± 1.97
Salinity (S—psu)	32.92 ± 0.49	33.57 ± 0.38	34.03 ± 0.41	33.7 ± 0.72	33.49 ± 0.68	33.61 ± 0.63	33.56 ± 0.66
Dissolved Oxygen (DO—mg.L−1)	4.5 ± 0.33	4.38 ± 0.2	4.23 ± 0.61	4.16 ± 0.6	4.34 ± 0.49	4.3 ± 0.48	4.32 ± 0.49
Nitrate (${\mathrm{NO}}_{3^{-}}$—µg-at.L−1)	1.28 ± 1.59	2.18 ± 1.81	3.13 ± 2.2	4.25 ± 2.74	3.53 ± 2.6	2.03 ± 1.97	2.71 ± 2.4
Nitrite (${\mathrm{NO}}_{2^{-}}$—µg-at.L−1)	0.12 ± 0.11	0.17 ± 0.13	0.24 ± 0.15	0.27 ± 0.19	0.25 ± 0.18	0.16 ± 0.12	0.2 ± 0.16
Phosphate (${\mathrm{PO}}_4^{3^{-}}$—µg-at.L−1)	0.59 ± 0.53	0.67 ± 0.52	0.85 ± 0.31	0.69 ± 0.37	0.91 ± 0.51	0.53 ± 0.31	0.7 ± 0.45
Silicate (${\mathrm{SiO}}_4^{4^{-}}$—µg-at.L−1)	6.19 ± 4.31	5.33 ± 3.62	9.62 ± 3.76	9.01 ± 3.07	8.85 ± 4.3	6.44 ± 3.65	7.54 ± 4.14

Table 2

Species variables (average ± SD).

Group	Species PhytoplanktonAbundance Cells L−1	Location (Both Seasons)				Season (Entire Coastal Profile)		Total
		North	Centre North	Centre South	South	Rainy	Dry
CD	Dactyliosolen fragilissimus (e1)	2.78 ± 2.14	3.88 ± 1.34	2.89 ± 2.22	1.94 ± 2.15	2.41 ± 2.21	3.26 ± 1.94	2.87 ± 2.11
CD	Guinardia striata (e2)	4.29 ± 1.42	4.23 ± 1.43	3.91 ± 1.9	2.84 ± 2.3	3.52 ± 1.82	4.07 ± 1.91	3.82 ± 1.89
CD	Rhizosolenia imbricata (e3)	3.76 ± 1.34	3.99 ± 1.41	3.17 ± 2.47	2.16 ± 2	2.37 ± 2	4.02 ± 1.65	3.27 ± 1.99
CD	Dactyliosolen antarcticus (e4)	2.76 ± 2.11	3.48 ± 1.74	2.07 ± 2.38	0.93 ± 1.97	1.26 ± 1.94	3.18 ± 2.15	2.31 ± 2.27
CD	Proboscia alata (e5)	3.46 ± 1.22	2.55 ± 1.97	2.11 ± 1.99	0.84 ± 1.78	1.61 ± 1.82	2.76 ± 2	2.24 ± 2.01
CD	Skeletonema costatum (e6)	1.97 ± 2.17	1.33 ± 2.21	1 ± 1.67	0.96 ± 2.08	1.12 ± 2.02	1.48 ± 2.12	1.32 ± 2.08
CD	Lauderia borealis (e7)	0.62 ± 1.33	0.4 ± 1.28	1.17 ± 1.92	1.17 ± 1.91	0.8 ± 1.62	0.88 ± 1.71	0.84 ± 1.67
PD	Nitzschia longissima (e8)	3.84 ± 1.29	4.48 ± 0.64	3.75 ± 1.83	4.46 ± 0.5	3.96 ± 1.74	4.28 ± 0.49	4.13 ± 1.24
PD	Nitzschia sp. (e9)	2.36 ± 1.8	3.41 ± 1.2	2.33 ± 1.88	2.97 ± 1.46	1.97 ± 2.04	3.43 ± 0.83	2.77 ± 1.67
PD	Pseudo-nitzschia pungens (e10)	3.25 ± 1.63	2.93 ± 1.82	2.05 ± 1.9	2.82 ± 2.23	2.85 ± 1.97	2.69 ± 1.95	2.76 ± 1.96
CD	Leptocylindrus danicus (e11)	4.38 ± 0.47	4.48 ± 0.77	3.25 ± 2.03	2.74 ± 2.11	3.41 ± 1.78	3.97 ± 1.59	3.71 ± 1.7
CD	Thalassiosira sp. (e12)	1.89 ± 1.75	2.73 ± 1.69	2.98 ± 1.84	2.77 ± 1.38	2.21 ± 1.85	2.91 ± 1.54	2.59 ± 1.73
CD	Chaetoceros affinis (e13)	3.51 ± 1.72	4.58 ± 0.41	3 ± 1.9	2.8 ± 2.18	3.29 ± 2.01	3.62 ± 1.65	3.47 ± 1.83
CD	Chaetoceros curvisetus (e14)	1.6 ± 2.14	1.44 ± 1.91	0 ± 0	0.87 ± 1.85	1.08 ± 1.89	0.89 ± 1.75	0.98 ± 1.82
CD	Chaetoceros didymus (e15)	0 ± 0	0.71 ± 1.51	0.34 ± 1.08	0.37 ± 1.16	0.78 ± 1.56	0 ± 0	0.35 ± 1.12
CD	Hemiaulus sinensis (e16)	3.05 ± 1.99	3.04 ± 1.9	2.87 ± 1.86	1.5 ± 2.08	2.32 ± 1.94	2.86 ± 2.13	2.62 ± 2.07
PD	Thalassionema nitzschioides (e17)	2.31 ± 1.77	2.75 ± 1.72	1.83 ± 2.02	2.06 ± 1.92	1.84 ± 1.86	2.56 ± 1.85	2.24 ± 1.89
D	Gymnodinium sp. (e18)	3.23 ± 1.53	3.5 ± 1.21	3.64 ± 1.23	3.22 ± 1.57	2.89 ± 1.71	3.82 ± 0.89	3.4 ± 1.41
D	Gyrodinium sp. (e19)	1 ± 1.66	0.91 ± 1.53	0.61 ± 1.29	0.55 ± 1.18	1.53 ± 1.73	0.13 ± 0.64	0.77 ± 1.44
C	Mesodinium rubrum (e20)	3.29 ± 1.07	2.92 ± 1.39	3.18 ± 1.13	3.72 ± 0.55	3.34 ± 1.21	3.23 ± 1.03	3.28 ± 1.12
CD	Ditylum brightwellii (e21)	1.46 ± 1.61	0.79 ± 1.29	0.63 ± 1.34	0.89 ± 1.45	0.79 ± 1.38	1.07 ± 1.52	0.94 ± 1.46
PD	Navicula sp. (e22)	1.69 ± 1.87	2.77 ± 1.36	1.89 ± 1.76	1.9 ± 1.79	2 ± 1.85	2.12 ± 1.68	2.06 ± 1.76
PD	Stauroneis membranacea (e23)	2.49 ± 1.91	2.55 ± 1.94	1.66 ± 1.82	1.42 ± 1.89	1.72 ± 1.92	2.28 ± 1.95	2.03 ± 1.96

CD: Centric Diatom, PD: Pennate Diatom, D: Dinoflagellate, C: Ciliate.