Symbiodiniaceae diversity varies by host and environment across thermally distinct reefs

Endosymbiotic dinoflagellates (Symbiodiniaceae) influence coral thermal tolerance at both local and regional scales. In isolation, the effects of host genetics, environment, and thermal disturbances on symbiont communities are well understood, yet their combined effects remain poorly resolved. Here, we investigate Symbiodiniaceae across 1300 km in Australia's Coral Sea Marine Park to disentangle these interactive effects. We identified Symbiodiniaceae to species‐level resolution for three coral species (Acropora cf humilis, Pocillopora verrucosa, and Pocillopora meandrina) by sequencing two genetic markers of the symbiont (ITS2 and psbAncr), paired with genotype‐by‐sequencing of the coral host (DArT‐seq). Our samples predominantly returned sequences from the genus Cladocopium, where Acropora cf humilis affiliated with C3k, Pocillopora verrucosa with C. pacificum, and Pocillopora meandrina with C. latusorum. Multivariate analyses revealed that Acropora symbionts were driven strongly by local environment and thermal disturbances. In contrast, Pocillopora symbiont communities were both partitioned 2.5‐fold more by host genetic structure than by environmental structure. Among the two Pocillopora species, the effects of environment and host genetics explained four times more variation in symbionts for P. meandrina than P. verrucosa. The concurrent bleaching event in 2020 had variable impacts on symbiont communities, consistent with patterns in P. verrucosa and A. cf humilis, but not P. meandrina. Our findings demonstrate how symbiont macroscale community structure responses to environmental gradients depend on host species and their respective population structure. Integrating host, symbiont, and environmental data will help forecast the adaptive potential of corals and their symbionts amidst a rapidly changing environment.


| INTRODUC TI ON
Tropical coral reefs rely on endosymbiotic algae, Symbiodiniaceae (LaJeunesse et al., 2018), that provide nutrition to support the function and survival of reef-building corals (Baker, 2003;Muscatine, 1990).Even a slight increase of 1-2°C above current sea surface temperatures (SST) can disrupt the delicate balance between corals and their symbionts (Weis, 2008;Weis et al., 2008), resulting in cascading negative effects on reef ecosystems (Baker et al., 2008;Graham et al., 2007).Symbiodiniaceae play a critical role in mediating coral susceptibility to thermal stress through adaptive and acclimatory mechanisms (Manzello et al., 2019;Marhoefer et al., 2021;Morikawa & Palumbi, 2019) that vary according to coral taxa (Quigley et al., 2022) and divergent symbiont lineages (Abrego et al., 2008;Butler et al., 2023).The increased frequency and duration of marine heatwaves have already impacted coral reefs globally (Hughes et al., 2017;Oliver et al., 2019).Thus, understanding the evolutionary relationships between corals and their symbionts can transform our understanding of adaptive traits in response to rapid environmental change.
The stability and flexibility of the coral symbiosis relies on the mode of symbiont transmission to coral offspring (Baker, 2003).
Corals which vertically transmit symbionts directly to eggs often mediate strong co-diversification over long evolutionary periods (Turnham et al., 2021), resulting in stronger host effects on symbionts and adaptive responses of corals to their environment (Seah et al., 2017).In contrast, horizontally transmitting coral species produce aposymbiotic larvae which acquire symbionts from the surrounding environment, promoting greater flexibility in symbiont selection and enabling rapid acquisition or exchange of symbionts better suited for their environment (Baird et al., 2007;Cumbo et al., 2013).In the face of a disturbance, horizontal transmitting corals harbour a greater ability to shuffle or switch dominant and background symbiont species (Elder et al., 2023;Jones & Berkelmans, 2010), though can be disadvantageous during periods of stress if the symbionts are not well-adapted to variable or extreme conditions due to a lesser ability for host-symbiont co-diversification (LaJeunesse, Smith et al., 2010).There is also mounting evidence that some vertical transmitting corals have a mixed-mode of transmission, complicating the symbiont-host relationship over the lifespan of a coral (Quigley et al., 2017, Quigley et al., 2018;Starko et al., 2024).Therefore, resolving the effect of host species on symbiont community structure requires parallel investigation of multiple species with similar demographics.
The transmission strategy of the coral host also influences the extent to which the environment structures symbiont communities (van Oppen & Medina, 2020).Vertically transmitting species often maintain a stronger host genotypic effect as the primary driver of symbiont community structure in comparison to strong environmental gradients (Hume et al., 2020;Voolstra et al., 2021).In horizontally transmitting species, Symbiodiniaceae community changes are generally structured by their environment, including variable thermal fluctuations or chronic thermal disturbances (e.g., mass bleaching events) (Howells et al., 2020;Kennedy et al., 2016;Quigley et al., 2022).In addition to strong environmental effects, there is evidence of secondary host effects on symbiont community structure in horizontally transmitting corals such as Acropora spp.(Cooke et al., 2020;Matias et al., 2023), alluding to both host and environment influences on symbiont community structure.Thus, it is likely that while vertical and horizontal transmission of the host is the dominant structuring predictor in symbiont communities, there are underlying influences of the environment which interact with host genotype.The analysis of both host and symbiont genetic diversity in parallel with contrasting environmental gradients can provide insight into how host species, genotypes, and environment co-regulate the impacts on Symbiodiniaceae community structuring.
To address this knowledge gap, we investigate the influence of coral host genetics and environment on symbiont macroscale community structure.We define symbiont macroscale community structure as the diversity of Symbiodiniaceae among conspecific hosts living along an environmental continuum (hereafter defined as 'symbiont communities') (Davies et al., 2023;LaJeunesse, Pettay et al., 2010).The diversity and distribution of Symbiodiniaceae in relation to host population structure and the environment, in particular, current and historical thermal stress, enables the identification of key genetic and environmental drivers that affect the capacity of corals and their symbionts to cope with climatic extremes into the future.To achieve this, we used ITS2 metabarcoding and a supporting genetic marker (psbA ncr ) to analyse patterns of symbiont community composition, in parallel with a genotype-by-sequencing approach to describe trends of host population structure and hostsymbiont co-phylogeny.Our approach targeted three coral species, two with vertical symbiont transmission, Pocillopora verrucosa and P. meandrina, and one with horizontal transmission, Acropora cf humilis.All coral species were sampled across a 1300 km gradient in Australia's Coral Sea Marine Park.Specifically, we tested the effects of biogeographic variables (latitude), historical thermal disturbances (maximum Degree Heating Weeks [maxDHW]), light attenuation in seawater (kd490), and individual-specific variables (depth and host genetic cluster).Samples were collected during the height of the mass bleaching event in 2020, therefore, we also accounted for accumulated heat experienced at each reef (DHW at collection) and the bleaching condition of each coral colony in relation to symbiont community structure.

| Sampling location and study species
Coral samples were collected on SCUBA from 13 reefs within the Coral Sea Marine Park (CSMP) between February 16th and March 12th, 2020 (SOM, Table S1).Reefs in the CSMP have historically been exposed to a range of disturbances including regular and severe marine heatwaves (Harrison et al., 2019).There are currently no records for Symbiodiniaceae genetic diversity in the CSMP, despite this region covering nearly 1 million km 2 and connecting fauna of the Great Barrier Reef to other western Pacific reefs (Payet et al., 2022;van Oppen et al., 2008).For this study, a total of 609 coral samples were collected from two genera, including 346 Pocillopora samples and 262 Acropora samples.Samples were collected at an average depth of 7.8 m (± 2.7), ranging between 2 and 16 m and sampled during the 2020 bleaching event where reefs had been exposed to between 5.7 and 10.0 DHW at the time of collection (Marzonie et al., 2023).Corals were sampled during or immediately after the peak of the marine heatwave (Figure S1).Within each reef, a range of bleaching phenotypes for each coral species was collected.The visual bleaching score of each coral colony was first scored against a six-point Coral Watch health chart (Siebeck et al., 2006) and collection depth was recorded in situ.Corals were photographed using an Olympus TG5 at the whole colony level and at a macro-scale to assist in species identification.Coral fragments were then collected using a hammer and chisel and each individual coral fragment was placed into a resealable, numbered Ziploc bag.After each dive, coral fragments were preserved in 100% ethanol and exchanged twice.

Symbiodiniaceae identification
DNA was extracted from coral samples of Pocillopora and Acropora using a modified version of Wayne's method (Wilson et al., 2002) (Appendix 1).Following extractions and after a minimum of 24 hours, DNA concentrations were quantified using a Qubit dsDNA High Sensitivity Assay Kit (Invitrogen) with a Qubit 3 Fluorometer.
Samples were then standardized to 10 ng μL −1 using an automated pipettor (QIAgility, QIAGEN).Normalized DNA was used to amplify two marker regions for Symbiodiniaceae species confirmation, the ITS2 and psbA ncr regions.Detailed PCR protocols for each marker can be found in Appendix 1.

Symbiodiniaceae identification
The Internal Transcribed Spacer 2 (ITS2) region was the primary marker used to distinguish Symbiodiniaceae.Amplicon sequencing of the ITS2 region was performed on an Illumina MiSeq platform at 2 × 300 bp paired-end V3 chemistry (Ramaciotti Centre for Genomics, UNSW).One MiSeq run was conducted for Acropora cf humilis and a separate run was done for Pocillopora spp., with each run conducted on one flow cell.Within each run, samples from each reef were randomly distributed across each plate.The raw sequence data are available under NCBI BioProject PRJNA1001407.
Demultiplexed forward and reverse Fastq files were submitted to SymPo rtal.org for analysis of 'Defining Intragenomic Variants' (DIVs) and ITS2 type profiles (Hume et al., 2019).The SymPortal analytical pipeline assesses patterns of intragenomic variation for the ITS2 marker as there can be multiple copies of each gene present within Symbiodiniaceae genomes.These sequences are then partitioned as 'defining intragenomic variants' (DIVS), from which 'ITS2 type profiles' are derived as discrete biological entities.Using SymPortal outputs, three samples with <1500 reads were removed from the data set prior to community analyses.Filtered datasets for DIVs and type profiles were used in downstream statistical analyses.
Coupled with the ITS2 marker, the non-coding region of the plastid, mini-circle psbA gene (psbA ncr ) often contains hyper-variable nucleotide sequences well suited for Symbiodiniaceae species confirmation, which can be used in tandem with the ITS2 marker (LaJeunesse & Thornhill, 2011;Moore et al., 2003) to validate ITS2 type profiles (Smith et al., 2017).Amplicons were sequenced in the forward and reverse direction using Sanger chemistry (Macrogen Inc.) and sequences were manually trimmed and aligned using DECIPHER (Wright, 2016).

| Construction of phylogenetic trees, GUniFrac distances, and ITS2-psbA ncr tanglegrams
All statistical analyses were conducted in the R statistical environment (R Core Team, 2022) and code is available online as a live GitHub repository (https:// github.com/ magen amarz onie/ Coral SeaSy mbiont) and as a static release on Zenodo (https:// doi.org/ 10. 5281/ zenodo.10820899).We used post-MED (Minimum Entropy Decomposition) sequences from SymPortal outputs to incorporate rare variants into the analyses (Hume et al., 2019).We opted for a k-mer based approach to produce pairwise sequence comparisons of SymPortal post-MED sequences (Fujise et al., 2021), given the challenges of multiple sequence alignment of Symbiodiniaceae ITS2 sequences.A k-size of 7 was chosen with default settings using 'kdistance' in the kmer package (Wilkinson, 2018).We then produced a UPGMA phylogenetic tree of ITS2 sequences from the pairwise k-mer distances with 'UPGMA' in Phangorn (Schliep, 2011;Schliep et al., 2017).This tree, along with post-MED sequence counts, was further analysed using Generalized UniFrac (GUniFrac, α = 0.5) intersample distances, a method considered to apply fair weighting to both rare and abundant sequences (Chen et al., 2012).The GUniFrac distance approach was selected as it partitioned samples into clusters congruent with SymPortal ITS2 type profiles, while allowing for the non-DIV/non-profile sequences to contribute to inter-sample distances.psbA ncr trees were then generated using the k-mer-based approach as above and untangled using the step2side method in dendextend (Galili, 2015) to investigate psbA ncr -ITS2 congruence.

| Identification of Pocilloporid species
Given the cryptic nature of Pocilloporids and their morphological plasticity in variable environments (Johnston, Wyatt, et al., 2022), specimens were genetically resolved to species-level after collection.
To confirm species identification, we amplified the mitochondrial open reading frame (mtORF region) of the coral host and applied a Restriction Fragment Length Polymorphism (RFLP) assay modified from Johnston et al. (2018) and detailed in Marzonie et al. (2023).
Of the 346 Pocillopora samples, 152 were identified as P. verrucosa, 134 as P. meandrina, and 61 ambiguous samples, which we refer to as 'unknown' Pocillopora spp.These remaining unknown samples of Pocillopora that were not identified as P. meandrina or P. verrucosa in the RFLP assay were then sequenced in both directions using Sanger Sequencing.Using the amplified mtORF region, samples underwent bi-directional Sanger Sequencing (Macrogen, South Korea).
We manually applied sequence trimming parameters based on sequence quality and generated consensus sequences using DECIPHER (Wright, 2016).These sequences were aligned with the mtORF reference sequences collated from Pinzón et al. (2013), Forsman et al. (2013), and Gélin et al. (2017) using the 'AlignSeqs' function (DECIPHER).Hamming distances were computed using the 'dist.hamming'function from Phangorn to identify any unknown Pocillopora samples to previously described haplotypes or species (Figure S2).Of the 61 unknown Pocillopora samples, 39 were identified as Pocillopora haplotype 8a, which is proposed to comprise a single species according to nuclear DNA (Johnston et al., 2022a).An additional 22 samples were either unable to be mapped to previously described Pocillopora species, haplotypes, or were removed due to their low representation (e.g., two samples of Pocillopora damicornis cf acuta).These 22 samples were removed from downstream analyses.

| Coral host genotype-by-sequencing
To assess the effects of host genotype and population structure on symbiont communities, host coral tissue was sequenced using a genotype-by-sequencing approach (DarT-sequencing; Kilian et al., 2012) at Diversity Arrays Technology (Canberra, Australia).
For Pocillopora, all 61 unknown samples from the RFLP analysis were removed from the genotype-by-sequencing analysis, including host samples of haplotype 8a.A separate run was conducted for each host species (A.cf humilis, P. verrucosa, and P. meandrina).To reduce technical bias, samples from each reef were randomly distributed across 96-well plates.Libraries were constructed using the PstI and HpaII compatible adaptors with two restriction enzyme overhangs (Sansaloni et al., 2011) and sequencing was performed on Illumina HiSeq2500 across two lanes.For Acropora cf humilis, samples were mapped to the Acropora tenuis (Liew et al., 2016) and A. millepora (Fuller et al., 2020) genomes.For Pocillopora verrucosa and P. meandrina, samples were mapped to the P. verrucosa reference genome (Buitrago-López et al., 2020).
Quality filtering of SNPs was conducted for each host species separately in R software using the dartR package (Gruber et al., 2018).Filtering parameters are detailed in Appendix 1.We first assessed genetic pairwise comparisons (F ST ) of Pocillopora using a DarT-seq co-analysis to ensure the absence of cryptic species.F ST values between species were greater than between host clusters within a species.The pairwise comparisons between P. verrucosa and P. meandrina resulted in an F ST = 0.501, while the two clusters of P. meandrina had an F ST range between 0.00-0.06(Figure S3).Host population structure was visually checked using unconstrained ordination (PcoA, vegan), then ADMIXTURE analyses were performed to assess group membership probability of each individual sample using LEA (Frichot & François, 2015).The optimal K was determined by running cross-entropy models using the 'snmf' function in LEA where K = 1:10 were tested.The K with the lowest cross-entropy value was selected for each host species.
The host PC1 and PC2 from each host species PCA were extracted and used as covariates to represent a host genetic predictor in downstream multivariate analyses of symbiont communities (dbRDA).In addition, a Neighbour-Joining tree was run on host individual pairwise genetic differences using the same genetic distances as admixture and PCAs to identify variation among clusters using the ape package (Paradis & Schliep, 2019).Clusters from the Neighbour-Joining trees were consistent with both PCAs and admixture group membership probability (Figures S4-S6).

| Host-symbiont cophylogenetic analysis
To measure the degree of co-phylogeny between coral host and symbiont phylogenetic structure, we implemented a Procrustean rotation analysis.This approach allows a quantitative measurement of the degree of alignment between two distance matrices, where two ordinations are overlaid and rotated to obtain the smallest distance between two points representing different community structure.Specifically, Generalized Procrustean Analysis allows for the rotation and transformation between two datasets with variable dimensions (i.e., host and symbiont genetic distances) to measure the distance between two observations of the same sample (Dray et al., 2003;Hutchinson et al., 2017).As dimensions of the host and symbiont distance matrices must match in a Procrustes analysis, new phylogenetic trees were constructed for symbionts using the same parameters and packages as per Section 2.6, with a subset of samples corresponding to the host DarT-seq genetic samples.
For each host genetic data set, a distance matrix of the filtered SNP data was constructed in parallel using the 'euclidean' distance method (dartR).A Generalized Procrustes rotation analysis was implemented to align the GUniFrac distances of the symbiont with the SNP-based Euclidean distances of the host ('procrustes', vegan).
Residuals were checked and the distance between each host and symbiont was plotted using RDA.Significance was tested using 'protest' with 999 permutations per analysis (vegan).

| Statistical analysis of environmental and host genetic predictors
We performed a principal coordinates analysis (PcoA; 'cmdscale' function, vegan [Oksanen et al., 2022]) on the symbiont ITS2 GUniFrac distances separately for each host species.We visualized  S2).
To run the distance-based RDA models, we ran two separate conditional models for each species: an 'Environment' and a 'Host Genetic' model.Each model contained a subset of individuals from the full Symbiodiniaceae dataset with matching host genetic data.
For the 'Environment' models, host genetic PC1 and PC2 (from host genetic PCAs) were defined as conditional effects to account for the confounding influence of host genetics on the environment.We then reversed the model parameters for the 'Host genetic' models, where significant environmental predictors were included as conditional effects.For both the 'Environment' and 'Host Genetic' models, we inspected correlation plots using corrplot (Wei & Simko, 2021) and variance inflation factors (vegan).First, we removed one of two collinear variables where r > |.80|.Model parameters were adjusted to keep the most significant factors, while concurrently removing any variables with high correlation (VIF score >5) to reduce the effects of high collinearity in the models (SOM, Figure S9, S10, and S11).A backward stepwise model was then run to further reduce parameters using the 'ordistep' function with 999 permutations until stopping criteria were met for each dbRDA model (vegan).We tested the marginal effects of each predictor in the ordination to control for the order of the included covariates.Significance was assessed by running a permutational analysis of variance (PERMANOVA) on the 'ordistep' output (vegan).Scores and vector information for each dbRDA were extracted and plotted in ggplot2 (Wickham, 2016).1).Detailed SymPortal profiles for each coral host species can be found in Tables S3-S6.
Most coral samples had one ITS2 type profile associated with an individual.For Pocillopora, 324 samples contained a single ITS2 type profile, compared to four samples containing two ITS2 type profiles within a single host sample.Similarly, for Acropora, 250 samples contained a single ITS2 type profile, while 10 samples hosted two ITS2 type profiles.Pocillopora sequences had a lower proportion of non-profile sequences within each sample, on average making up 16.1% of each sample (Figure 1a) compared to 24.6% for Acropora samples.Within Pocillopora, P. verrucosa had a lower proportion of non-profile sequences (13.9%), while P. meandrina and haplotype 8a had a higher composition of non-profile sequences at 17.9% and 17.6%, respectively.

| Pocilloporid symbiont communities
SymPortal analysis of ITS2 sequences from Pocillopora yielded a total of 23 ITS2 profiles collectively and varied by host species (Table 1).
Across all Pocillopora species, the majority sequences, 'C1' and 'C42.2', were present in most samples, though the mean relative abundance varied by host species.Some detected sequences were highly specific to a coral host species.For example, 'C1d' was exclusive to P. verrucosa, found in 92.7% of samples, and averaged 37.2% relative abundance of reads per sample.There was a major split in the right side of the UPGMA tree associated with two distinct symbiont lineages in P. meandrina.Of the two lineages, 'C42u' was a diagnostic sequence for the left branch associated with P. meandrina, found in 34.3% of samples, but was present in high relative abundance, averaging 24% of reads per sample (Figure 1a).In contrast, the right branch for P. meandrina was driven more by the presence of 'C42g', found in 64.0% of samples but in a lower mean proportion of 15.7% of reads per sample (Figure 1a).Hereafter, the symbiont lineage associated with P. verrucosa will be referred to as C. pacificum (Turnham et al., 2021), while the two symbiont lineages associated with P. meandrina will be referred to as C. latusorum (Turnham et al., 2021) 'north' and 'south' to reflect the biogeographic pattern in their distribution (Figure 1a and read below).

| Acroporid symbiont communities
Sequencing of A. cf humilis symbiont communities yielded 23 ITS2 type profiles across all samples, the most common being C3k/C3-C50a-C29-C21ab-C3b. The co-dominant majority sequences 'C3k/C3' were present in 99.2% of all Acropora samples, representing on average 74.8% of reads in these samples.There was a major split in the UPGMA tree for Acropora samples (Figure 1b).The two branches within the Acropora UPGMA tree were generally congruent with the proportion of non-profile sequences.The left branch had a relative abundance of 20.7% non-profile sequences averaged across samples (Figure 1b, n = 122), significantly lower than 29.0% in the right branch (Figure 1b, n = 123) (t-test: p < .001,df = 231, t = −16.5).Hereafter, these lineages will be referred to as C3k 'max' and C3k 'min'.There were eight outlier samples that did not conform to the two main groups of symbionts, characterized by a high relative abundance of 'C1' and 'C1c' and found only in northern reefs, Osprey and Moore (Figure 1b).

| psbA and ITS2 marker alignment indicates presence/absence of novel symbiont species
Sequencing of the psbA ncr region allowed validation of ITS2 profiles as putative species (Smith et al., 2017).Here, a low entanglement score (0.012) highlighted the congruency between the ITS2 and psbA ncr gene regions for C. latusorum in P. meandrina/haplotype 8a and C. pacificum in P. verrucosa as previously reported (Turnham et al., 2021) (Figure 2a).meandrina and haplotype 8a samples (according to mtORF identification) were distributed across the C. latusorum north and south lineages, therefore, these samples were grouped for further analysis (hereafter: 'P. meandrina' sensu Johnston et al., 2022a).In contrast, there was a higher entanglement score between the ITS2 and psbA ncr markers of C3k hosted by A. cf humilis (entanglement score = 0.199), and the ps-bA ncr marker did not resolve differences between the two phylogenetic clusters of A. cf humilis for the ITS2 marker (Figure 2b).indicating a link between host genetic and symbiont cophylogeny (SOM, Figure S13).
Of the two genetic clusters detected within P. verrucosa with admixture analysis, (PVCL1 and PVCL2), PVCL1 was dominant across all reefs except Wreck reef, making up between 57%-96% of individuals per reef population.PVCL2 was detected in highest abundance from Wreck and Marion reefs, found in 61% and 31% of individuals, respectively (Figure 4b).We did not detect multiple lineages of the symbiont C. pacificum, which was ubiquitously found across the ten reef populations where P. verrucosa was sampled (Figure 4b), in contrast to patterns observed for symbionts hosted by P. meandrina.

| Biogeographical distribution of Acropora hosts and their symbionts
The three genetic clusters of A. cf humilis (referred to as AHCL1, AHCL2, and AHCL3) were confirmed by admixture analysis.While the dominant cluster AHCL1 was detected across each of the 11 sampled reefs, the range in relative abundance varied between 20% and 92% of sampled individuals per population.The two less common clusters were more centralized to specific reefs, where AHCL2 accounted for 30% of population structure within Marion reef.
The two symbiont lineages detected in A. cf humilis, C3k 'max' and 'min' were both found in all 11 reefs where Acropora was collected and did not appear to be associated with latitude or shelf position (Figure 4c).C3k 'max' was detected in the highest relative proportion of samples (78%-85%) in both low latitude reefs (Osprey, Bougainville, and Moore) and the highest latitude reefs (Wreck: 67% abundance).In contrast, C3k 'min' was found in highest proportion in central latitude reefs including Chilcott and Flinders reefs, detected in 93%-94% of samples (Figure 4c).

| Environmental and host genetic drivers of Pocillopora symbiont communities
A distance-based RDA evaluated symbiont community structure for each host species, in relation to environmental conditions and host genetic structure (Table 2).Both Pocillopora species had significant symbiont structuring attributed to host genetics, while accounting for the conditional effects of environmental variables (Figure 5a,b).
Overall, the strongest predictors of P. meandrina symbionts were driven by host genetic structure, explaining 14.9% of total model variance, compared to environment which explained 6.0% of total model variance (Figure 5a).In contrast, Pocillopora verrucosa symbionts had a lesser, but significant influence of host genetic structure, accounting for 3.4% of total variance, while all environmental factors accounted for 1.3% of explained variance (Figure 5b; Table 2).
verrucosa were significantly associated with host genetic clustering Note: The Environmental model incorporates all symbiont samples, while the Environment + Host Genetic model incorporates a reduced number of symbiont samples with matching DArT-sequencing host data.The type of environmental predictor and their relative strength are comparable between the two model types.(df = 1, F = 5.055, p = .001).Of all environmental and thermal metrics, depth was the only significant covariate of P. verrucosa symbiont community structure (df = 1, F = 2.142, p = .043),when accounting for host genetics as a conditional effect in the model (Table 2).Our data is consistent with the hypothesis that corals with horizontal transmission of symbionts have a greater capacity for adaptation via environmental sorting of symbionts adapted to the environment, compared to corals with vertical transmission driven by co-phylogeny (Chakravarti et al., 2017;van Oppen & Medina, 2020).

| Environmental and host genetic drivers of Acropora symbiont communities
The ability to validate such differences in Symbiodiniaceae and their respective environmental drivers can be attributed to the coanalysis of two genetic markers tested in this study.The congruence between the psbA ncr and ITS2 markers for Pocillopora symbionts, but not Acropora, supports the evolutionary stability of the symbiotic partnership in vertically transmitting coral species but not horizontally transmitting corals.There is evidence that this oceanographic current acts as a defining barrier to gene flow for organisms with larval dispersal (Payet et al., 2022), providing an ecological context to study rates of cophylogenetic speciation rates in organisms with host-symbiont diversification.

| Acropora symbiont communities are primarily structured by thermal history
The symbionts detected in A. cf humilis were nearly exclusive to the C3k radiation of Cladocopium, distinct from other C3-like species commonly found in Acropora within Australia's Great Barrier Reef and other Pacific islands (Butler et al., 2023).While symbionts were overall partitioned by thermal history, the effect of historical thermal stress (DHW between 1985 and 2020) explained twice as much variation as the concurrent bleaching event in 2020, suggesting that symbiont communities predominantly shift after, not during, a bleaching event (Silverstein et al., 2015).
The two main lineages of Acropora symbionts (C3k 'max' and 'min') associated with reefs specifically experiencing high and low historical heat stress, respectively.Thus, C3k 'max' may foster higher thermal tolerance given its detection in reefs which have experienced historically higher DHW and highlights a potential source of increased resilience for A. cf humilis in response to rapid environmental change.
The predominant effect of environment and minimal effect of host genetic structure has been observed in A. tenuis symbiont partitioning on the GBR, where communities were structured primarily by reef shelf position (Matias et al., 2023) and their proximity to freshwater plumes (Cooke et al., 2020), with a lesser effect of host genetic structure.Further, the conspicuous absence of Durusdinium, a common genus found in horizontally transmitting species such as Acropora in the GBR, highlights the niche requirements of Durusdinium to live in marginal, in-shore environments (Hoadley et al., 2019;LaJeunesse et al., 2014), leaving Cladocopium to dominate oligotrophic, offshore reefs such as the Coral Sea, explaining their detection in over 99.98% of samples quantified in this study.

| Effects of bleaching on immediate symbiont community restructuring varies by host taxa
The concurrent bleaching event during the collection period exposed corals to between 5. In the case of Pocillopora, P. meandrina symbiont communities were in alignment with host bleaching condition, but the communities harboured by P. verrucosa were not.Variable physiological limits of the host may explain the disparities in symbiont responses to the concurrent bleaching event.P. meandrina has demonstrated a 0.3°C lower thermal threshold in heat stress experiments compared to P.
verrucosa across nine of the reefs measured in this study (Marzonie et al., 2023).Further, Acropora exhibited the greatest sensitivity to heat stress evidenced by a 0.4-0.7°Clower thermal threshold compared to both Pocillopora spp.(Marzonie et al., 2023), corroborating the strong influence of multiple thermal disturbance and irradiance metrics (i.e., maxDHW, DHW2020, kd490, bleaching condition) for symbiont structure in Acropora observed here.The degree of thermal and light sensitivity of the host may influence the biogeography of symbiont communities, where lower tolerance to heat stress incurs higher bleaching and consequently a more pronounced symbiont community restructuring in response to recurrent thermal disturbances.However, due to the high collinearity between light attenuation and other environmental variables (e.g., PAR), we acknowledge there may be other environmental variables that contributed to, though are unaccounted for, as predictors of symbiont diversity for both Acropora and Pocillopora.

| CON CLUS ION
This study provides evidence that the diversity of coral-symbiont as- wrote the original draft of the manuscript.All authors provided critical edits to the manuscript.

ACK N OWLED G EM ENTS
We acknowledge the Meriam, Bindal, and Wulgurukaba Traditional Owners and their role as the first ecologists and custodians of the land and sea country where we conduct our work.Funding for the 2. Secondary loci were removed to retain a single SNP per sequencing tag at random.
3. SNPs were then filtered for reproducibility (0.98) and call rate (0.85) to confidently call by genotype.
4. Minor allele frequencies were filtered to remove rare allelic variants (MAF <0.05) and missing data was imputed.
5. Admixture models were run using the 'snmf' algorithm in the package LEA.First, cross-entropy criterion was plotted to estimate the optimal number of K ancestral populations to run.Ten repetitions were performed for each K, and the repetition with the lowest cross entropy was selected (snmf).
variation in symbiont communities by reef (Figure S7) and by host bleaching category (Figure S8) as the first signals of underlying environmental drivers relating to symbiont distribution.We then proceeded to test a suite of environmental variables predicted to influence symbiont communities given signals of reef-level structuring among PcoAs.To test additional environmental and host genetic covariates, ITS2 GUniFrac distances of the symbionts were incorporated into a constrained, distance-based Redundancy Analysis ('dbRDA', vegan).Thermal environmental variables were obtained for each reef from NOAA Coral Reef Watch Operational Daily Near-Real-Time Global 5-km Satellite Coral Bleaching Monitoring Product Version 3.1 from 1986 to 2020.As samples were collected during a bleaching event, a range of metrics were included to represent both climatic and environmental trends.The thermal metrics of the 'Environment' dbRDA model included accumulated Degree Heating Weeks (DHW) at the time of field collection (DHW2020), maximum historical DHW at each reef from 1986 to 2020 (maxDHW), and the bleaching condition of each coral individual at the time of collection (bleaching category).We then obtained other non-thermal variables from eReefs (Australian Institute of Marine Science) including kd490 (a measure of light attenuation to represent water turbidity), chlorophyll a, and secchi depth.Finally, we included spatial metrics collected during the expedition including collection depth, latitude, and longitude (SOM, Table

|
Characteristics of symbiont DIVs and type profiles among coral genera Sequencing of the ITS2 locus yielded 15,336,586 total sequences of Symbiodiniaceae across all host species and environments.The average per-sample depth of post-MED sequences was 52,523 reads (Acropora cf humilis: 50,979 reads; Pocillopora spp.: 53,761 reads).Over 99.98% of sequences were from the genus Cladocopium (15,335,016 sequences), with only 1552 sequences in the genus Symbiodinium and only 13 sequences from the genus Durusdinium (Table

F
I G U R E 1 Analysis of SymPortal post-MED ITS2 sequence data in Pocillopora (a) and Acropora (b).Each vertical bar represents the symbiont structure of one coral sample.The top third of each plot shows the UPGMA tree of the between-sample Generalized UniFrac distances.The middle third represents ITS2 sequence composition, with the proportion of individual sequences represented by different colours.The bottom third of each plot is coloured according to an assigned SymPortal ITS2 type profile, with grey representing the proportion of non-profile sequences found in each sample.For Pocillopora species (a), samples for P. verrucosa, P. meandrina and haplotype 8a are represented by tree leaf colours (green, orange, and blue, respectively) derived from the RFLP and mtORF molecular assays.Arrow 1 = C. pacificum/C.latusorum branch, 2 = C. latusorum 'north' and 'south' branch, and 3 = C3k 'max' and 'min' branch.[Colour figure can be viewed at wileyonlinelibrary.com] The two lineages belonging to C. latusorum in the ITS2-based UPGMA F I G U R E 2 Tanglegram of ITS2 and psbA ncr sequences.ITS2 sequences vs. psbA ncr sequence UPGMA trees for Pocillopora spp.(a) and Acropora cf humilis (b).ITS2 sequence UPGMA trees are according to k-mer = 7 + GUniFrac (0.5) distances and psbA ncr sequence UPGMA trees use k-mer = 7 distances.Blue and red arrows correspond to the phylogenetic split in each species tree shown in Figure 1a,b, respectively.Arrow 1 = C. pacificum/C.latusorum branch, 2 = C. latusorum 'north' and 'south' branch, and 3 = C3k 'max' and 'min' branch.[Colour figure can be viewed at wileyonlinelibrary.com] tree of GUniFrac distances (Figure 2a; C. latusorum 'north' and C. latusorum 'south') were completely congruent with the psbA ncr marker.P.
Admixture, PcoAs, and neighbour-joining trees revealed two host genetic clusters for P. meandrina, two clusters for P. verrucosa, and three clusters for A. cf humilis.A Procrustean Rotation Analysis examined the relationship between the population structure of each coral host species and their respective symbiont communities, revealing a significant degree of cophylogeny for each host-symbiont partnership, the strength of which varied among host taxa.P. meandrina demonstrated the highest degree of cophylogeny, indicated by a strong, positive procrustean correlation (n = 66, Sum of Squares = 0.724, r = .52,p < .001)forming two distinct clusters of host plus symbiont (indicated by short distances between two given points of the same sample, Figure3a).In contrast, the degree of cophylogeny was weaker, yet significant, for P. verrucosa (n = 130, Sum of Squares = 0.882, r = .34,p = .02),characterized by one primary cluster and less partitioning within the ordination (Figure3b).Patterns of cophylogeny in A. cf humilis were significant when all three host genetic clusters were included (n = 257, Sum of Squares = 0.954, r = .21,p = .01)(Figure3c).Symbionts were less organized than their respective host and formed one central cluster (SOM, FigureS12).A reduced Procrustes Rotation Analysis with only the numerically dominant host genetic cluster (AHCL1) did not yield a significant result(n = 201, Sum of Squares = 0.948, r = .22,p = .14) patterns of host and symbiont structuring among the three coral taxa.Pocillopora meandrina was more abundantly sampled in high latitude reefs of the Coral Sea, while P. verrucosa was sampled abundantly in lower latitude reefs, likely reflecting biogeographic patterns in species distribution (Figure 4a,b).Pocillopora meandrina comprised 77%-100% of samples collected in the three highest latitude reefs positioned at 21-22° and was sampled in high abundance within low latitude, most offshore reefs, making up 72% and 81% of samples collected from Moore and Chilcott reefs, respectively (Figure 4a).In contrast, P. verrucosa made up 74% of samples collected from the lowest latitude reef, Osprey, as well as 67%-77% F I G U R E 3 Procrustean rotation analysis of host and symbiont genetic diversity.Pocillopora verrucosa (a), Pocillopora meandrina (b), and Acropora cf humilis (c).Each point represents a sample, with host samples (triangle) and symbiont samples (circle) from the same individual linked by a line segment.Each point is coloured by reef.Inset shows symbiont points without host points overlaid to indicate symbiont distribution by reef.[Colour figure can be viewed at wileyonlinelibrary.com] of samples collected from mid-latitude reefs in closest proximity to Australia's coast (Holmes and Flinders) (Figure 4b).Intraspecific trends in host population structure were detected, with two sympatric genetic clusters per species of Pocillopora.For P. meandrina, the two genetic clusters confirmed by admixture analysis (PMCL1 and PMCL2) were divergent across a latitudinal gradient (Figure 4a), with PMCL2 accounting for 55% (range: 35%-71%) of population structure in low latitude populations, compared to PMCL1, comprising 89% (range: 76%-96%) of individuals in high latitude reefs.The two symbiont lineages, C. latusorum 'north' and 'south' mirrored trends of host genetic structure, where C. latusorum 'north' was detected in reefs where the host cluster PMCL2 was found.In three of the highest latitude reefs, corals hosted 100% C.

F
Schematic diagram depicting biogeographic patterns and relative proportion of host population structure and symbiont lineages detected throughout the Coral Sea Marine Park at 13 reefs.(a) Pie charts represent host genetic clustering (admixture proportion) for P. meandrina (top) and symbiont lineage clustering (bottom) include C. latusorum 'north' and C. latusorum 'south'.(b) Host genetic clustering for P. verrucosa (top) and one symbiont lineage of C. pacificum (bottom).(c) Pie charts show relative abundance of three host genetic clusters of A. cf humilis (top) and symbiont lineage clustering (bottom) include C3k 'max' and C3k 'min'.[Colour figure can be viewed at wileyonlinelibrary.com] Distance-based redundancy analysis (dbRDA) ordinations of GUniFrac (0.5) ITS2 sequence composition.Pocillopora meandrina (a), Pocillopora verrucosa (b) and Acropora cf humilis (c).Each circle is coloured by latitude.Blue arrows represent significant environmental vectors derived from a backward stepwise model.[Colour figure can be viewed at wileyonlinelibrary.com] Communities of C3k symbionts hosted by A. cf humilis were structured predominantly by a suite of environmental factors, accounting for 15% of total model variance, compared to only 0.2% of total model variance explained by the host genetic model.For Acropora symbionts, the significant environmental predictors included historical thermal stress (maxDHW; df = 1, F = 15.198,p = .001),Degree Heating Weeks at the time of collection (DHW2020; df = 1, F = 6.046, p = .001),light attenuation (kd490; df = 1, F = 6.829, p = .001),latitude (df = 1, F = 5.042, p = .001),depth (df = 1, F = 3.009, p = .015),and host bleaching condition (df = 1, F = 2.377, p = .048)(Figure5c).Additionally, a separate model was run using only the dominant host genetic cluster of A. cf humilis (Host cluster: AHCL1, n = 198).Comparable to the full Acropora model, AHCL1 symbionts were consistently driven by the same four environmental predictors, except for bleaching category of the host, which was only significant in the full model but not the AHCL1 model (p = .116),indicating an interaction between host genotype and bleaching condition.4 | DISCUSS ION 4.1 | Host-specific trends of symbiont community structure among Pocillopora and Acropora Symbiodiniaceae community structure can be driven by a suite of host and environmental factors and varies by the host-symbiont partnership.Here, symbiont specificity was most strongly driven by host genera and their mode of symbiont transmission (i.e., Acropora and Pocillopora), with concomitant impacts of the environmental factors shaping their distribution.Both Pocillopora species maintained symbionts congruent with host genetic structure, 2.5-fold higher than structuring by their environment.The host genetic structure of P. meandrina was considerably stronger in partitioning symbiont communities than P. verrucosa, demonstrated by 4-fold higher environment and host genetic influence, as well as high alignment between the two host genetic clusters and two symbiont clusters.In contrast, symbiont communities in A. cf humilis maintained an inverse trend, structured 75 times more by environmental predictors than by host genetics.Specifically, symbionts were partitioned by thermal history experienced at each reef.Both DHW experienced at the time of coral collection (DHW2020), and historical thermal stress patterns over 35 years (maxDHW) were influential, indicating the importance of historical and contemporary thermal stress in structuring symbiont communities in Acropora.Similarly, the light attenuation coefficient (kd490) was a stronger factor in driving Acropora symbionts, compared to either Pocillopora spp., further alluding to the higher symbiont structuring in Acropora in response to both thermal and environmental conditions.

4. 2 |
Host genetic structure and cophylogeny as key traits for symbiont communities in PocilloporaInTurnham et al. (2021), regionally distinct lineages of C. latusorum were presented across the Indo-Pacific.Here we show that C. latusorum lineages in P. meandrina arise at much smaller spatial scales (<2000 km), and with multi-gene congruence, these regionally differentiated lineages may warrant descriptions as mutually exclusive species.Further, the stronger co-phylogeny in P. meandrina and symbionts indicates that host population structure attributes to the species-specific differences in symbiont community structure and the degree of host influence.The geographic boundaries of Pocillopora have previously been shown to vary among species, with P. verrucosa localized to low latitude reefs and P. meandrina to high latitude reefs(Johnston et al., 2018), similar to the species distribution ranges of Pocillopora in this study.A myriad of biological factors can shape species ranges, including depth distribution and the rate of larval development.In addition, abiotic factors including habitat heterogeneity or geological features mediate a species' ability to occupy space over a geographic break(Keith et al., 2013), possibly explaining the two distinct host and symbiont groups in high and low latitude populations observed in P. meandrina.The two host genetic clusters and two lineages of C. latusorum correspond to the strong bifurcation in the South Equatorial Current at 16° S(Ceccarelli et al., 2013).
7 and 10.0 DHW, dependent on the sampling location.Correlations between the extent of bleaching and the symbiont communities varied across the three host taxa and with respect to the accumulated heat stress in the environment.Acropora was the most thermally sensitive species, and it is thus consistent with Acropora symbiont community structure being predominantly influenced by all thermal and irradiance parameters in the model, including DHW at the time of collection, historical thermal stress (maximum DHW between 1986 and 2020), and light attenuation.Host bleaching condition strongly dictated symbiont structure in A. cf humilis when accounting for the three host genetic clusters but not for the main host genetic cluster alone (AHCL1), indicating an interaction between host genetic lineages and bleaching condition on symbiont community structure.InQuigley et al. (2022), repeat mass bleaching events resulted in variable levels of symbiont restructuring among three species of Acropora in the Great Barrier Reef, highlighting a host taxa-specific effect of bleaching on symbionts.We found similar patterns of symbiont structuring at the host population level, applying not only to species-specific, but possibly to the level of host populations or individual genotypes(Rose et al., 2021).
semblages is influenced by a range of host-specific and environmental factors.The ability of the coral-symbiont partnership to adjust to fluctuating environmental conditions may rely on a suite of both host factors (i.e., host co-phylogeny, symbiont mode of transmission, population structure) and the environment (i.e., thermal history, climatology, and irradiance).Contextualizing both host and symbiont genetic structure in response to environmental drivers is critical in the search for the mechanisms of coral adaptation.As climate change continues to impact coral reefs in the coming decades, there will undoubtedly be shifts in symbiont community structure.Tracking these changes is essential to quantifying coral species and populations with both high and low adaptive potential against a backdrop of rapid ocean warming and should be done with a combination of molecular and environmental data for both hosts and their symbionts.AUTH O R CO NTR I B UTI O N S M.R.M., H.B.H. and L.K.B. designed the experimental collection.M.R.M. and H.B.H. conducted the field work.M.R.M. conducted the laboratory work.M.R.N. and M.R.M. analysed the data.M.R.M.