Structure and membership of gut microbial communities in multiple fish cryptic species under potential migratory effects

The animal gut microbiota evolves quickly towards a complex community and plays crucial roles in its host’s health and development. Factors such as host genetics and environmental changes are regarded as important for controlling the dynamics of animal gut microbiota. Migratory animals are an important group for studying how these factors influence gut microbiota because they experience strong environmental perturbations during migration. The commercially important grey mullet, Mugil cephalus, is a cosmopolitan species complex that display reproductive migration behaviour. There are three cryptic species of M. cephalus fish distributed across the Northwest Pacific, and their spawning sites overlap in the Taiwan Strait. This extraordinary natural occurrence makes the grey mullet an ideal model organism for exploring the nature of wild animal-gut microbiota relationships and interactions. This study investigates the diversity and structure of the gut microbial community in three cryptic M. cephalus species using 16S rRNA amplicon sequencing. Gut microbial compositions from adult and juvenile fish samples were analysed. Our results indicate that gut microbial communities within the grey mullet share a core microbiome dominated by Proteobacteria, Firmicutes and Actinobacteria. However, the structures of gut microbial communities were more distinct between adult mullet groups than they were between juvenile ones. Intriguingly, we found that adult fish that migrate to different geographical tracts harbour gut microbiota similar to historical records of seawater microflora, along their respective migration routes. This observation provides new insights into the interaction between aquatic animal gut microbial communities and the environments along their hosts’ migratory routes, and thus warrants future study.

The animal gastrointestinal tract harbours vast and dynamic populations of microbes known as intestinal microbiota. Intestinal microbiota plays an important role in host development and immunity; its structure and some of its special functions can be shaped by multiple endogenous factors in the host and exogenous influences from the environment 1,2 . It has been suggested that the symbiotic microbial community is largely determined by its host's genetic background, and gut microbial compositions between hosts of the same species or family are more similar than between those of different host species or families 3 . Growing evidence suggests that the hosts' genetics influence their gut microbial compositions, which can maintain the stability of the gut microbial community membership throughout the life of an organism, or change as a result of strong perturbations 4 . Moreover, this supports the hypothesis that gut microbial composition plasticity, especially in special external environmental conditions, facilitates host acclimation and adaptation 5,6 . Yet, even though we have a considerable understanding of the interaction between animals and their gut microbiota, the extent to which the environment and host genetics contribute to shaping the microbiota across the host's life remains unaccounted for.
Animal migration occurs seasonally-mainly due to food availability, or for breeding purposes-and has been observed in various groups including birds, mammals, reptiles, amphibians, insects, and fish 7 . Accordingly, migration is considered an ideal event for investigating not only the long-term influence of environmental conditions and food sources in an organism's original habitat, but also the temporal impacts of external factors and environments and adjacent areas in the East China Sea (ECS), whereas NWP2 and NWP3 migrate following the Kuroshio Currents from the SCS in the middle of the winter. Consistent with other migratory fish species, grey mullet adults are usually observed to have empty stomachs, which suggests that the grey mullet also does not consume food during migration. Young grey mullets of all three cryptic species develop in the waters around Taiwan until they reach the appropriate size to migrate to their ideal nursing areas. The grey mullet M. cephalus has three known cryptic species that use the Taiwan Strait for spawning and have fascinating migratory histories, making them ideal model species for investigating whether the divergence in gut microbiota is triggered by the influence of host genetic background and migration.
Here, we studied the gut microbiota of both juveniles and adults of three cryptic species in the Mugil cephalus species complex that spawned in Taiwanese waters to investigate the importance of host genetics and migration history. We hypothesized that migration impacts the gut microbiota. We predicted that adult fishes that migrated from different areas in the Northwest Pacific into Taiwanese waters would have distinct gut microbiota but their juveniles that spawned in Taiwanese waters would share similar gut microbiota. Alternatively, if host genetics were an important factor, then all fish from the species complex-including juveniles and adults-would share the same core microbiota that are resistant to perturbations from migration histories. In addition, we compared the gut microbial community of the grey mullet with historical records of the microbiomes of seawater along their suggested geographical migration tracts in the Northwest Pacific to explore the potential impact of seawater on the gut microbiota of migratory fish.

Results
General features of collected samples and characteristics of sequenced data. We analysed a total 16 adult samples (NWP1 n = 8; NWP2 n = 4; NWP3 n = 4) and 12 juvenile samples (NWP1 n = 6; NWP2 n = 6) from three cryptic species (NWP1, NWP2, and NWP3) in the grey mullet M. cephalus species complex in the Taiwan Strait during their spawning season. Fish samples were identified as specific cryptic species based on methods from previous studies 15,16 (see Materials and Methods). The diversity and structure of the gut microbial composition were determined by performing deep sequencing of bacterial 16S rRNA genes from all fish gut samples. After quality filtering and removing chimeras and single reads, a total of 3,704,851 high-quality reads were obtained from 16 adult (hereafter AD) and 12 juvenile (JV) grey mullet samples. Evaluation of the microbial complexity in the M. cephalus species complex. In this study, 750,613 reads were rarefied in all samples. Rarefaction curves approached the saturation plateau and adequately represented the gut microbial community diversity in all samples ( Supplementary Fig. 1). Although there was no difference between the results of the analyses using rarefied and unrarefied data, it is believed that rarefied data can ignore the presence of rare species, which leads to false positives 17 . Therefore, we prioritized the results from the unrarefied data. At the species level, wild AD grey mullets had an average of 123 ± 9, 177 ± 39, and 166 ± 57 OTUs in NWP1, NWP2, and NWP3, respectively. Meanwhile, we detected 566 ± 64 OTUs in NWP1 and 489 ± 95 OTUs in NWP2 juveniles. Inter-specific variation among individuals in sampled species were examined using a permutational test for homogeneity of multivariate dispersions. Although the box plot showed higher inter-variation between individual samples of NWP3 juvenile adult ( Supplementary Fig. 2), the p value was not significant (p = 0.292). The alpha diversity-including observed richness (sob), Chao diversity, Shannon index (H') and inverse Simpson diversity (S) indices-were calculated. In addition, permutation T-test was conducted to test whether there was any difference in alpha diversity among adult groups, among juvenile groups, and between adults and juveniles of the same cryptic species. The results showed that, while there are similarities between the alpha diversities of AD_NWP1 and AD_NWP3 and between JV_NWP1 and JV_NWP2 in both indices, there was significantly more diversity in juvenile samples compared to adults of the same cryptic species (AD_NWP1 -JV_NWP1: P sob < 0.01, P chao < 0.01, P H' < 0.01, P S < 0.01; AD_NWP2 -JV_NWP2: P sob < 0.05, P chao < 0.05, P S < 0.05). At the same time, AD_NWP2 had the lowest alpha diversity (except in the Chao estimation) compared to AD_NWP1 in the observed richness and Shannon and Inverse Simpson indices and AD_ NWP3 in the Inverse Simpson index (Fig. 2, Supplementary Table 2). Gut microbial composition. The structure of the gastrointestinal microbiota was characterized using the relative abundance of representative sequences from each OTU from each sample to assign a taxonomic classification against the Greengenes database. A total of 45 different cultured and candidate bacterial phyla were detected from all fish gut samples. In adult samples, M. cephalus gut microbiota across the phylum were dominated by Proteobacteria (65.4-89.5%) in cryptic species 1 (NWP1) while the most abundant phylum of bacteria was Actinobacteria (74-94.5%) in cryptic species 2 (NWP2) and Spirochaetes (14.8-66.6%) in cryptic species 3 (NWP3). Furthermore, the gut microbiota of juvenile samples NWP1 were dominated by Proteobacteria (37.9%), Cyanobacteria (33.9%), and Firmicutes (12.5%), while the NWP2 juvenile gut microbial community consisted of high abundances of Firmicutes (46.4%), Proteobacteria (32.3%), and Actinobacteria (13.7%) (Fig. 3A).
At the family level, Moraxellaceae was the most abundant bacterial group in NWP1 adults (accounting from 43.7-77.7%), while one family in the order Actinomycetales dominated 70.2-98.9% of the gut microbial community in NWP2 adults. Samples belonging to cryptic species NWP3 had high abundances of Brevinemataceae and Mycoplasmataceae bacteria (14.9-92.3%). In juvenile samples, Streptococcaceae, Clostriaceae, Vibrionaceae, and one unclassified family belonging to order Stramenopiles were the most representative of the gut microbial compositions, contributing to 30-94.5% of the total abundance ( Supplementary Fig. 3).
In contrast to the relative abundances of gut bacteria, our analysis detected similar proportions of community membership among all grey mullet fish, including adult and juvenile samples. Community membership is calculated as the number of bacteria taxa at the phylum level. Specifically, Proteobacteria, Firmicutes and Actinobacteria accounted for an average of 44, 19, and 18% of identified OTUs in all fish samples, respectively (Fig. 3B).
We next used beta diversity analysis to approximate the variation in gut microbial community structure inside the M. cephalus species complex. Non-metric multidimensional scaling (NMDS) plots showed that grey mullet gut microbial compositions were distinctively separate from different cryptic species in adult samples, and were strongly clustered between juvenile fishes based on Bray Curtis and weighted UniFrac distance matrixes (Fig. 4). PERMANOVA analyses also supported our delineated results (Supplementary Table 3). In particular, there were differences in microbial composition among adults of the three cryptic species (p < 0.05). The microbial communities in adults and juveniles of the same cryptic species were also significantly different (p < 0.001). On the other hand, there were no species-specific differences observed between juvenile samples of NWP1 and NWP2. core gut microbiome of M. cephalus. To investigate the existence of a core gut microbiome that is maintained across the M. cephalus life history, we focused on 122 genus-level OTUs shared in all the fish groups in this study. Consistent with membership of gut microbiota in grey mullets, most of the core microbes were Proteobacteria (56 OTUs, 45.9%), Firmicutes (22 OTUs, 18.03%), and Actinobacteria (18 OTUs, 14.75%). The fact that the relative abundance of the core gut microbes can be as high as 95.8% in the communities might reflect the importance of bacteria inside the gut microbial community, as well as their potential functions in and benefits to the hosts. Within these 122 OTUs, 35 were present in over 90% of observed samples. Permutation T-test results indicated that many core microbes inside the gut microbial community in AD_NWP1 were significantly different with not only the other adult cryptic species (NWP2 and NWP3; 22 and 13 OTUs, respectively), but also the juveniles of this cryptic species (31 OTUs). OTU499 (assigned to a genus of Brevinemataceae) was the only OTU with a notably different abundance among the gut microbiota of all three cryptic species in the adult stage (P NWP1-NWP2 < 0.05, P NWP1-NWP3 = 0.01, P NWP2-NWP3 < 0.05) ( Table 1).
clustering of juvenile gut microbial composition. From the 827 OTUs detected at the genus level in the gut microbial communities of 12 juvenile grey mullets (NWP1 n = 6 and NWP2 n = 6), we identified 598 OTUs shared between the two cryptic species, which accounted for 99.75% and 99.87% of the relative abundance in the gut microbial communities of NWP1 and NWP2 juveniles, respectively (  Table 3). High similarities in libraries between juveniles is consistent with the results of the NMDS plots from the Bray-Curtis distant and weighted UniFrac (Fig. 4). In particular, juvenile samples' gut microbial communities clustered together, and the PERMANOVA test found that they were not significantly different. T-test results showed that only 17 of the 598 shared OTUs were significantly different between NWP1 and NWP2 juveniles (Fig. 5B).
Cryptic species-specific microbiome. In addition to the core microbiome at the species complex level, we clarified unique shared OTUs in specific cryptic species. Fourteen OTUs were only observed in NWP1, while six were only found to be shared between adults and juveniles of NWP2. We ran sequences of unique and shared OTUs in two cryptic species through NCBI BLAST, which identified species names with NCBI IDs. Next, we Significant p values (*p < 0.05; **p < 0.01) were obtained using a permutation T-test. X axis shows the group IDs; JV is juvenile, AD is adult, NWP1 is cryptic species 1, NWP2 is cryptic species 2, and NWP3 is cryptic species 3. www.nature.com/scientificreports www.nature.com/scientificreports/  www.nature.com/scientificreports www.nature.com/scientificreports/ generated phylograms between identified OTUs using Maximum Likelihood analysis based on the Kimura 2-parameter model ( Supplementary Fig. 4). Interestingly, even though its gut was dominated by Actinobacteria in the adult stage, four of the six unique OTUs in NWP2 were Proteobacteria (i.e. Neisseria sp., Aquabacterium sp., Rheinheimera sp.and Chujaibacter sp.). At the same time, three Actinobacteria species (Corynebacterium sp., Agrococcus sp. and Janibacter sp.) were only detected in NWP1. Moreover, we detected two closely related Neisseria spp. in the gut microbiota of these two different cryptic species.
Potential impact of seawater microflora on the fish gut microbial community. The influence of host genetics on the gut microbiota was partly reflected in the similarities in communities among juvenile fishes of two closely-related cryptic species, as well as the observation that compositions harboured almost identical membership in all the fish gut samples. However, the clear distinction between structures of gut microbial compositions in adults belonging to three different cryptic species of grey mullet, suggests that exogenous factors may also influences the grey mullet gut microbiota. To test whether migration history is associated with differences in the structure of grey mullet gut microbial composition, we compared our fish gut microbiota with historical seawater microflora of 33 seawater samples from previous studies (Supplementary Table 3 www.nature.com/scientificreports www.nature.com/scientificreports/ seawater data were collected in the East and South China Seas along the migratory tracts of three cryptic species in the grey mullet species complex suggested by Shen et al. 15 . Our analysis at the phylum level found that the gut microbial composition of NWP1 and seawater-microbiota in the ECS had similar compositions; in particular, most were dominated by Proteobacteria, but not Actinobacteria or Firmicutes (Fig. 6). A dendrogram was used to show the Pearson correlation coefficient results, which were calculated to elucidate the correlation between the fish microbiota and seawater microbes. Clear clusters were identified between mullet gut microbial communities of AD_NWP1 samples and many sea water microbial compositions in ECS. In contrast, the gut microbiota of fish belonging to cryptic species NWP2 were correlated vigorously with seawater samples SCS_P1, SCS_P2, SCS_A3, SCS_A4, SCS_S12 and ECS_T24, the bacterial communities of which were all dominated by Actinobacteria (Fig. 6, Supplementary Fig. 5). The dendrogram did not show any close clustering between gut microbial compositions in adult NWP3 and historical seawater data.

Discussion
The gut microbiota has a fundamental influence on host fitness by contributing to host metabolic capabilities, immunity levels and development. The gut microbial community is, in turn, influenced by several factors, such as genetics, developmental stages, environment and food intake 1-3 . Migratory animals experience extreme conditions during transit that may strongly impact gut microbiota, but this remains poorly understood. Here, we characterize the structure and membership of gut microbial communities in three cryptic species (NWP1, NWP2, and NWP3), of the migratory fish species complex Mugil cephalus. Multiple comparisons among adult and juvenile groups of different cryptic species were conducted and showed that, while distinct bacterial communities were observed in the adult fish of different cryptic species, there was still a core microbiome in all fish samples, which was dominated by Proteobacteria, Firmicutes, and Actinobacteria. This demonstrates that, despite the dynamics of bacterial relative abundance, the gut microbial communities of M. cephalus maintain an inflexible proportion of specific bacteria groups inside their gastrointestinal tracts. This finding is consistent with observations in humans 4 and other animals such as crickets 22 . We suggest that the similarities in membership, not the structure, of the gut microbiota in different cryptic species of Mugil cephalus might be highly correlated with their close genetic backgrounds. This supports the hypothesis that fish hosts selectively filter ideal/beneficial microbes from specific exogenous species pools to shape their own gut microbiota so as to keep them healthy and increase their survival ability 23 .
At the family level, the gut microbiota in NWP1 and NPW2 samples were dominated by Moraxellaceae (Phylum Proteobacteria) and a member of a family in order Actinomycetales (Phylum Actinobacteria), respectively ( Supplementary Fig. 3). Moraxellaceae and Actinomycetales are popular bacteria groups found in the gut microbiota of many fish species based on analyses of high throughput sequencing data [24][25][26][27] . Surprisingly, this study found that Brevinemataceae (phylum Spirochaetes) and Mycoplasmataceae (phylum Tenericutes) were the most abundant bacterial families in the gut microbiota of NWP3 adults. A diverse range of Spirochaetes are regularly found in marine sediment, deep within soil and in the digestive tracts of arthropods and several mammals. Spirochaetes is commonly detected at low abundances in fish gut microbial communities; however, Givens 28 reported that one of the three barracuda examined contained a community comprised of 99% Spirochaetes and the three Mahi mahi (common dolphinfish) studied contained 64-98% Spirochaetes. We believe that the gut microbiota of marine fish is extremely diverse due to complex influence from various genetic factors and differences in gut habitats. Clostridiaceae, Streptococcaceae (phylum Firmicutes) and Vibrionaceae (phylum Proteobacteria) were represented in the juvenile grey mullets' gut microbial communities. Vibrionaceae, representing Vibrio bacteria, has been observed as a pathogen group in many animals [29][30][31] . Vibrionaceae, represented by the genera Vibrio and Photobacterium in the mullet intestinal tract, is one of the most important bacteria families in aquaculture, and both its pathogenic and probiotic species are commonly associated with marine animals 32 .
Our comparison also showed that there is significantly higher diversity in juvenile samples than adult ones. This finding contrasts with observations of the gut microbiota in humans and other animals 33 . However, our results are consistent with the conclusion presented by Stephens et al. 34 , who suggested that the gut microbiota in www.nature.com/scientificreports www.nature.com/scientificreports/ zebrafish had significantly lower richness and diversity in adults than younger fish. In addition, it was reported that juveniles of wild Atlantic salmon, Salmo salar, contained a more diverse gut microbiota than those in older life stages 35 . The diversity of food sources of juvenile fish might be one of the reasons for the similarities observed across fish gut microbiota. Juvenile grey mullets (M. cephalus) are believed to use several food sources in natural environments for a range of organisms, from copepods and small zooplankton to detritus particles and microalgae [36][37][38] . We observed high abundances of Cyanobacteria-represented by the family Synechococcaceae and genus Synechococcus-in juvenile gut microbiota, especially in NWP1 juveniles. Synechococcus is distributed widely in ocean ecosystems as one of the most common picoplanktonic marine Cyanobacteria genera 39 . Our findings suggest that the presence of these bacteria in grey mullet juveniles may be part of not only the fish gut microbiota, but also bacterioplankton and be ingested along with fish hosts' food items.
Our finding that gut microbial communities in adult fishes of different cryptic species diverged supports our hypothesis about the impact of migration history on the gut microbiota. Therefore, to test whether seawater microflora along the suggested migratory tracts of the three cryptic species are associated with fish gut microbiota, we compared our fish gut data with historical seawater microbial communities [18][19][20][21] . It is believed that microorganisms in seawater change over different timescales according to various biological and non-biological forces in the surrounding environment 40 . However, seawater microbial community data from long-term time series in several previous projects from (sub)tropical to polar regions (BATS (Bermuda Atlantic time-series study) 41,42 , HOT (Hawaiian Ocean Time-series) 43 , LMO (Linnaeus Microbial Observatory) 44 , and SPOT (San Pedro Ocean Time-series) 45,46 ) suggested that seasonal variation in ocean microbial communities is weakest in sub-tropical regions 40 . Our target species, the grey mullet (Mugil cephalus), is mostly distributed in tropical and subtropical zones (between latitudes 42°N and 42°S) 47 , and all fish samples in this study were observed in the Taiwan Strait (22.73°N, 119.60°E). In addition, it was also suggested that there are "internal feedback mechanisms" inside the microbial community-including mutualism, competition, parasitism and commensalism-that can maintain community's stability over several years 40 . Therefore, although our fish gut microbiota were observed at different times with historical marine water microflora, the association between our grey mullet gut microbial community results and published seawater microbial compositions is reliable. Nevertheless, our analysis was only at the phylum level. Our comparison showed that gut microbial communities in AD_NWP1 were highly correlated with seawater samples from the ECS, whereas compositions in AD_NWP2 were highly correlated with some specific seawater samples near the shore in the SCS, a proposed M. cephalus nursing area. This finding contrasts with results suggested by previous studies that fish gut microbial compositions are rarely influenced by environmental conditions such as water, especially in the mature/adult stage 23,34 . Starvation and osmoregulation in migratory fish could explain this observation. It was proposed that starvation is an opportunity for some specific bacterial groups to increase their competitive ability by utilizing alternative energy sources and then enhancing their relative abundances, as demonstrated with the bacteria inside the intestinal tract of the Asian seabass 48 .
In conclusion, the gut microbiota in different cryptic species of the species complex Mugil cephalus are heavily impacted by host genetics, as shown by the presence and high abundance of a core group of microbes (i.e. all sample groups had similar gut microbiota compositions), but strongly structured through life history such as long-distance migration (i.e. distinct gut microbial community structures in the adult stage of different migratory cryptic species). Although the fish sample size is quite small and samples were not taken at the same time as ocean samples were, the observed association between fish gut microbial community and seawater microflora along its migration tracts provides new insights into the host-gut microbiota interaction together with the potential interaction between aquatic animals and their environment in extreme conditions or during perturbations. Further studies with a broader interspecific comparison are highly recommended to investigate in more detail the effects of this exogenous factor on shaping the structure of gut microbiota in aquatic animals. Lastly, clarification of specific invasion microbes and their functions inside the gut microbial community is needed to determine how gut microbiota support their fish host during long-distance migration.

Materials and Methods
Sample collection. Fertile female (i.e. fish with fully ripe internal egg masses in their ovaries) migratory grey mullets, estimated to be 3-4 years old fish, were collected in Kaohsiung offshore during their spawning season from mid-November 2016 to mid-January 2017. We specifically chose this site to collect specimens of all three grey mullet cryptic species migrating to the Taiwan Strait, following the migratory tracts suggested by Shen et al. 15 . In addition, mullet juveniles were collected from the Danshui estuary in February 2017 (n = 12). Collected fish samples were immediately preserved on dry ice and stored at −20 °C. In the laboratory, the fish surface was sterilized with 70% ethanol and sterile water, then dissected under a hood. The gut contents were kept in 2 mL Eppendorf tubes and preserved at −80 °C until DNA extraction. For adult samples, the gut contents of the last 10 cm of the hindgut were collected for further analysis. For juvenile samples, the whole intestine was used, including both intestinal contents and epithelial associated microbes. Standard length and body weight were also recorded. We could not detect sex by visualization in young grey mullets. All animal procedures were approved by the Academia Sinica Institutional Animal Care and Utilization Committee (AS IACUC) (approved protocol #13-09-576). All experimental methods were performed in accordance with the relevant approved guidelines and regulations.
DNA extraction and 16S rRNA amplicon sequencing. DNA was extracted from the gut contents of adult and juvenile mullet samples using the DNeasy Blood & Tissue Kit (Qiagen, Venlo, Netherlands) according to the manufacturer's protocol. DNA concentration was measured with a NanoDrop ND-1000 spectrophotometer (NanoDrop. Technologies, Wilmington, DE). A universal eubacterial primer pair including the forward primer 27F (5′-AGAGTTTGATCMTGGCTCAG-3′) and the reverse primer 335R (5′-GCTGCCTCCCGAGGAGT-3′) 6,49 combined with an adapter sequence and sample barcodes were used to target specific 16S rRNA gene regions (V1 and V2) by PCR amplification. Extracted DNA samples were diluted to 20 µg/mL directly before PCR processing. The total volume of the PCR mixture was 40 µL, and contained 24 µL 2X Phusion Flash High-Fidelity PCR Master Mix (Finnzymes Oy, Finland), 4 µL of each primer (2 µM), 4 µL template DNA and 4 µL nuclease-free water. The cycling conditions consisted of an initial denaturation at 95 °C for 5 min; followed by 30 cycles of denaturation at 95 °C for 1 min, annealing at 55 °C for 30 s, extension at 72 °C for 1 min; and a final extension at 72 °C for 10 min. The purified amplified PCR products were observed by gel electrophoresis before being extracted with a NucleoSpin Gel and PCR Clean-up Kit (Machrey-Nagel GmbH & Co) and quantified using the Nanodrop-1000 Spectrophotometer (Thermo Scientific, Wilmington, DE, USA). Equal amounts of PCR products were pooled together following the Illumina standard protocol for 16S rRNA sequencing library preparation. Sequencing was performed using the Illumina MiSeq platform (Illumina, San Diego, CA, USA) with the reagent kit v3 at the NGS High Throughput Genomics Core Facility at Academia Sinica. Because collected samples were randomised and high quality samples were refined, we observed and analysed NGS data from a total 16 adult samples (NWP1 n = 8; NWP2 n = 4; NWP3 n = 4) and 12 juvenile samples (NWP1 n = 6; NWP2 n = 6). Detailed information on all samples in this study is provided in Supplementary Table 1.
Cryptic species identification. We used the methods described by Shen and his research group 15,16 to identify cryptic species among the observed samples. In general, we extracted DNA from fish tissue using the DNeasy Blood & Tissue Kit (Qiagen, Venlo, Netherlands) according to the manufacturer's protocol and conducted a multiplex COI haplotype-specific PCR (MHS-PCR). Cryptic species were detected based on the length of PCR products run on gel agarose (2%). Six PCR product samples were sequenced using Sanger sequencing after being cloned into plasmid vectors. Nucleotide BLAST on Genbank was used to identify sequences. The cloned genes showed more than 97% identity with the COI gene of grey mullet. In addition, a phylogeny tree (not shown) was constructed, and cryptic species samples identified in our study were clustered with corresponding cryptic species from previous studies. As we used the method from previous studies 15, 16 , we decided to keep the same abbreviations for the three grey mullet cryptic species(i.e. NWP1, NWP2, and NWP3). The list of fish samples is presented in Supplementary Table 1 and is made up of 16 adult fish samples-eight NWP1, four NWP2, and four NWP3-and 12 juvenile fish samples-six each of NWP1 and NWP2. We did not observe any juvenile NWP3 samples, probably because NWP3 juveniles are scarce in Taiwanese waters (unpublished data). However, comprehensive data on juvenile and adult samples of NWP1 and NWP2 suggested that more information is needed on the differences in gut microbial diversity among the different life stages, as well as the potential influence of migratory tracts on the intestinal microbiota.
Historical ocean water microflora data. The three cryptic species used in this study are believed to migrate from the East and South China Seas to the Taiwan Strait during their spawning season 15 . In particular, fertile NWP1 were proposed to migrate to the Taiwan Strait following the China coastal current, while NWP2 and NWP3 individuals reach their spawning area by following the Kuroshio Current. Therefore, to investigate the potential impacts of migration on the structure of gastrointestinal microbiota in migratory grey mullets, historical seawater microbial community data on grey mullet migration routes were adapted from four previous studies [18][19][20][21] . The list of cited data and detailed information on the samples are listed in Supplementary Table 4. We specifically chose data from those studies because their methods were similar, using 16S rRNA amplicon sequencing to clarify the microbial community structure from environmental samples. Raw 16S rRNA sequences from seawater samples were combined with our raw fish gut microbiota data and analysed together following the sequence analysis process.
Taxonomic assignment and statistical analysis. All raw paired-end Illumina amplicon sequences in FASTQ format were processed with MOTHUR v.1.36.1 50 . During pre-processing, raw paired-end sequences (including forward and reverse sequences) were merged using the fastq_mergepairs command in USEARCH. Next, the sequences were quality-trimmed using command trim.seqs in Mothur (www.mothur.org/wiki/Trim. seqs) with a quality score of 30, number of errors in the barcode (maximum of 1), number of errors in the primer (maximum of 2), and length from 250 to 500 bp 51 . After picking a reference sequence using the command unique. seqs in Mothur, the partial 16S rRNA sequences were hierarchically clustered at 97% using command cluster_otus in USEARCH, which also removes chimeras from all libraries. Final OTUs were taxonomically classified using the classify.seqs script in Mothur against GreenGenes 13_8 database as a reference 52 . The final step was to remove OTUs derived from the Chloroplast, Mitochondria, Eukaryota, or unknown kingdoms and to generate an OTU table.
Measurements of alpha diversity-including observed number of OTUs, richness (Chao), and community diversity (Shannon and inverse Simpson's index)-were calculated using the collect.single script in Mothur. The permutation T-test in R using the RVAideMemoire package 53 tested for differences in alpha diversity between adults of different cryptic species as well as between adults and juveniles of the same cryptic species. A p value < 0.05 was considered statistically significant.
Microbiota β-diversity was studied using the Bray-Curtis dissimilarity index and weighted UniFrac distance. Later, β-diversities of grey mullet gut microbial communities were visualized by generating non-metric multidimensional scaling (NMDS) plots in R using the vegan package 54 . A permutational multivariate analysis of variance (PERMANOVA) test was conducted using the pairwiseAdonis package 55 to compare microbial community structures among the cryptic species as well as inter-specific cryptic species at different life stages.
Next, the phylogenetic tree of shared OTUs within cryptic species was constructed using Maximum Likelihood analysis based on the Kimura 2-parameter model in MEGA 7 56 . Box plots, bar charts and heatmaps were created using R and STAMP 57 . Geospatial biodiversity mapping analysis of fish gut microbial compositions and historical seawater microflora was performed using GenGIS II 58 . This method was used to determine the