Methods
Oceanographic data
Monthly CTD profiles of temperature, salinity, density, oxygen, and
chlorophyll fluorescence paired with discrete bottle sampling, provided
hydrographic context. The details of sampling, quality control and
post-processing are available on the BATS website
(http://bats.bios.edu/index.html). Each profile was assigned to
one of four seasons (mixed, spring, stratified, fall; Table 1) based on
retrospective evaluation of the full time series for each year. Vertical
layer boundaries (Table 1) were determined on a profile-by-profile basis
from the fluorometer and density criteria.
Biological samples were acquired monthly from February 2016 to December
2018. March 2016 was sampled twice (March 8th and
25th); March 2017 sampling was cancelled due to rough
weather. The protocol included 12 fixed depths between surface and 1000
m. At each depth, 4 L of water were filtered using 0.22µM Sterivex
filters (MERCK, MA, USA). Filtering pressure was kept at 5 “ Hg. Both
ends were sealed with parafilm or a sterile plug/cap. Filters were
wrapped in foil and stored at -80°C.
Molecular methods
Sucrose lysis buffer (1mL) was added to the Sterivex filters upon
thawing. DNA was extracted using the method by Giovannoni et al. (1990).
Amplification and sequencing was carried out in the Center for Genome
Research and Biocomputing at Oregon State University. The 18S rRNA V4
hypervariable region was amplified using the primers
5-CCAGCA[GC]C[CT]GCGGTAATTCC-3 and
5-ACTTTCGTTCTTGAT[CT][AG]A-3 (Stoeck et al. 2010), using
the KAPA HiFi HotStart ReadyMix (Kapa Biosystems; Wilmington, MA, USA).
The PCR thermal cycler program consisted of a 98°C denaturation step for
30 s, followed by 10 cycles of 10 s at 98°C, 30 s at 53°C and 30 s at
72°C, and 15 cycles of 10 s at 98°C, 30 s at 48°C and 30 s at 72°C, and
a final elongation step at 72°C for 10 min. After PCR, dual indices and
Illumina sequencing adapters were attached using the Illumina Nextera XT
Index Kit. Amplicon sequencing was conducted on an Illumina MiSeq using
2x250 PE V2 chemistry.
Bioinformatics and analyses
Quality control and read merging were carried out with MeFiT (Parikhet al. 2016), using Jellyfish Ver. 2 (Marçais & Kingsford 2011),
with strict parameters to compensate for the lack of complete overlap
between forward and reverse: CASPER (Kwon et al. 2014)K -mer size of 5, threshold for difference of quality score 5,
threshold for mismatching ratio 0.1, and minimum length of overlap of
15. The Merge and Quality Filter Parameters (using the meep score
(Koparde et al. 2017)) were set at a filtering threshold of 0.1,
discarding non-overlapping reads. Assembled reads passing the QC were
fed into MOTHUR ver. 1.43 (Schloss et al. 2009). In short, after
trimming to the V4 region (by aligning to the SILVA 138 release trimmed
to the V4 region) and removal of chimeras with VSEARCH (Rognes et
al. 2016), single variants were obtained using UNOISE2 (Edgar 2016) as
implemented in MOTHUR (diffs=1). Taxonomic assignment was done using the
PR2 database (Guillou et al. 2013). After
excluding metazoans and “Eukaryota_unclassified”, samples were
rarefied to 20,000 reads and all ASVs with less than two reads (global
singletons) were equaled to zero. Scripts are available at the GitHub
repository https://github.com/blancobercial/Protist_Time_Series.
Counts per sample and environmental metadata were loaded into PRIMER
Ver. 7 (Clarke & Gorley 2015). Diversity indices (S, H’, J’) were
calculated in PRIMER; Faith’s phylogenetic diversity PD (Faith 1992) was
calculated in MOTHUR using the same input table as for the other
indices. Samples were then standardize by the total, and square-root
transformed (Hellinger transformation) (Legendre & Gallagher 2001). A
Bray-Curtis distance similarity matrix was built, and a Principal
Coordinates Analysis (PCoA; Gower 1966) was carried out. The
relationship between environmental variables and the similarity between
samples was analyzed with BioEnv in PRIMER, using Spearman Rank and up
to five variables for the model. The variables included day of year,
depth, temperature, salinity, oxygen, fluorescence, turbidity and
density. A second analysis included the factor “hydrographic layer”
(Table 1), to assess how this partitioning compared to the protist
community. All environmental variables were normalized (subtracting the
mean and dividing by the standard deviation) before analyses to remove
biases associated to different variable scales.
Non-hierarchical clustering based on group-average ranks was used to
find the best grouping reflecting the self-organization of the
community, from number of clusters K =2 until two of the clusters
were not statistically different (K =14). At each K level,
an analysis of similarities (ANOSIM) was run (with 999 permutations) to
assess the R -value at each clustering level. The K levels
showing the maximum R -values were further analyzed. Clustering
was mapped to the hydrographic layers and evaluated within the PCoA
ordination, and their taxonomic composition at PR2“Clade” level analyzed.
The molecular lineages (single variants) with the largest contribution
to the similarity in each cluster, and to dissimilarity between adjacent
clusters, were assessed using a Similarity Percentage (SIMPER) analyses.
Fine taxonomic assignment was done using BLAST(Altschul et al.1990), analyzing several of the top scores to understand the precision
of assignments (species/genus/family/etc.). Interpretation the changes
was contextualize within the oceanographic landscape and transitions
between clusters, and with the functional profile (based in the
literature) of each lineage.
PCoA analyses were further performed at each discrete depth to study the
interaction between seasonality and depth. Same procedure was repeated
following those samples that overlapped with the Chl-a maximum peak for
each profile (which did not happen every month because a set of fixed
depths were sampled each time). SIMPER analyses were run to detect taxa
related to seasonal succession at each depth.
For all analyses, a Holm-Bonferroni correction was applied to the
significance threshold when multiple comparisons were done (Holm 1979).