Biodiversity patterns of epipelagic copepods in the South Pacific Ocean: Strengths and limitations of current data bases

communities to determine biogeographical regions, as a community that is a subset of a larger one (nestedness) has no exclusive species, and thus such difference should not be taken into account to delimit regions; whereas by using indexes not affected by nestedness it is possible to delimit regions with unique biological communities [58]. We therefore calculated alpha diver- sity as species richness and beta diversity as species composition while distinguishing between turnover and nestedness. Alpha diversity was calculated in Biodiverse 3.1 software [60], whereas beta diversity was estimated using the packages betapart [61], CommEcol [62] and letsR [63] in R software (R Core Team, 2021), and using the following equation: bsor ¼ bsim þ bsne � b þ c 2a þ b þ c ¼ b b þ a þ c b 2a þ b þ c � � a b þ a � � ; ð1Þ where βsor is Sørensen dissimilarity, βsim is Simpson dissimilarity (i.e., turnover component of Sørensen dissimilarity), βsne is the nestedness component of Sørensen dissimilarity, a is the number of shared species between two cells, b the number of species unique to the poorest site, and c the number of species unique to the richest site. As biodiversity index, we used the Shannon-Wiener index which was calculated in Biodi- verse 3.1 software [62], according to Eq 2 [64] in the following way: H ¼ Xn i¼1 pi � ln pi ð2Þ where pi is the number of samples (in this case, occurrences) of the i species as a proportion of the total number of occurrences in the neighborhoods (1˚ x 1˚ sampling cells). This proportion is estimated as: pi ¼ ni N ð3Þ

where ni is the number of records of the i species and N is the total number of records at spe- cies level in the sampled cell. In order to account for the low sampling coverage and explore whether cells with low values of richness were causing biases in the analyses, we discarded the cells with values of alpha diversity < 10, and later performed the spatial analysis using an inverse distance weighted (IDW) technique that interpolates a raster surface from points, limiting to the range of the val- ues used to interpolate and giving more (or less) weight where there is more (or less) data. In addition, to explore further at higher taxonomic levels, we performed the same analysis con- sidering genus and family levels. Both analyses were plotted through Kriging interpolation using ArcMap 10.4.1 software (ESRI 2016). Since species richness in a sample is highly dependent on sample size or sampling efforts [65], and the ocean has under-sampled areas due to spatial biases in the distribution of sam- pling locations [6], a skew on the apparent patterns of marine species richness can be pro- duced. Hence, redundance statistics were used to count the total number of species, and an assessment of sample completeness was done for species richness through the R package iNEXT, which provides functions to compute and plot the seamless rarefaction and extrapola- tion sampling curves [65]. This latter analysis was performed for each main surface current system of the SPO, namely: Cape Horn current, East Australian Current, Humboldt Current

System, South Equatorial Current, South Pacific Subtropical Gyre, and West Wind Drift, to evaluate whether some areas had a greater or lower coverage than others. The previously men- tioned main surface current systems were downloaded from Harvard Dataverse [66]. We generated a biological-environmental spatial database, wherein each cell had the bio- geographic information of copepod alpha and beta diversity together with oceanographic data PLOS ONE Biodiversity of copepods in the South Pacific PLOS ONE | https://doi.org/10.1371/journal.pone.0306440 July 11, 2024 6 / 28