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).