Phylogeographic patterns in the Philippine archipelago influence symbiont diversity in the bobtail squid–Vibrio mutualism

Abstract Marine microbes encounter a myriad of biotic and abiotic factors that can impact fitness by limiting their range and capacity to move between habitats. This is especially true for environmentally transmitted bacteria that cycle between their hosts and the surrounding habitat. As geologic history, biogeography, and other factors such as water temperature, salinity, and physical barriers can inhibit bacterial movement to novel environments, we chose to examine the genetic architecture of Euprymna albatrossae (Mollusca: Cephalopoda) and their Vibrio fischeri symbionts in the Philippine archipelago using a combined phylogeographic approach. Eleven separate sites in the Philippine islands were examined using haplotype estimates that were examined via nested clade analysis to determine the relationship between E. albatrossae and V. fischeri populations and their geographic location. Identical analyses of molecular variance (AMOVA) were used to estimate variation within and between populations for host and symbiont genetic data. Host animals demonstrated a significant amount of variation within island groups, while symbiont variation was found within individual populations. Nested clade phylogenetic analysis revealed that hosts and symbionts may have colonized this area at different times, with a sudden change in habitat. Additionally, host data indicate restricted gene flow, whereas symbionts show range expansion, followed by periodic restriction to genetic flow. These differences between host and symbiont networks indicate that factors “outside the squid” influence distribution of Philippine V. fischeri. Our results shed light on how geography and changing environmental factors can impact marine symbiotic associations at both local and global scales.


| INTRODUC TI ON
The dispersal of marine species across suitable habitats can be affected by physical barriers (temperature, distances across oceans, island formations) as well as life history strategies (e.g., dispersal method of larvae and adult motility; (Kool, Paris, Barber & Cowen, 2011). Biogeographic barriers, as reported by floral and faunal separations, occur worldwide and provide an opportunity to study how physical barriers coupled with other abiotic factors may be affecting species dispersal and ultimately distribution (Lohman et al., 2011;Tonon et al., 2015). Analysis of population structure and physical orientation of the distribution of taxa across these barriers has given us clues to the factors that fragment available habitat . Although previous work has provided evidence for several causes for speciation among closely related populations in areas where distinct barriers exist, there is less known about species that coexist with one another and whether rules that govern distribution patterns via allopatric speciation influence such associations (Hellberg, Burton, Neigel & Palumbi, 2002;Palumbi, 1994).
One region that has been studied extensively for its unique patterns of biogeography and geologic history is the Indo-Pacific barrier (IPB), which was created by the uprising of the Indonesian archipelago separating the Indian and Pacific oceans (Gaither, Toonen, Robertson, Planes & Bowen, 2010). Interestingly, dispersal mechanisms and rapid adult motility have allowed certain taxa in the region to cross the IPB due to various dispersal strategies and larval residence time prior to metamorphosis compared to other taxa which are geographically restricted (Horne, 2014;Liu, Chang, Borsa, Chen & Dai, 2014;Sorenson, Allen, Erdmann, Dai & Liu, 2014). As part of the IPB, the Philippine island archipelago is a "hotspot" for species diversity and endemism and has warranted investigation of the distribution of taxa across the region (Roberts et al., 2002). In the Philippines, the current research has focused on the phylogeographic distribution of some fishes, bent-toed geckoes, as well as bivalves across established biogeographic margins that limit some other terrestrial and marine taxa (Carpenter & Springer, 2005;Esselstyn et al., 2010;Gaither & Rocha, 2013;Huxley, 1868;Lemer et al., 2016;Siler, Oaks, Esselstyn, Diesmos & Brown, 2010;Wallace, 1860Wallace, , 1863. Local analysis of the distribution and connectivity of some marine taxa across the Philippines has also been investigated in western populations of the sea star Linckia laevigata and the giant clam Tridacna crocea near the island of Palawan, as well in the western portion of the Central Visayas (Alcazar & Kochzius, 2016;Juinio-Menez, Magsino, Ravago-Gotanco & Yu, 2003;Magsino, Ravago & Juinio-Menez, 2002; Ravago-Gotanco, Magsino & Juinio-Menez, 2007) Across the globe, sepiolid squids (Cephalopoda: Sepiolidae) form mutualistic associations with bioluminescent bacteria from the genera Vibrio and Photobacterium (γ-Proteobacteria: Vibrionaceae (Herring, 1977). Vibrio bacteria are housed within a specialized internal organ called the light organ (LO), where the host provides a nutrient-rich habitat for the symbiont, and in return, Vibrio bacteria provide luminescence to the squid to be used in a behavior termed counterillumination (Jones & Nishiguchi, 2004). Squid hosts use the Vibrio-produced light to reduce their silhouette during the evening, which enhances their survivability and predation success (McFall-Ngai, Heath-Heckman, Gillette, Peyer & Harvie, 2012). After each nightly foraging session, approximately 95% of the Vibrio bacteria are vented out into the surrounding seawater, seeding the local area with symbiotically viable Vibrios (Boettcher, Ruby & McFall-Ngai, 1996). Local cycling of symbiotic V. fischeri exposes these bacteria to a wide range of abiotic and biotic factors outside the host that can affect their fitness and ability to infect new hosts. This also allows for symbiotically competent free-living bacteria to migrate to new host habitats, where they can invade and colonize different populations of sepiolids (Nyholm, 2004;Nyholm & Nishiguchi, 2008;Nyholm, Stabb, Ruby & McFall-Ngai, 2000).
Earlier work on sepiolid squids has focused on the influence of geographic distance on symbiont prevalence and genotype in both sympatric and allopatric populations (Jones, Lopez, Huttenburg & Nishiguchi, 2006;Kimbell, McFall-Ngai & Roderick, 2002;Zamborsky & Nishiguchi, 2011). Allopatric and sympatric populations for both squids and Vibrio bacteria show distinct population breaks that are not necessarily driven by host specificity. Additionally, host-mediated factors along with abiotic variables such as water temperature and salinity have been known to shape these mutualist assemblages (McFall-Ngai, 2014;McFall-Ngai et al., 2012;Nishiguchi, 2000;Soto, Gutierrez, Remmenga & Nishiguchi, 2009). Collectively, either genomic comparison of closely related populations (Bongrand et al., 2016) or haplotype comparisons of allopatric populations of Indo-west Pacific squid and their vibrio symbionts (Jones et al., 2006) do not address the connectivity of populations across physical and biogeographic barriers like those in the Philippines or across the IPB. Therefore, we examined the genetic architecture of Euprymna albatrossae (Cephalopoda: Sepiolidae) and their V. fischeri symbionts in the Philippine archipelago using a combined phylogeographic approach to determine whether host specificity or geographic location influence the distribution of symbiotic Vibrios in this region. The unique geographic origin of the Philippines, its proximity to deeper and colder water, as well as currents that move through the area allow for the investigation of what roles geography and host specificity have in the distribution of mutualistic associations.
F I G U R E 1 Sampling location where host squid were collected during the months of May, June, July, and August of 2010, 2012, 2013, and 2015 TA B L E 1 Geographic location, sampling years, and sites where Euprymna albatrossae and Vibrio fischeri were collected in the Philippines during the months of May, June, July, and August of the years listed below

| Specimen collection and bacterial isolation
Squids were collected in the months May, June, July, and August during the years 2010, 2012, 2013, and 2015 at eleven different sites around the Philippine islands ( Figure 1,

| DNA extraction and amplification
Euprymna albatrossae DNA was extracted using approximately 25 mg of ethanol preserved tissue that was dissected from the gill or mantle of each squid. Dissected tissues were washed with 100 μL of nuclease-free water to remove any residual ethanol. E. albatrossae DNA was extracted using the DNeasy© blood and tissue protocol for animal tissues (Qiagen, Valencia, CA). All genomic DNA extractions were visualized on a 1% agarose gel and quantified using a Nanodrop 9600 (ThermoFisher Scientific, Waltham, MA).
Isolation of DNA from V. fischeri light organ isolates was completed using the Qiagen DNeasy © blood and tissue kit (Valencia, CA) Gram-negative bacterial protocol. Approximately 2 × 10 9 cells were transferred from each log-phase culture to the extraction tube for centrifugation. After, the remaining pellet was used for extraction using the Qiagen protocol. Purified V. fischeri DNA was visualized on a 1% agarose gel and quantified using a Nanodrop 9600 (ThermoFisher Scientific, Waltham, MA). Isolated DNA extracted from each V. fischeri isolate was used to amplify a portion of the glyceraldehyde phosphate dehydrogenase (gapA) locus (~900 bp) by PCR, using previously described Vibrio-specific primers (Table 2; (Jones et al., 2006;Nishiguchi & Nair, 2003). The gapA locus has been used reliably to estimate deep phylogenetic connections between bacterial families (Nelson, Whittam & Selander, 1991) within the Vibrionaceae (Thompson, Gomez-Gil, Vasconcelos & Sawabe, 2007) as well local population structure of mutualist V. fischeri (Jones et al., 2006;Nishiguchi & Nair, 2003

| Haplotype networks, nested clade analysis, and molecular variance
Haplotype networks for squid and symbiont were generated using TCS v1.12 using statistical parsimony methods outlined by Templeton (Templeton & Sing, 1993). Nested clade analyses were performed using Templeton's nesting algorithm as implemented in GEODIS (Posada, Crandall & Templeton, 2000). Analysis of molecular variance (AMOVA) was executed using the population genetics software platform ARLEQUIN (Excoffier & Lischer, 2010). Analyses were run for measures of within-and among-population variation along with a separate analysis assessing variation by island for both host and symbiont data.
Concurrently, theta (Θ), a base-pair-by-base-pair measure of polymorphism was calculated for each mutualist population at each sample site.    (Figures 2 and 3). The largest haplotype (haplotype 1) contains a significant number of Cebu haplotypes coupled with F I G U R E 2 TCS nested haplotype network generated from Euprymna albatrossae molecular data acquired from animals captured in the Philippines during the years 2010-2015. Each line in the diagram represents one-base-pair mutational step between haplotypes. Black circles represent unsampled mutational steps connecting haplotypes. The size of each circle is indicative of the number of sequences that make up that haplotype, with the color of each circle representing the geographic origin of the sequence data and its proportion of the total haplotype. Each haplotype is represented by a two-digit indicator with the dotted line enclosures indicating the nesting hierarchy. Each nesting level is labeled with a dashed two-digit label populations from Negros and Palawan. Each of the major haplotypes listed requires a minimum of one-base-pair change, with the largest number of changes needed (6) to go from haplotype 1 to haplotype 8.

| Nested clade and molecular variance analysis of V. fischeri
Contingency analyses of symbiotic V. fischeri nesting revealed significant evidence for restricted gene flow with isolation by distance in clades 2-1 and 2-8 (Table 5) (Table 5).
An identical AMOVA of symbiont genetic data revealed that a significant portion of the variance exhibited by these populations exists within and among populations (14.40%, 80.08%; Table 3).

| Host genetic architecture
The genetic structure of the E. albatrossae sampled for this study indicates that geographic location impacts host distribution. Island  Table 3).
The genetic fixation observed in the host genetic data, reported as F ST , indicates that genetic flow is limited throughout the region and Notes. Location of significance is indicated by (D n ), nested clade distance, and/or (D c ), the within-clade distance. I-T indicates the average distance between a tip clade and an interior clade. S or L indicates that the distance measure is significantly smaller or larger at the 5% inference level. Inference steps were performed using the automated inference key in GEODIS, part of the AneCA v1.2 population genetics analysis software platform (Posada et al., 2000).
populations are genetically isolated from each other (F ST = 0.70897, p < 0.0001; Table 3). This is reflected in the three distinct host networks detected in the Philippines (Figure 2) and suggests that geography may be influencing host genetic exchange and distribution. This could be because the benthic lifestyle adult Euprymna squid (~2-4 cm) lead as adults rarely migrate; however, the semipelagic nature of newly hatched squid (~3-5 mm) would allow water flow to transport juveniles to novel locations (Kimbell et al., 2002;Villanueva, Vidal, Fernández-Álvarez & Nabhitabhata, 2016). Two of the distinct networks occur in the central island chain and, while having similar geologic origin, show no genetic connection in habitats that are homogeneous . Clade 4-7 is distinctly made up of mostly squid haplotypes detected around the island of Panay, with small contributions from populations around Negros (haplotypes 21, 22, 23, and 33, Figure 2). This indicates that populations may have been fragmented due to the result of geologic activity in the region, seasonal changes in currents, or even modern-day commercial fishing management, which have all been shown to influence fragmentation of marine habitats in the area (Abesamis, Russ & Alcala, 2006;Huang, Wu, Zhao, Chen & Wang, 1997;Savina & White, 1986;Wyrtki, 1961;Zhou, Ru & Chen, 1995). Habitat fragmentation was also inferred within clades 2-11 and 3-1 (Table 4) (Table 4), suggesting that the alternating direction of prevailing currents in the region is one mechanism of dispersal as well as isolation for these squids.
F I G U R E 3 TCS haplotype network generated from Vibrio fischeri molecular data acquired from isolates harvested from squid light organs in the Philippines during the years 2010-2015. Each line in the diagram represents one-base-pair mutational step between haplotypes. Black circles represent unsampled haplotypes. The size of each circle is indicative of the number of sequences that make up that haplotype, with the color of each circle representing the geographic origin of the sequence data and its proportion of the total haplotype. Each haplotype is represented by a two-digit indicator The divergence of the equatorial current (EC) as it approaches the Philippines from the east influences the directional flow of water through the central island chain, particularly north of the island of Panay and south of the islands of Negros and Cebu ( Figure 5; (Huang et al., 1997;Wyrtki, 1961). The amount and speed of water that is funneled around or through the central islands depends on the time of year and is reflected in the patterns detected in host haplotype networks (Wyrtki, 1961). In late winter/early spring, water flows from the east to the west from the San Bernardino Strait, the south of Masbate, and finally around the north of Panay toward the Sulu Sea (Figure 5a). This flow pattern coupled with the flow of the southern divergence of the EC allows for genetic exchange between geographically isolated populations from Panay, Cebu, and Negros While clade 4-1 had no valid inference (Table 4), clade 3-1, which is fully nested within 4-1 (Figure 2), had an inference of F I G U R E 4 Nested Vibrio fischeri haplotype network generated from molecular data acquired in the Philippines during 2010, 2012, 2013, and 2015. Each haplotype is represented by a two-digit identifier (see Figure 3), with each hierarchical nesting level represented by a twoto three-digit dashed identifier and enclosed within dashed and dotted lines. Lines between haplotypes represent the mutational steps required to transition from one genetic station to another, with the small black dots representing unsampled haplotypes restricted gene flow with isolation by distance. This indicates that these habitats may have changed over time and influenced the genetic flow between once connected populations of host squid. A third separate network (clade 4-2, Figure 2) consisted solely of samples collected around the western island of Palawan provides evidence that the unique geologic origin of Palawan may have fragmented a once continuous population of host animals (Sathiamurthy & Voris, 2006;Wallace, 1863;Zhou et al., 1995).

| Symbiont genetic architecture
Symbiont genetic data indicate that Vibrio bacteria seem to be able to mitigate the barriers that restrict host genetic exchange. Analysis (AMOVA) of total symbiont genetic data reveals that most of the genetic variation observed lies within each population in contrast to the partitioning of variation from host genetic data (df = 186, SS = 438.538, VC = 2.6704, PV = 80.08%; Table 3). The level of genetic diversity of the total symbiont population in the region is also Notes. Location of significance is indicated by (D n ), nested clade distance, and/or (D c ), the within-clade distance. I-T indicates the average distance between a tip clade and an interior clade. S or L indicates that the distance measure is significantly smaller or larger at the 5% inference level. Inference steps were performed using the automated inference key in GEODIS, part of the AneCA v1.2 population genetics analysis software platform (Posada et al., 2000). NS, not significant.
indicative of symbiont gene flow between populations of hosts that are isolated from one another (F ST = 0.1991, Table 3). One contiguous haplotype network was detected in the symbiont genetic data (Figures 3 and 4) revealing several genetic connections between symbionts collected from host squid that are genetically and geographically isolated from one another (Figure 2). The predominant haplotypes found within symbiont genetic data (haplotypes 1, 3, 8, 57, and 58;

| Geography, geologic history, and environmental conditions
Results from this study indicate that geography plays a role in host squid distribution, without demonstrating a significant influence on symbiont distribution. The disparity in these patterns may be a result of differences in dispersal methodology between mutualist partners; that is, host squid have a limited range as adults and rarely travel far from their birthplace due to the limited dispersal ability of direct developing, benthic hatchlings (Kimbell et al., 2002;Villanueva et al., 2016). Conversely, symbiotically viable Vibrio bacteria are cycled out of the host daily exposing them to environmental factors (i.e., currents) that allow for movement into novel areas where they are able to recruit into a new host. While bacteria alone cannot cross great expanses of ocean, the use of rafting has been shown to be an effective dispersal mechanism for marine bacteria like V. fischeri (Jones et al., 2006;Theil & Gutow, 2005). The ability for vibrios to cross great expanses of oceans has been previously reported in other marine bacteria and undoubtedly will allow symbiotically viable vibrios to be shuttled to new areas and novel hosts (González-Escalona, Gavilan, Brown & Martinez-Urtaza, 2015).
Prevailing currents in and around the central Philippine island chain vary in direction and magnitude seasonally (Wyrtki, 1961). As E. albatrossae breed all year long, this change in directionality may provide newly hatched squid the opportunity to be carried to new areas, despite their otherwise limited dispersal ability, while being cut off from other available habitats when the prevailing currents change (Hanlon, Claes, Ashcraft & Dunlap, 1997). The pattern of direction in the symbiont genetic data presented here indicates that introgression across the Sulu Sea, which appears to be a biogeographic margin for host animals, is facilitated by the directional flow of water during the monsoon season (Huang et al., 1997). Euprymna hatchlings are known to be "pelagic"; that is, they linger in the water column before settling to their benthic lifestyle (Moltschaniwskyj & Doherty, 1995). This also might heighten the ability of host populations to move to new localities.
Geologic changes and the physical oceanography of this region may also explain the patterns detected in the genetic data. Glacial maximum sea levels exposed portions of what was host native range within the central island chain, creating a disconnect between populations in the west and central island squid assemblages (Gaither & Rocha, 2013;Gordon, 2005;Zhou et al., 1995). During glacial norms, hosts are restricted by a deep-water thermocline that has persisted since before the Holocene, between Palawan and the central island chain (Miao et al., 1994). Fluctuating sea level during glacial cycles as well as Cenozoic volcanic uprising of the central Visayas may also explain the disjunctive distribution of host animals across this region (Miao et al., 1994;Zhou et al., 1995). Given that many of the more abundant haplotypes examined have prevalence in localities that are geographically distinct, this provides additional evidence that host populations have been established well beyond the geologic history of the Philippines (e.g., Palawan).
While previous research has shown that symbiont gene flow can be restricted by temperature, symbionts in this region seem to be able to mitigate environmental barriers which hosts cannot, crossing geographic and biologic barriers with apparent ease (Nishiguchi, 2000). Symbiont gene flow demonstrates a current dependent directionality of introgression by vibrios from the central islands west to Palawan in the winter and west to east in the summer months (Huang et al., 1997). The El Niño Southern Oscillation has also been shown to influence not only sea surface temperatures, wind direction, and rainfall in this region, but also the position of this deep-water thermocline, further isolating local populations of host squid while not restricting symbiont distribution (Chen et al., 2015;Stuecker et al., 2013). While other Indo-west Pacific and Mediterranean populations of Vibrio demonstrate that some degree of host specificity, geography, or other environmental factors can impact symbiont genetic architecture, findings from this study indicate that geography alone cannot explain symbiont distribution and that physical factors (e.g., currents) are important drivers of microbial diversity in the region (Jones et al., 2006;Zamborsky & Nishiguchi, 2011).
Beneficial associations like the sepiolid squid-Vibrio mutualism will undoubtedly be impacted by the reduction of available habitat, highlighting the importance of investigating the influence geography has on symbiont prevalence and distribution. Findings from this study point to a need to better understand the mechanisms that will impact symbiotic associations across a changing landscape and what factors will influence the fitness of beneficial microbes when they are moved to a novel habitat. Our findings have provided clues as to how established populations of host squids are the foundation for symbiont population structure, yet abiotic factors still influence where vibrios can move and establish new populations.