Automatic Segregation of Pelagic Habitats

It remains difficult to segregate pelagic habitats since structuring processes are dynamic on a wide range of scales and clear boundaries in the open ocean are non-existent. However, to improve our knowledge about existing ecological niches and the processes shaping the enormous diversity of marine plankton, we need a better understanding of the driving forces behind plankton patchiness. Here we describe a new machine-learning method to detect and quantify pelagic habitats based on hydrographic measurements. An Autoencoder learns two-dimensional, meaningful representations of higher-dimensional micro-habitats, which are characterized by a variety of biotic and abiotic measurements from a high-speed ROTV. Subsequently, we apply a density-based clustering algorithm to group similar micro-habitats into associated pelagic macro-habitats in the German Bight of the North Sea. Three distinct macro-habitats, a “surface mixed layer,” a “bottom layer,” and an exceptionally “productive layer” are consistently identified, each with its distinct plankton community. We provide evidence that the model detects relevant features like the doming of the thermocline within an Offshore Wind Farm or the presence of a tidal mixing front.

Recent advances in marine remote sensing technology (Wedding et al., 2011) enabled scientists to separate benthic structures into mosaic-like patterns of different habitat classes (Hinchey et al., 2008;Pittman et al., 2011) following the role model of terrestrial ecosystems. However, what is well known and trivial in landscape ecology can be quite challenging in seascape ecology. While it remains difficult to segregate pelagic habitats, which exhibit no clear boundaries (Hinchey et al., 2008;Pittman et al., 2011;Wedding et al., 2011) and can be quite dynamic on a wide range of scales, benthic habitat maps can give an impression of physically distinct areas that consistently occur together with particular species communities (Harris and Baker, 2012). Some effort has been undertaken to characterize fish habitats (e.g., Bellido et al., 2008;Giannoulaki et al., 2011;Tugores et al., 2011;Laman et al., 2017;Amorim et al., 2018;Friedland et al., 2020;Funk et al., 2020), but fewer studies focused on zooplankton (e.g., Labat et al., 2009;Alvarez-Berastegui et al., 2014;Espinasse et al., 2014). Thus, mechanisms contributing to the enormous diversity of plankton, a fundamental component of pelagic food webs, are still not fully understood (Sano et al., 2013;North et al., 2016). Understanding the processes shaping plankton communities is essential to improve our knowledge of existing ecological niches (Houliez et al., 2021). Despite the growing awareness of the importance of spatial structure for ecology and management Wedding et al., 2011), there is still a lack of concepts and techniques applicable to characterize the spatial structure of the seascape in pelagic environments (Alvarez-Berastegui et al., 2014). Mainly, because traditional oceanographic methods are inadequate for observing the submesoscale (Baschek and Maarten Molemaker, 2010) due to insufficient resolution and range (Marmorino et al., 2018). Recent advances in instrumentation partially closed this gap, but there still is a need for novel analysis methods to take advantage of the existing data (North et al., 2016). Some machine learning techniques are specifically designed to identify and characterize features in a "sea of data, " which makes it very promising to apply them also in this challenging field of research.
Autoencoders (AE) are a common tool in the machine learning community which consist of an encoding and a decoding part (Hinton, 2006). Initially devised to reduce (Encoder) and recover (Decoder) the dimensionality of their inputs (Hinton, 2006), they have been soon applied to a wide range of tasks like denoising (e.g., Vincent et al., 2010) or anomaly detection (e.g., Zhao et al., 2017;Chen et al., 2018).
Autoencoders do not classify or detect specific elements or objects in their inputs, but learn meaningful low dimensional representations, i.e., relevant high-level abstractions, of their inputs (Bengio et al., 2006) so that the original data can be reconstructed as similar as possible by the decoder part. The input data don't need any pre-processing, e.g., labeling of subsets, by humans, since the target the network aims to reconstruct is basically the original input. The compressed representations of the encoder can also be used as input for subsequent modeling, e.g., in a Convolutional Neural Network (CNN) application. In that case the unsupervised pre-training of a CNN embedded in an AE can help to capture more intricate dependencies (Erhan et al., 2009) and better initialize the weights of the extended model (Bengio et al., 2006). Thus, the (local) minimum in the loss surface of the AE corresponds to a good transformation of a high dimensional input to a lower dimensional intermediate output (output of the Encoder-part) (Bengio et al., 2006), which would become the input for the classifier in a CNN. In this setting, the final output of the AE, the reconstructions, are secondary. However, a low reconstruction error of the AE ensures that the compressed signal incorporates the important features of the original high dimensional input data.
In this study we take advantage of this specific application of AEs. Instead of substituting the decoder part with a classification or regression network we use the compressed signal of the encoder as input for a subsequent clustering algorithm. We use a fully connected AE to reduce a high dimensional input consisting of a variety of abiotic and biotic oceanographic measurements to a lower dimensional meaningful representation (intermediate output), skip the decoding part after the training is completed and cluster the encoded features to macro-habitats. Similar microhabitats lead to similar representations and therefore regions with different characteristics are segregated as different macrohabitats. These macro-habitats correspond to distinct pelagic habitats in the southern North Sea, whose plankton communities are compared and analyzed.

Data Acquisition and Preparation
Physical and biological oceanographic measurements were recorded on a North Sea summer cruise with the RV Heincke (HE429, July 19-24, 2014) with a MacArtney TRIAXUS Remotely Operated Towed Vehicle (ROTV). For a detailed description of the device see Plonus et al. (2021). The ROTV transects were located in the direct vicinity of two Offshore Wind Farms (OWF) BARD Offshore 1 (BARD) and Global Tech I (GTI; Figure 1). The map was generated using QGIS v3.18 (QGIS, 2021) with bathymetric metadata and Digital Terrain Model data products from the EMODnet Bathymetry portal (15.7.21) 1 . The ROTV was towed at a speed of 8 knots (4.1 m s −1 ) with a three-degree lateral offset to lessen any disturbance from the vessels wake. During most transects the ROTV was undulating with a vertical speed of 0.1 m s −1 from ∼ 4 m below the sea surface to ∼ 8 m above the sea floor. The horizontal resolution between two surface peaks was ∼ 560 m, while the vertical resolution was ∼ 0.3 m. The ROTV measured water temperature, salinity, oxygen, and chlorophyll-a at a frequency of 1 Hz and was equipped with a Video Plankton Recorder (VPR, Seascan Inc., Falmouth, MA, United States) which provided zoo-and phytoplankton densities on the taxonomic family-, and sometimes even genus-level. For a detailed description of the VPR plankton image classification see Floeter et al. (2017). We used a similar summer cruise with the RV Heincke 5 years later (HE534, June 16-21, 2019) as a test data set. For our analyses we selected the following variables: temperature ( • C), salinity (PSU), oxygen (µmol l −1 ), density (kg × m −1 ), and chlorophyll-a (RFU). For each of the variables, we calculated the horizontal (grid cell to the left, i.e., ∼ 25 m) and vertical (grid cell above, i.e., 1 m) gradient. Furthermore, we had sufficient density data (N l −1 ) available for the taxa "Appendicularia, " "Copepoda, " "Dinoflagellates, " "Gastropoda, " "Jelly, " "Marine snow, " "Nauplii, " "Ophiuroida, " "Pilidium, " "Pluteus, " and "Polychaeta." Transect diagrams were generated using Ocean Data View (ODV, Schlitzer, 2020) with the embedded spatial interpolation software DIVA (Troupin et al., 2012) and exported as grids with a resolution of ∼ 25 m length × 1 m depth. Abiotic measurements as well as density values were normalized and rescaled to range from −1 to 1. This was necessary since deep learning models generally perform better with homogeneous, small values (Bishop, 1995).
To check for multicollinearity between our variables we calculated the variance inflation factor (VIF) in R (R Core Team, 2020) using functions provided by Zuur et al. (2009). A threshold of VIF > 3 was applied to identify highly collinear variables and exclude them from further analyses (Zuur et al., 2010). The exported grids for each selected parameter were stacked and transformed into feature-vectors where each grid cell became one vector with four features (1 parameter = 1 feature). In our definition, a pelagic micro-habitat with a spatial extent of ∼25 m × 1 m corresponds to one of those feature-vectors (Figure 2).
Based on these feature-vectors the AE was trained to reconstruct the original micro-habitats and thereby learn relevant abstractions that represent important patterns in the pelagic environment. We used a GPU supported TensorFlow backend (Abadi et al., 2015) for Keras (Chollet, 2015) under Python 3.7 (Van Rossum and Drake, 2009) to build and train our AE.

Model Description
The AE consisted of two fully connected layers in the Encoder and Decoder, respectively. The Decoder used the transposed weights of the Encoder in reversed order, e.g., the weights of the first Encoder-Layer were shared with the last Decoder-Layer. The first layer of the Encoder inflated the 4-dimensional featurevector to a 100-dimensional feature-vector, which was reduced to a 2-dimensional feature-vector by the second layer (4 -100 -2). The Decoder did the same in reverse (2 -100 -4). The batch size (number of inputs that are processed simultaneously) was set to 38 and the learning rate followed a sawtooth-like scheme, initialized at 5e −8 . Each input feature-vector corresponds to one micro-habitat and includes 1 measurement of each parameter selected for the analyses. The model was trained using the data from HE429 exclusively. Approximately ∼ 13% of the data was separated to validate the training process based on the remaining 87%. Data from HE534 was used as a final test set. As an AE is a gradient-based method, the chosen starting point may be crucial for the final fit of the model (Hinton, 2006), and one way of assessing and reducing the effect of start conditions are multi-start approaches (Subbey, 2018). Therefore, we trained multiple models and selected the one with the smallest final validation RMSE.

Habitat Segregation
By applying the trained Encoder only, we projected all microhabitats into a xy-coordinate system using the 2-dimensional intermediate output. We will refer to the encoded outputs as "Encoded Components" in the following. Micro-habitats with similar characteristics were projected closer to each other than micro-habitats with different characteristics. We used the Euclidean distance to calculate the dissimilarity matrix for the Encoded Components of the micro-habitats, which was clustered by the HDBSCAN algorithm (McInnes et al., 2017). HDBSCAN uses a density-based linkage function, defining clusters by the size of the area in which a certain number of neighbors is found. Micro-habitats in "sufficiently dense" regions were assigned to a macro-habitat. Obviously, the parameters "size of the neighborhood" (Epsilon) and the "critical number of neighbors" (min_samples) are determining the resulting clusters FIGURE 2 | Schematic data processing from measurements to feature-vector.
(dense areas) with HDBSCAN. Thus, we checked the resulting macro-habitats for multiple different combinations of these two parameters as well as "min_cluster_size." The parameter "min_cluster_size" is the threshold that separates "sufficiently dense" regions (clusters) from the random background noise. All micro-habitats that were not assigned to a specific macrohabitat by HDBSCAN got the label "−1." Homogeneous regions in the transect produced more dense regions in the 2-dimensional surface that were more likely to trespass the "min_cluster_size" threshold and were separated from other homogeneous water masses by less dense regions. We used the silhouette method (Rousseeuw, 1987) to select the best segregation of microhabitats. The silhouette score ranges from "−1" to "1" and indicates how well each point fits into the assigned cluster (macro-habitat) and is one of the best performing indices available (Arbelaitz et al., 2013). "−1" is probably wrong labeled, "0" is close to the decision boundary of two clusters and "1" means this specific point is far away from points of other clusters. The silhouette scores were calculated using the scikitlearn module (Pedregosa et al., 2011) for python.

Analyses
We used ODV to add a transect plot of the identified macro-habitats to the original measurements and plankton densities. Isolines of selected parameter measurements were overlayed on the macro-habitat plots to investigate which feature characteristics contributed to the segregation and to assess the associated plankton communities.
We furthermore described the macro-habitat plankton communities by modified Species-Abundance-Plots (SAP). We calculated the relative number of micro-habitats by plankton density and taxonomic group for each cluster. As is common for SAPs we used a log2 scale for density. That way we visualized the shift in specific species densities between the macro-habitats of different segregations of the same ROTV survey transect.
Pelagic submesoscale features often are highly productive areas and aggregate particles (Levy and Martin, 2013;Lévy et al., 2018). Therefore, we calculated Lloyd's mean crowding (Lloyd's MC) and Lloyd's index of patchiness (Lloyd's IP) with the R-function "agg_index" from the "epiphy" package (Gigot, 2018) and compared the results for different segregations of the same transects. Lloyd's index is >1 if species were aggregated, 1 if the distribution is random and <1 indicates an overdispersed distribution compared to a homogeneous distribution. The Index of aggregation proposed by Bez (2000) (Bez's IoA) was calculated in addition to Lloyd's IP.

RESULTS
In the initial VIF analysis with the full dataset a couple of parameters had VIF > 3. After removing "density" which had the highest score, no further parameter exceeded this threshold (Supplementary Table 1). After a detailed analysis of model sensitivities and reconstruction quality we decided to limit the final parameter selection to (1) vertical temperature difference to the grid cell above, (2) salinity, (3) oxygen, and (4) chlorophylla concentration.

Model Training
The root mean squared error (RMSE) after the first epoch ranged roughly between 0.7 and 1.0. Each training epoch took 10-15 s using a graphic card with 768 gpu-cores and we trained each model for 15 epochs until a plateau was reached ( Supplementary  Figure 1). The final training and validation RMSE of our selected model were RMSE Tr ∼ 0.33 and RMSE Val ∼ 0.35 (Figure 3).

Clustering
Depending on the HDBSCAN parameter selection, microhabitats were grouped into 2-20 macro-habitats that ranged in size from <0.1 to 97% of all micro-habitats in a transect. We present exemplary the results for the segregation of T3 into different numbers of macro-habitats. Different parameter combinations could lead to an identical number of segregations. We chose an inverse size-cluster-relationship in the figure since more macro-habitats were usually ecologically less plausible (Figure 4). Mostly, "epsilon" had a great impact on the segregation with specific combinations of "min_cluster_size" and "min_samples" but less influence with other tested combinations of those two parameters, indicating that segregations changed discontinuously with slopes and plateaus (Supplementary Figure 2).
Even though the label "−1" is used by HDBSCAN to indicate the lack of belonging to a specific cluster, we observed a close relationship between micro-habitats labeled "−1" and exceptionally strong chlorophyll-peaks throughout all transects. Therefore, we decided to treat "−1" as a macro-habitat of its own instead of unclassified micro-habitats. Micro-habitats labeled as "−1" were also frequently located between the BL and the SL.

Projections
While cluster-labeling was not consistent in that cluster "0" always referred to, e.g., the "surface mixed layer, " the projections of the "surface mixed layer" micro-habitats were always located in a similar position throughout all projection plots. Thus, while the cluster denotations related to a macro-habitat were not consistent, the position indicated the affiliation to a specific macro-habitat (Figure 5).

Silhouette Method
The segregation into three macro-habitats gave the highest average silhouette-scores in most cases: notably high chlorophyllpeaks were merged into one macro-habitat (1) and two further macro-habitats were separated at around 17 • C in an upper surface mixed layer (2) and a lower bottom layer (3). There was only one exception from this rule in T1 where in the northern, deeper area the bottom layer (3) was replaced with the layer including the chlorophyll-peaks (1). Another anomaly occurred in T2, where one of the basic macro-habitats was further separated into two "sublayers" so that a total of four macrohabitats were segregated. The highest silhouette-scores ranged from 0.35 (T1) to 0.59 (T5) ( Table 1).
In the following we will use abbreviations for the three main layers and their sublayers, namely "PL" for the productive layer with the high chlorophyll values, "BL" for the bottom layer and "SL" for the surface mixed layer. A segregation into more than one layer is indicated using numbers, e.g., "SL1"/"SL2" instead of "SL."

Habitat Maps
We present T2 exemplarily for all transects of HE429 (Figure 6). Segregating the output of the Encoder into 3 macro-habitats, we got the typical scheme of a SL with temperatures above 17 • C, a macro-habitat which was strongly associated with extraordinary high chlorophyll-peaks (PL) and a BL as a third macro-habitat. The average silhouette-score for the entire transect was 0.54 (Table 1). However, this clustering did not account, e.g., for the intrusion of marine snow particles into the SL in the eastern half of the transect. When accepting 4 different macro-habitats, the BL and PL macro-habitats were mostly unaffected, while the SL was further separated into 2 different macro-habitats. One corresponded to the area where marine snow particles were predominant while the second macro-habitat corresponded to the area where pluteus larvae were observed in high densities. Segregating characteristics of the two macro-habitats were a salinity difference of 0.2 and a shallowing of the thermocline from 10 m to 5 m water depth. Notably, this change around section distance 18-20 km was located at the entry point of the transect into the Offshore Wind Farm BARD. This segregation increased the average silhouette-score for the entire transect to 0.56 (Table 1).

Species Abundance Plots
The segregation into four macro-habitats was further supported by the modified SAPs. The relative amounts of PL and Frontiers in Marine Science | www.frontiersin.org FIGURE 5 | Projections for different segregations of all transects from HE429. The number of segregated macro-habitats is indicated by the last digit in the panel header. BL, Bottom layer (blue); SL, Surface mixed layer (yellow-red); TC, Thermocline Layer; PL, Productive layer (black). Numbers were used to indicate that, e.g., more than 1 "Bottom layer" was segregated (e.g., B1/B2 instead of BL).
BL did not change much between 3 and 4 macro-habitats. However, SL1 included all micro-habitats with copepod densities >8 N l −1 and basically all micro-habitats where pluteus larvae occurred. SL2 instead included micro-habitats with The numbers in "Segregation" give the number of micro-habitats by macro-habitat, e.g., X_Y_Z indicates 3 macro-habitats with X, Y, and Z micro-habitats, respectively. copepod densities <8 N l −1 and generally less chlorophyll, but most micro-habitats where Appendicularia occurred. Thus, the SL1 and SL2 plankton communities were clearly distinct (Figure 7).

Lloyd
Lloyd's mean crowding underpinned the SAP results. Patchiness in PL and BL did not change for "marine snow" and "pluteus" but differed clearly between SL1 and SL2, indicating a higher pluteus aggregation in SL1 and a higher aggregation of marine snow in SL2 ( Table 2).

Test Dataset HE534
The temperature maximum during HE534 was around 15 • C, i.e., 2 • C lower than the threshold that separated SL and BL in HE429. Consequently, no thermal stratification was detected by the Encoder trained with HE429 measurements. However, this model segregated an oxygen-rich layer that, based on the projections of the Encoder, resembled a similar habitat as the SL in HE429. This oxygen driven stratifications were not consistent over the entire range of a transect and in some areas the macrohabitat with projections similar to BL in HE429 comprised the entire water column, indicating a mixed water column closer toward the coast. Notably, plankton aggregations were commonly located at the border between oxygen-stratified and mixed water columns (Figure 8). The highest average silhouette-scores were reached with three segregated macro-habitats. However, the scores were much lower compared to HE429 with a maximum between 0.26 and 0.40.   Presence gives the relative number of micro-habitats with density >0. ms, marine snow; plu, pluteus.

Selection of Parameters
When training the model, we got the best results with a limited selection of parameters compared to the entire set of available data. The selected parameters are, however, in accordance with previous findings that physical properties contribute most to differences in habitat utilization by plankton organisms (Schulz et al., 2012;Friedland et al., 2020). In contrast to Alvarez-Berastegui et al. (2014), we did not benefit from the combination of gradients with the original measurements. However, a prior wavelet analysis as in North et al. (2016) could help to identify relevant spatial scales for the derivation of gradients. It is also possible, that the architecture of the model limited the amount of compressed information accessible to the clustering algorithm. In convolutional AEs, the size of the bottleneck (intermediate output) limits the generalization of the model (Manakov et al., 2019). This is also true for the fully connected AE architecture of this model and might limit the potential of including more variables like species densities and environmental gradients.

Reconstruction Loss
The loss for the optimization of an AE is based on the difference between the reconstruction and the original input. However, driving forces behind habitat partitioning vary with study region and season and specific parameters have a higher contribution than others (Schulz et al., 2012;Espinasse et al., 2014;Friedland et al., 2020). Thus, we deemed it more important to accurately reconstruct specific features (parameters) compared to entire vectors (micro-habitats). Accordingly, we calculated the sum of the batchwise RMSE between the specific featurevalues (e.g., temperature) of each input and the corresponding feature-values of the reconstructions and not the RMSE of an entire feature-vector (micro-habitat) and its reconstruction. This forced the AE to learn all parameters individually and furthermore made it possible to give specific parameters a higher priority if appropriate.

Aggregation
Lloyd's IP is an area-related quality measure for Lloyd's MC and thus sensitive to zeros. As can be seen in our example (Table 2), the "spillover" from a crowded to an empty macro-habitat in the area of the decision boundary leads to misleadingly high Lloyd's IP, and to a lesser degree, misleadingly high Bez's IoA, even though this Index is supposedly insensitive to zeros. In accordance with the recommendation by Bez (2000) we therefore suggest Lloyd's MC as a measure of aggregation within a macrohabitat. Lloyd's IP and Bez's IoA might still be informative if the overall colonialization of the macro-habitat is considered.

Pelagic Habitats
The model segregated three (four) distinct pelagic habitats in HE429: (1) a SL (SL1/SL2) mainly characterized by temperatures >17 • C, (2) a BL on the other side of that threshold, and (3) a PL dominated by high chlorophyll concentrations. In contrast to SL and BL, PL was not a true cluster by the definition of HDBSCAN, which indicates a great variability within the microhabitats belonging to PL. That makes them "special" or at least "different" from common micro-habitats in SL and BL. Microhabitats of PL were usually located around the 17 • C isoline and at the occurrence of exceptionally strong chlorophyll-peaks.
In the North Sea, peaks of primary production following the spring bloom were observed in subsurface layers  (Richardson et al., 1998(Richardson et al., , 2000. The PL most likely resembles such areas of subsurface productivity. Furthermore, the model detected an upward doming of the thermocline within an OWF, probably caused by enhanced vertical mixing (Segtnan and Christakos, 2015;Floeter et al., 2017;Schultze et al., 2020). The upward doming and the resulting temperature differences are comparable to those observed within cyclonic eddies (Dong and McWilliams, 2007;Marmorino et al., 2018), indicating that OWF's can influence the pelagic realm in the same order of magnitude as natural (sub-) mesoscale processes like eddies. The doming of colder, nutrient-rich water can produce chlorophyll peaks (Munk et al., 1999), indicating the potential for an enhanced primary production in this area.
Cumulative effects of single foundations might lead to a blocking effect around OWF's, similar as observed for islands (Simpson et al., 1982), which has the potential to produce submesoscale eddies (Dong and McWilliams, 2007) in addition to local upwelling fronts (Floeter et al., 2017). Common properties that are used to describe hydrographic eddies and fronts include water velocity, vorticity and the Rossby number (e.g., Lévy et al., 2012;Marmorino et al., 2018), all of which were not available to us, which makes it less likely to detect such features.
The situation during HE534 was fundamentally different from HE429, most likely due to the weather conditions prior to the cruise that dispersed a thermal stratification. However, the projections indicated that similar SL and BL as in HE429 were detected. In case of HE534, segregations occurred along an oxygen isoline (>235 µmol l −1 ) instead of temperature (>17 • C) as during HE429. This is in accordance with findings of Friedland et al. (2020) and references within that the predictive power of variables might change. This variable nature inherently present in pelagic data (Hinchey et al., 2008;Thompson et al., 2016) makes it so challenging to accurately predict pelagic habitats. The temperature isolines in T1 and T2 (Supplementary Figure 3) clearly indicate the presence of a tidal mixing front (see Hill et al., 1993). A convergence slick, which is typically associated with such tidal mixing fronts (Hill et al., 1993), would also explain the observed aggregation of plankton particles at the intersection of the two macro-habitats (Figure 8).
There exists plenty of evidence that physical properties also structure the marine plankton communities (e.g., Swalethorp et al., 2015;Van Leeuwen et al., 2015;Lindegren et al., 2020). However, the generated habitat maps have only limited explanatory power considering the observed plankton communities. This is not unexpected since physical properties are merely incomplete predictors for the community structure which is most likely further shaped by niche-based processes and interactions (Houliez et al., 2021).
Top predators aggregate in areas with the highest preypatch densities (not to be confused with the area of highest prey densities!) (Benoit-Bird et al., 2013) and peak abundances of zooplankton and fish larvae are frequently observed in the direct vicinity of frontal convergence zones (Munk et al., 1995(Munk et al., , 2002Höffle et al., 2013;Munk, 2014;Swalethorp et al., 2015). In addition to the horizontal agglomerations, thermo-, and haloclines can produce further vertical structuring (Höffle et al., 2013;Lindegren et al., 2020). Thereby, more pronounced differences lead to a stronger niche separation and less overlap between different species (Lindegren et al., 2020). Changes in nitrate (Scharfe and Wiltshire, 2019) and silicate (Wiltshire et al., 2015) availability produce a temporal succession of different dominant taxa in the tidal advected phytoplankton community. Especially the plankton community is thus shaped by complex spatio-temporal dynamics and local prey patches have the potential to shape the distribution of higher trophic levels (Pope et al., 1994;Burkhard et al., 2011;Benoit-Bird et al., 2013;Defriez et al., 2016), even though this might be of less importance for ecosystem services in a highly diverse and partly functionally redundant plankton communities like that of the North Sea (Atkinson et al., 2015).

COMMENTS AND RECOMMENDATIONS
Future work should aim to include species densities and water current related measurements in order to accurately predict not only physical habitats but also realized ecological niches and hopefully improve our understanding of the complex dynamics shaping the pelagic realm.
Our approach offers beneficial properties to solve this challenge: the AE is a highly non-linear tool to reduce the dimensionality of a nearly unlimited amount of data that can be extended as needed. Additionally, HDBSCAN is a cluster algorithm that makes as few assumptions as possible, i.e., regarding number or shape of clusters. HDBSCAN can also handle outliers on it's own in opposite to, e.g., k-means, and even enables to treat them in our case as an own macro-habitat. While machine learning might not give insight into the underlying mechanistic, it can give a starting point from which to begin future investigations (Friedland et al., 2020).

DATA AVAILABILITY STATEMENT
The raw data supporting the conclusions of this article will be made available by the authors, without undue reservation.

AUTHOR CONTRIBUTIONS
SV and R-MP constructed the Autoencoder used to process the data. R-MP was responsible for the data analyses. JF and R-MP wrote the submitted manuscript. All authors contributed to the article and approved the submitted version.

FUNDING
The RV Heincke Research Cruises were supported by grant numbers AWI_HE429_00 and AWI_HE534_00.