*3.1. 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.

### *3.2. 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-*<sup>a</sup>*, 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).

**Figure 2.** Canonical correspondence analysis (CCA) used to evaluate (**A**) the relationships between abundance of operational taxonomic units (OTUs) 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)). (**B**) Sites grouping based on CCA analysis.

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

**Figure 3.** Abundance of OTU sequences along the study sites (**A**) with site description of physical features and zone classification (**B**) based on genus identification level.

### *3.3. 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.

**Figure 4.** Comparison of the habitat types ((**A**),marine-freshwater-estuarine), trophiclevel ((**B**), consumer 1 [filter feeder], consumer 2 [carnivore], producer and symbiotic) and natives ((**C**), non-indigenous and indigenous) among five zoning types based on genus identification level.


**Table 1.** Linear relationship (Pearson correlation) between 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)) and classified OTU sequences based on habitat types, feeding habits and natives (Bold: significant relationship).
