Discrimination of Spatial Distribution of Aquatic Organisms in a Coastal Ecosystem Using eDNA

Featured Application: Complement monitoring tool. Abstract: The nonlinearity and complexity of coastal ecosystems often cause di ﬃ culties when analyzing spatial and temporal patterns of ecological traits. Environmental DNA (eDNA) monitoring has provided an alternative to overcoming the aforementioned issues associated with classical monitoring. We determined aquatic community taxonomic composition using eDNA based on a meta-barcoding approach that characterizes the general ecological features in the Gwangyang Bay coastal ecosystem. We selected the V9 region of the 18S rDNA gene (18S V9), primarily because of its broad range among eukaryotes. Our results produced more detailed spatial patterns in the study area previously categorized (inner bay, main channel of the bay and outer bay) by Kim et al. (2019). Speciﬁcally, the outer bay zone was clearly identiﬁed by CCA using genus-level identiﬁcation of aquatic organisms based on meta-barcoding data. We also found signiﬁcant relationships between environmental factors. Therefore, eDNA monitoring based on meta-barcoding approach holds great potential as a complemental monitoring tool to identify spatial taxonomic distribution patterns in coastal areas.


Introduction
Biological monitoring contributes to the understanding of complex ecosystem structures and functions in targeted systems [1]. Accordingly, it is crucial in detecting and assessing environmental changes in order to ensure proper management and conservation of complex ecosystems [2]. Coastal environments are among the most complex ecosystems due to tidal activity, and typically retain high economic and environmental values in light of aquatic resources and biodiversity [3]. Coastal ecosystems are often severely affected by anthropogenic activities such as industrial fishing [4,5], marine transport and leisure activities [6], aquaculture and the aquarium trade [7,8], living seafood and lure fisheries [9,10], and non-indigenous species (NIS) which induce greater pressures on endemic ecosystems and often drive native species to extinction in the resulting habitats [11][12][13]. The nonlinearity and complexity of coastal ecosystems due to the aforementioned activities often causes difficulties when analyzing spatial and temporal patterns of ecological traits.
Environmental DNA (eDNA) as described by Ogram et al. [14] who extracted microbial DNA from the sediment, has increasingly been used in recent years for biological monitoring purposes. Recently, a large number of papers have reported the use of eDNA monitoring in analyses of soil, water and even air [15]. Andersen et al. [16] examined the possibility of monitoring large mammals using eDNA in soil samples, and eDNA from water monitoring of fish [17][18][19][20] and amphibians [21,22] has been successful. Furthermore, Hawkins et al. [23] demonstrated that a complete taxonomic list of functional feeding

Study Area
Gwangyang Bay is part of the Korean National Archipelagos located off the south coast of the Korean peninsula ( Figure 1). The bay receives an annual mean discharge of 2298 × 10 6 m 3 yr −1 from the Seomjin River [26]. A significant amount of nutrients drains into the system from the watershed (~5000 km 2 ). The water depth varies from 10 m at the Seomjin River estuary to 50 m at the outer Gwangyang Bay. The tidal cycle appears to be semi-diurnal. Compared to other Korean river estuaries which have barriers, the Seomjin River estuary remains open, and thus the water mass is exchanged between the river and ocean more actively. The natural condition of Gwangyang Bay is apt to increase primary productivity as well as biological diversity. In this respect, Gwangyang Bay (~450 km 2 from the estuary to the outer bay) is the most economically productive coastal ecosystem in the Korean peninsula [25].

Sampling, Data Collection and Primer Selection
A survey permitted by Ha-dong local government (permission number: 2010-0165) was conducted in June 2018. We sampled the surface water (approximately the top 50 cm) in this study. In total, nineteen sampling sites were selected, which covered an extensive area from the Seomjin River estuary to the outer Gwangyang Bay (Figure 1). According to Kim et al. [25], at Gwangyang Bay, higher water temperatures corresponded to lower salinity, and vice versa. Phosphorus and nitrogen concentrations were spatially similar across the Bay. Based on the divisions indicated in Kim et al. [25], Zone I covered the inner Bay (sites 3-7, and 9), Zone II represented the main channel of the Bay (sites 8, 10, and 11), and Zone III mostly belonged to the outer Bay (sites 12-21) (Figure 1). Water samples for meta-barcoding analysis (more than 1 L per sample) were obtained at the same time and moved with dry ice to laboratory to filter (0.45 µm pore-size membrane; Advantec MFS membrane filter, Dublin, USA) and stored at −80 • C before NGS analysis. Negative controls were included for every study site to prevent cross contamination. Water temperature and salinity were measured on-site using portable equipment (Model: YSI Professional Plus, OH, USA), while the nutrient and chlorophyll-a concentrations were analyzed in the lab. Phosphorus and nitrogen concentrations were measured using an UV spectrophotometer based on standard analytical methods proposed by the Korean Ministry of Oceans and Fisheries. Chlorophyll-a measurements were also based on UV spectrophotometry. In contrast to nutrient measurements, chlorophyll-a samples were filtered through a 0.45 µm pore-size membrane (Model: Advantec MFS membrane filter). The filter membrane was then homogenized after acetone extraction prior to spectrophotometry. Organic and inorganic carbon concentrations were measured using a carbon analyzer (Model: vario TOC cub, Langenselbold, Germany) using 850 • C combustion catalytic oxidation methods.
Appl. Sci. 2019, 9, x FOR PEER REVIEW 3 of 12 every study site to prevent cross contamination. Water temperature and salinity were measured onsite using portable equipment (Model: YSI Professional Plus, OH, USA), while the nutrient and chlorophyll-a concentrations were analyzed in the lab. Phosphorus and nitrogen concentrations were measured using an UV spectrophotometer based on standard analytical methods proposed by the Korean Ministry of Oceans and Fisheries. Chlorophyll-a measurements were also based on UV spectrophotometry. In contrast to nutrient measurements, chlorophyll-a samples were filtered through a 0.45 μm pore-size membrane (Model: Advantec MFS membrane filter). The filter membrane was then homogenized after acetone extraction prior to spectrophotometry. Organic and inorganic carbon concentrations were measured using a carbon analyzer (Model: vario TOC cub, Langenselbold, Germany) using 850 °C combustion catalytic oxidation methods. We selected the V9 region of the 18S rDNA gene (18S V9), primarily because of its broad range among eukaryotes [27,28]. NGS approaches using the 18S V9 region have recently allowed the characterization of marine planktonic biodiversity in the oceans [29] and prompted biomarker establishment initiatives [30].

DNA Extraction and Metagenomic Sequencing
Genomic DNA was extracted using PowerSoil ® DNA Isolation Kit (Cat. No. 12888, Qiagen, Düsseldorf, Germany) according to the manufacturer's protocol. DNA extracted for sequencing was prepared according to Illumina 18S Metagenomic Sequencing Library protocols (San Diego, CA, USA). DNA quantity, quality, and integrity were measured by PicoGreen (Thermo Fisher Scientific, Waltham, MA, USA) and VICTOR Nivo Multimode Microplate Reader (PerkinElmer, Akron, OH, USA). The 18S rRNA gene was amplified using 18S V9 primers. The primer sequences are as follows: 18S V9 primer including adaptor sequence (Forward Primer: 5′ TCGTCGGCAGCGTCAGATGTGTATAA GAGACAGCCCTGCCHTTTGTACACAC 3′/Reverse Primer: 5′ GTCTCGTGGGCTCGGAGATG TGTATAAGAGACAGCCTTCYGCAGGTTCACCTAC 3′). To amplify the target region attached with adapters, as a first PCR process, the extracted DNA was amplified by 18S V9 primers with one cycle of 3 min at 95 °C, 25 cycles of 30 s at 95 °C, 30 s at 55 °C, 30 s at 72 °C, and a final step of 5 min at 72 °C for amplicon PCR product. As a second process, to produce indexing PCR, the first PCR product was subsequently amplified with one cycle of 3 min at We selected the V9 region of the 18S rDNA gene (18S V9), primarily because of its broad range among eukaryotes [27,28]. NGS approaches using the 18S V9 region have recently allowed the characterization of marine planktonic biodiversity in the oceans [29] and prompted biomarker establishment initiatives [30].

DNA Extraction and Metagenomic Sequencing
Genomic DNA was extracted using PowerSoil ® DNA Isolation Kit (Cat. No. 12888, Qiagen, Düsseldorf, Germany) according to the manufacturer's protocol. DNA extracted for sequencing was prepared according to Illumina 18S Metagenomic Sequencing Library protocols (San Diego, CA, USA). DNA quantity, quality, and integrity were measured by PicoGreen (Thermo Fisher Scientific, Waltham, MA, USA) and VICTOR Nivo Multimode Microplate Reader (PerkinElmer, Akron, OH, USA). The 18S rRNA gene was amplified using 18S V9 primers. The primer sequences are as follows: 18S V9 primer including adaptor sequence (Forward Primer: 5 TCGTCGGCAGCGTCAGATGTGTATAA GAGACAGCCCTGCCHTTTGTACACAC 3 /Reverse Primer: 5 GTCTCGTGGGCTCGGAGATG TGTATAAGAGACAGCCTTCYGCAGGTTCACCTAC 3 ). To amplify the target region attached with adapters, as a first PCR process, the extracted DNA was amplified by 18S V9 primers with one cycle of 3 min at 95 • C, 25 cycles of 30 s at 95 • C, 30 s at 55 • C, 30 s at 72 • C, and a final step of 5 min at 72 • C for amplicon PCR product. As a second process, to produce indexing PCR, the first PCR product was subsequently amplified with one cycle of 3 min at 95 • C, eight cycles of 30 s at 95 • C, 30 s at 55 • C, 30 s at 72 • C, and a final step of 5 min at 72 • C. The final products were normalized and pooled using PicoGreen (Thermo Fisher Scientific, Waltham, MA, USA), and the size of libraries were verified using the LabChip GX HT DNA High Sensitivity Kit (PerkinElmer, Akron, OH, USA).
The library was sequenced using the MiSeq™ NGS platform (Illumina, San Diego, CA, USA) provided as a commercial service (Macrogen Inc., Seoul, Korea). Raw reads were trimmed with CD-HIT-OTU and chimeras were identified and removed using rDnaTools. For paired-end merging, FLASH (Fast Length Adjustment of SHort reads) version 1.2.11 was used. Merged reads were processed using Qiime version 1.9 [31] and were clustered into operational taxonomic units (OTUs) with UCLUST [32], using a greedy algorithm with OTUs at a 97% OUT cutoff value. Taxonomic classifications were assigned to the obtained representative sequences using BLASTn [33] and UCLUST [32].

Data Analysis and Statistics
OTUs assigned by meta-barcoding were classified into habitat types (marine, freshwater and estuarine), trophic level (first consumer [filter feeder], second consumer [carnivore], producer and symbiotic) and ISR (indigenous, non-indigenous). Canonical correspondence analysis (CCA) of OTUs exhibiting >1% relative abundance was performed using the PAST 3.0 program [34] to evaluate relationships between abundance of OTU sequences (based on genus level identification) and environmental factors (temperature (Temp.), salinity, total phosphorus (TP), total nitrogen (TN), chlorophyll-a (Chl-a), total organic carbon (TOC), total inorganic carbon (TIC), total carbon (TO), elemental carbon (EC)). All CCA results were constructed using relative abundance data, with natural logarithms transformation (ln 1 + X) used for sample normalization. Linear relationships (Pearson correlations) were calculated between the above environmental factors and classified OTU sequences, based on habitat using XLSTAT version 2018.6.54467 (64 bit) as a plug-in for the Microsoft Excel program [35].

Meta-Barcoding
In total, 3,066,013 paired-end reads from the 19 samples were generated on the Illumina MiSeq™ platform, of which 98.5% passed Q30 (Phred quality score > 30) for improving accuracy of sequences in this study (Supplementary Materials). Each sample yielded paired-end reads ranging from 21,101-299,305 reads (mean: 161,369 reads), which was similar to the amount of reads in the previous study [36], and all samples exhibited saturation of the number of OTUs by rarefaction curve analysis. Gamma-diversity was 352 OTUs produced with a cutoff of 97% similarity. The resulting 352 OTUs were classified into 19 genus-level taxonomic groups (those representing <0.04% abundance were not plotted). Uncultured and non-assigned reads were discarded.

Spatial Distributions of Aquatic Organisms Based on Meta-Barcoding
We carried out a survey in June 2018. Relationships among the environmental variables and aquatic organisms based on meta-barcoding were explored by CCA (Figure 2). The first axis explained 33.3% of the variance and distinguished Zones III-2 and III-3 with higher salinity, EC, TIC, and lower Chl-a and TN from the other zones. The second axis (19.6% of variance explained) distinguished the sites from Zone III-1 and the other zones, by having higher TOC and TC, lower temperature, and TP. Nutrient-related factors (Chl-a, TN, and TP) were correlated with Acropora sp. (small polyp stony coral) and Megabalanus sp., (barnacle) and salinity was correlated with Skeletonema sp. (diatoms) and Centropages sp. (copepods).  Abundance of assigned OTU sequences showed different patterns among the study sites ( Figure 3). The dominant OTU was assigned to Acropora sp. and the sub-dominant OTU was Acartia sp. (copepods). Most common OTUs within the study sites were comprised of phytoplankton, followed by zooplankton and amphipods ( Figure 3B). Interestingly, the abundance of OTUs showed different compositions among the three different zones (I, II, III). Zone I showed similar patterns to Zones II and III-1, whereas Zones III-2 and 3 showed different compositions of aquatic organisms which distinguished them from other zones. In particular, Zone III-3 exhibited an entirely different pattern, comprised of marine organisms such as Euchaeta sp. (copepods) and Tessarabrachion sp. (krill) compared with other zone divisions ( Figure 3B).
Appl. Sci. 2019, 9, x FOR PEER REVIEW 6 of 12 Abundance of assigned OTU sequences showed different patterns among the study sites ( Figure  3). The dominant OTU was assigned to Acropora sp. and the sub-dominant OTU was Acartia sp. (copepods). Most common OTUs within the study sites were comprised of phytoplankton, followed by zooplankton and amphipods ( Figure 3B). Interestingly, the abundance of OTUs showed different compositions among the three different zones (I, II, III). Zone I showed similar patterns to Zones II and III-1, whereas Zones III-2 and 3 showed different compositions of aquatic organisms which distinguished them from other zones. In particular, Zone III-3 exhibited an entirely different pattern, comprised of marine organisms such as Euchaeta sp. (copepods) and Tessarabrachion sp. (krill) compared with other zone divisions ( Figure 3B).

Functional Features and Non-Indigenous Species
When divided into categories (habitat types, trophic level, and ISR), averages of abundant sequences showed different patterns among the zones (I, II, III-1-2-3). The most dominant habitat type of OTU was the marine type for the three different zones ( Figure 4A), while the types of feeding habit and indigenous rate showed complex response patterns. The dominant trophic level was first consumer and second consumer, followed by producer and symbiotic trophic levels ( Figure 4B). The indigenous rate of Zone III-3 was significantly higher than other zones, but Zones II and III-1 showed opposite patterns ( Figure 4C), while Zones I and III-2 presented no significant difference in indigenous rates. The three divisions of Zone III from CCA results demonstrated similar patterns within the three zone types (Zones I, II-and III).
We found significant relationships between environmental factors (Table 1). Specifically, temperature had positive relationships with nutrient-related factors such as TP and PN, and showed a negative relationship with TOC. Salinity showed significant negative relationships with TN and TP, but also displayed a positive relationship with carbon-related factors such as EC and TC. However, there were no significant relationships between environmental factors and divisions of categories except for habitat types (marine, estuarine, freshwater). The marine habitat of OTUs revealed a significant positive relationship with salinity (r = 0.879) and EC (r = 0.878), respectively, whereas freshwater and estuarine types of OTUs showed negative relationships with salinity (r = −0.888 and −0.856) and EC (r = −0.884 and −0.856) respectively.

Functional Features and Non-Indigenous Species
When divided into categories (habitat types, trophic level, and ISR), averages of abundant sequences showed different patterns among the zones (I, II, III-1-2-3). The most dominant habitat type of OTU was the marine type for the three different zones ( Figure 4A), while the types of feeding habit and indigenous rate showed complex response patterns. The dominant trophic level was first consumer and second consumer, followed by producer and symbiotic trophic levels ( Figure 4B). The indigenous rate of Zone III-3 was significantly higher than other zones, but Zones II and III-1 showed opposite patterns ( Figure 4C), while Zones I and III-2 presented no significant difference in indigenous rates. The three divisions of Zone III from CCA results demonstrated similar patterns within the three zone types (Zones I, II-and III).
We found significant relationships between environmental factors (Table 1). Specifically, temperature had positive relationships with nutrient-related factors such as TP and PN, and showed a negative relationship with TOC. Salinity showed significant negative relationships with TN and TP, but also displayed a positive relationship with carbon-related factors such as EC and TC. However, there were no significant relationships between environmental factors and divisions of categories except for habitat types (marine, estuarine, freshwater). The marine habitat of OTUs revealed a significant positive relationship with salinity (r = 0.879) and EC (r = 0.878), respectively, whereas freshwater and estuarine types of OTUs showed negative relationships with salinity (r = −0.888 and −0.856) and EC (r = −0.884 and −0.856) respectively. Appl. Sci. 2019, 9, x FOR PEER REVIEW 7 of 12    Temp.), Salinity, total phosphorus (TP), total nitrogen (TN), chlorophyll-a (Chl-a), total organic carbon (TOC), total inorganic carbon (TIC), total carbon (TO), elemental carbon (EC)) and classified OTU sequences based on habitat types, feeding habits and natives (Bold: significant relationship).

Effectiveness of eDNA Monitoring
Our results showed that eDNA monitoring based on NGS holds great potential as a complementary monitoring tool to identify spatial taxonomic distribution patterns in coastal areas. We characterized detailed zonation patterns of the outer bay of Gwangyang Bay previously categorized by Kim et al. [25]. Specifically, Zones III-2 and III-3 were clearly resolved by CCA using genus-level identification of aquatic organisms based on meta-barcoding data ( Figure 2). When divided into categories (habitat types, trophic level, and ISR), averages of abundant sequences showed different patterns among the zones (I, II, III-1-2-3). The NIS detection of Zone III-3 was significantly higher than the other zones. Therefore, eDNA based on meta-barcoding is a promising ecological tool for monitoring, such as biodiversity assessment or NIS management [24,[37][38][39][40].
This approach by meta-barcoding in complex coastal ecosystems confers three distinct advantages over other methods. First, the range of aquatic organisms in coastal ecosystems open to study is widened. Second, identification of aquatic organism to the species or genus level (Appendix A), which can be converted into ecological values such as FFC and NIS, is possible. Finally, sensitivity of the meta-barcoding approach using small volumes of water (1 L) is a promising alternative to traditional methods (i.e., dredges and surber samplers) for biodiversity detection in coastal areas. Using these advantages, eDNA monitoring can provide a useful tool for use in environmental management.
However, some taxa were unable to be identified to the species or genus level due to the incompleteness of reference databases (i.e., NCBI GenBank). More accurate target regions such as the cytochrome oxidase I region, which is called the standard region of barcoding, are needed to identify target species at the species level and are required for identification of aquatic organisms [41]. Another challenge associated with eDNA monitoring is the risk of false-positive and false-negative detections [42]. The reliability of the eDNA sampling method should be demonstrated using in silico, in vitro and in situ validation tests in the coastal area. Even though this approach has limitations, taxonomic expertise is not required and it can supplement observational records and field surveys to obtain marine ecosystem samples and information, as NGS reveals biodiversity and the number of NIS in one specific area along the their temporal and spatial distribution.

Ecological Values of eDNA Monitoring
Our results showed that Acropora sp. was dominant in terms of abundance of OTU sequences across all study sites, even though it is not a planktonic species. The Acropora colonies post-recruitment by larval recruitment demonstrates higher efficiency for survival than coral-colony growth [43]. It was concluded that larval recruitment largely determines species composition, and that reduced larval recruitment is responsible for the sparse distribution of fragmenting species [44]. Therefore, the larval stage could be easily detected among the study sites by our eDNA monitoring approach. We also found that nutrient-related factors (Chl-a, TN, and TP) were significantly correlated with Acropora sp. and Megabalanus sp. (Figure 2). It is therefore possible that the strong positive correlation between the two species could be the result of coral-inhabiting barnacles [45].
The comprehensive understanding of feeding characteristics of aquatic organisms has not been well elucidated in comparison to their importance [46][47][48][49]. Our results showed that composition of aquatic organisms have different patterns of feeding habits among the three different zone divisions designated by Kim et al. [25] (Figure 4), and we also found similar patterns using CCA based on genus-level identification resulting from our division of Zone III ( Figure 2). However, previous research findings were attributed to the absence of adequate information for analysis obtained from the field sites due to difficulties in culture, handling, and identification of eDNA samples derived from bulk sediment and filtered water. In this sense, our results overcame the aforementioned limitations by using a broader detection range for aquatic organisms.
Although, Zhan et al. [50] has described the increased sensitivity of meta-barcoding for NIS, its application for monitoring biological invasion has only recently been demonstrated [51]. In the present work, we identified the utility of meta-barcoding for detection of NIS and their spatial distribution patterns ( Figure 4). The reason behind this is the capacity to detect the presence of individuals at early life stages, such as eggs or nauplius larvae, whose identification is difficult with traditional methods [52]. The sensitivity of meta-barcoding, combined with the relatively low time and cost associated with this technique [53], makes it a promising alternative approach for the rapid and accurate detection of biodiversity shifts in aquatic organisms, allowing its potential implementation in environmental policies.