Phillip A Morin

Brenna R Forester

Karin A Forney

Carla A Crossman

Brittany L Hancock-Hanser

Kelly M Robertson

Lance G Barrett-Lennard

Robin W Baird

John Calambokidis

Pat Gearin

M Bradley Hanson

Cassie Schumacher

Timothy Harkins

Michael C Fontaine

Barbara L Taylor

Kim M Parsons

Population structure in a continuously distributed coastal marine species, the harbor porpoise, based on microhaplotypes derived from poor quality samples

Keywords: sPCA, dbRDA, seascape genetics, mtDNA, SNP, microhaplotype, GT-seq, cetacean

ou non, émanant des établissements d'enseignement et de recherche français ou étrangers, des laboratoires publics ou privés.

Introduction

Harbor porpoises (Phocoena phocoena) are small cetaceans found in temperate and sub-Arctic coastal waters, typically in less than 200m depth [START_REF] Fontaine | Harbour porpoises, Phocoena phocoena, in the Mediterranean Sea and adjacent regions: Biogeographic relicts of the last glacial period[END_REF][START_REF] Read | Harbour porpoise (Phocoena phocoena)[END_REF] (but see [START_REF] Nielsen | Oceanic movements, site fidelity and deep diving in harbour porpoises from Greenland show limited similarities to animals from the North Sea[END_REF]. Their range is considered continuous along the continental coasts, but there is high variability in density and apparent gaps in suitable habitat [START_REF] Evenson | Disappearance and return of harbor porpoise to Puget Sound: 20 year pattern revealed from winter aerial surveys[END_REF]Forney, Moore, Barlow, Carretta, & Benson, in press). They are susceptible to entanglement in gillnets and disturbance from construction activities for wind farms [START_REF] Carstensen | Impacts of offshore wind farm construction on harbour porpoises: acoustic monitoring of echolocation activity using porpoise detectors (T-PODs)[END_REF][START_REF] Reeves | Marine mammal bycatch in gillnet and other entangling net fisheries, 1990 to 2011[END_REF]. In addition, the linear coastal distributions with gaps raises the question of how harbor porpoise will respond to climate change. Areas with sufficient data indicate fine-scaled population structure is common [START_REF] Chivers | Small-scale population structure of eastern North Pacific harbour porpoises (Phocoena phocoena) indicated by molecular genetic analyses[END_REF][START_REF] Crossman | Population structure and intergeneric hybridization in harbour porpoises Phocoena phocoena in British Columbia, Canada[END_REF][START_REF] Fontaine | Harbour porpoises, Phocoena phocoena, in the Mediterranean Sea and adjacent regions: Biogeographic relicts of the last glacial period[END_REF][START_REF] Lah | Spatially explicit analysis of genome-wide snps detects subtle population structure in a mobile marine mammal, the harbor porpoise[END_REF][START_REF] Rosel | Genetic structure of harbour porpoise Phocoena phocoena populations in the northwest Atlantic based on mitochondrial and nuclear markers[END_REF][START_REF] Tiedemann | Mitochondrial DNA sequence patterns of harbour porpoises (Phocoena phocoena) from the North and the Baltic sea[END_REF][START_REF] Walton | Population structure of harbour porpoises Phocoena phocoena in the seas around the UK and adjacent waters[END_REF][START_REF] Wang | Mitochondrial DNA analysis of harbour porpoises (Phocoena phocoena) in the Baltic Sea, the Kattegat-Skagerrak Seas and off the west coast of Norway[END_REF].

Along the U.S. west coast, early genetic studies based on mitochondrial DNA (mtDNA) provided data to delineate management stocks [START_REF] Chivers | Small-scale population structure of eastern North Pacific harbour porpoises (Phocoena phocoena) indicated by molecular genetic analyses[END_REF], but data gaps made boundary placement difficult. Because harbor porpoises avoid vessels, which makes dartbiopsy an unviable method to obtain sufficient samples, most samples are from either beach-stranded animals or entangled animals from areas with gillnet fisheries. As a result, sampling gaps existed in areas without gillnet fisheries and samples were often obtained days after death and in moderate to advanced stages of decomposition, and hence had degraded DNA.

Although gillnetting has declined along the U.S. west coast, mitigation of other potential threats, such as development of wind farms, requires an understanding of population structure, often in areas poorly sampled in previous studies. The U.S. Marine Mammal Protection Act (MMPA) manages at the scale of demographically independent populations (DIPs) where allele (or mtDNA haplotype) frequencies will differ but evolutionary differences are not expected. Because our area of interest includes the southernmost part of harbor porpoise distribution in the Pacific where ocean temperatures are already rising, an understanding of evolutionary barriers to dispersal is also of interest to evaluate potential impacts of such temperature shifts.

In our study area along the U.S. West Coast and inland waters of Washington, samples have slowly accrued to fill most gaps over a period of 30 years, but most samples are of degraded quality and sample size remains low in some areas. An additional complication requiring consideration for marker choice is that there is evidence of intergeneric hybridization with Dall's porpoise (Phocoenoides dalli) in the eastern North Pacific [START_REF] Baird | An intergeneric hybrid in the family Phocoenidae[END_REF][START_REF] Crossman | Population structure and intergeneric hybridization in harbour porpoises Phocoena phocoena in British Columbia, Canada[END_REF][START_REF] Willis | Natural hybridization between Dall's porpoises (Phocoenoides dalli) and harbour porpoises (Phocoena phocoena)[END_REF]), which could affect population analyses, especially if hybridization has been regionally restricted and/or resulted in introgressive gene flow between species.

Genomic methods such as genome re-sequencing (e.g., [START_REF] Foote | Killer whale genomes reveal a complex history of recurrent admixture and vicariance[END_REF] and reduced representation sequencing (e.g., [START_REF] Andrews | Harnessing the power of RADseq for ecological and evolutionary genomics[END_REF][START_REF] Maisano Delser | Population genomics of C. melanopterus using target gene capture data: demographic inferences and conservation perspectives[END_REF] can produce thousands to tens of thousands of genetic markers, providing unprecedented power to detect population differences (e.g., [START_REF] Candy | Population differentiation determined from putative neutral and divergent adaptive genetic markers in Eulachon (Thaleichthys pacificus, Osmeridae), an anadromous Pacific smelt[END_REF][START_REF] Emerson | Resolving postglacial phylogeography using high-throughput sequencing[END_REF][START_REF] Leslie | Using genome-wide SNPs to detect structure in high-diversity and low-divergence populations of severely impacted eastern tropical Pacific spinner (Stenella longirostris) and pantropical spotted dolphins (S. Attenuata)[END_REF] and identify candidate loci under selection [START_REF] Ahrens | The search for loci under selection: trends, biases and progress[END_REF]. However, for studies restricted to the use of poor quality archived samples, the optimal strategy involves targeting a reduced number of genetic markers that can efficiently and reproducibly be obtained from a large number of samples, such as "genotyping in thousands by sequencing" (GT-seq; [START_REF] Campbell | Genotyping-in-Thousands by sequencing (GT-seq): A cost effective SNP genotyping method based on custom amplicon sequencing[END_REF]. Typically targeting a few hundred loci, GT-seq relies on multiplexed short amplicon sequencing that requires only ~20ng of DNA, and locus variability can be increased by targeting highly variable regions with multiple single nucleotide polymorphisms (SNPs) genotyped as microhaplotypes [START_REF] Baetscher | Microhaplotypes provide increased power from short-read DNA sequences for relationship inference[END_REF][START_REF] Mckinney | Managing mixed-stock fisheries: genotyping multi-SNP haplotypes increases power for genetic stock identification[END_REF].

To increase power to detect population structure from larger numbers of SNPs, while maximizing our ability to use an existing collection of poor-quality samples collected over three decades, we employed GT-seq multiplex sequencing combined with microhaplotype analysis, and traditional mtDNA control region sequencing to incorporate new and previously published data. We describe population structure at both evolutionary and demographic scales through a continuous range in the eastern North Pacific by applying a combination of Bayesian assignment to detect more substantial divergence of evolutionarily significant units (ESUs; [START_REF] Moritz | Defining 'Evolutionarily Significant Units' for conservation[END_REF], and traditional divergence metrics and spatially explicit methods to detect demographically independent populations (DIPs), the basis of marine mammal population stocks under the U.S. Marine Mammal Protection Act [START_REF] Martien | The dip delineation handbook: A guide to using multiple lines of evidence to delineate demographically independent populations of marine mammals[END_REF][START_REF] Waples | What is a population? An empirical evaluation of some genetic methods for identifying the number of gene pools and their degree of connectivity[END_REF]. Evolutionarily significant units are expected to have very low gene flow (on the order of a few successful dispersers per generation) and significant differences in both nuclear and mitochondrial markers. In contrast, a demographic independence means that the population dynamics of the affected group is more a consequence of births and deaths within the group (internal dynamics) than immigration or emigration (external dynamics) [START_REF] Martien | The dip delineation handbook: A guide to using multiple lines of evidence to delineate demographically independent populations of marine mammals[END_REF]. DIPs could have gene flow on the order of 1-2% per year. Besides the large difference in gene flow, DIPs could be based solely on mtDNA because birth and death rates depend on females, though use of larger numbers of nuclear markers may provide higher statistical power to detect population differences.

We also use seascape genotype-environment association to investigate patterns of localized adaptation. We anticipate that these genomic scale data will not only provide increased resolution to detect spatial structure among populations, but also allow us to link genotypes to environmental variables that vary both spatially and temporally. This emerging field of seascape genomics [START_REF] Riginos | Navigating the currents of seascape genomics: how spatial analyses can augment population genomic studies[END_REF] is still in its early stages relative to landscape genomics, due to both the difficulty of broad scale sampling of many marine species and the limited availability of relevant environmental predictors. However, the increasing accessibility of remotely-sensed oceanographic variables has improved our ability to identify spatial, temporal, and ecological factors that promote population structure and local adaptation in complex and dynamic seascape environments (reviewed by [START_REF] Riginos | Navigating the currents of seascape genomics: how spatial analyses can augment population genomic studies[END_REF][START_REF] Selkoe | A decade of seascape genetics: contributions to basic and applied marine connectivity[END_REF].

Materials and Methods

Samples

Harbor porpoise (N=441), Dall's porpoise (N=9) and putative Harbor/Dall's porpoise hybrids (N=13) skin samples were collected from beach-cast carcasses, carcasses recovered as fisheries bycatch, or from animals live-captured for tagging. Tissue samples were preserved in salt-saturated 20% DMSO or 100% ethanol and subsequently stored at -20°C in the U.S. National Marine Fisheries Service (NMFS) Marine Mammal and Sea Turtle Research (MMASTR) Collection at the Southwest Fisheries Science Center (SWFSC). Sample information is in supplemental Table S1 and locations are shown in Figure 1. DNA was extracted from tissue samples using a variety of common extraction methods, including silica-based filter membranes (Qiaxtractor® DX reagents, Qiagen, Valencia, CA, USA), standard phenol/chloroform extraction (modified from [START_REF] Sambrook | Molecular Cloning: A Laboratory Manual[END_REF], lithium chloride [START_REF] Gemmell | An efficient method for the extraction of DNA from vertebrate tissues[END_REF] and sodium chloride protein precipitation [START_REF] Miller | A simple salting out procedure for extracting DNA from human nucleated cells[END_REF], and quantified using the Quant-iT PicoGreen dsDNA assay kit (Invitrogen, Paisley, UK) with a Victor X3 fluorospectrometer (Perkin Elmer, Waltham, MA, USA).

mtDNA control region sequencing

New 395bp mitochondrial DNA control region sequences (N=176) were generated according to previously published methods for harbor porpoise studies [START_REF] Chivers | Small-scale population structure of eastern North Pacific harbour porpoises (Phocoena phocoena) indicated by molecular genetic analyses[END_REF][START_REF] Chivers | Additional genetic evidence for population structure of Phocoena phocoena off the coasts of California[END_REF][START_REF] Crossman | Population structure and intergeneric hybridization in harbour porpoises Phocoena phocoena in British Columbia, Canada[END_REF]. Electropherograms of sequences from previous studies (N=224 from [START_REF] Chivers | Additional genetic evidence for population structure of Phocoena phocoena off the coasts of California[END_REF]; N=91 from [START_REF] Crossman | Population structure and intergeneric hybridization in harbour porpoises Phocoena phocoena in British Columbia, Canada[END_REF] were evaluated and compared to newly generated sequences, and regions with ambiguities or poor quality were either re-evaluated by a single person to ensure consistent interpretation (N=224), or the DNA was sequenced again from new PCR reactions using current Sanger sequencing chemistry (N=91). Haplotype IDs were assigned using the LabelHaplotypes function in the R package strataG (v2.4.905;[START_REF] Archer | STRATAG: An R package for manipulating, summarizing and analysing population genetic data[END_REF] and associated with previously published sequences and haplotypes (supplemental Table S2).

SNP discovery

DNA from 12 North Pacific harbor porpoises was used for indexed genomic library preparation using the Accel-NGS 2S PCR-free genomic library preparation kit (Swift Biosciences, Ann Arbor, MI, USA). Either 100 ng or 500 ng of genomic DNA, determined by Qubit (Thermo Fisher Scientific, Waltham, MA, USA) quantification, was used for the library prep to obtain sufficient library product to pool 15 nmol/l of each sample prior to sequencing. For 12 North Atlantic samples, genomic DNA was used for indexed genomic library preparation at BGI (Hong Kong) following their proprietary protocol and pooled for next generation sequencing (NGS). The North Pacific pooled library was first sequenced in two 150 bp, paired-end Illumina, Inc. (La Jolla, CA, USA) MiSeq lanes, then both pooled libraries were sequenced in two Illumina HiSeq-4000 lanes each. NGS read data were trimmed and filtered using Trimmomatic (v. 0.36;[START_REF] Bolger | Trimmomatic: a flexible trimmer for Illumina sequence data[END_REF] to remove Illumina adapters, reads <50bp and low-quality bases (Q<20) at the beginning and end of reads. The sliding window approach was implemented to change internal bases with Q<15 (sliding window of 4bp) to N's. Because there was some bias in the nucleotide frequencies, the first 4 bp of all sequences were also removed.

Genome alignment and SNP discovery were conducted as previously described [START_REF] Morin | SNP discovery from single and multiplex genome assemblies of nonmodel organisms[END_REF]. The repeat-masked killer whale genome (accession GCA_000331955.2; [START_REF] Foote | Convergent evolution of the genomes of marine mammals[END_REF] was used as the reference for genome assembly. Briefly, the paired-end MiSeq reads from the 12 North Pacific samples were assembled de novo using CLC Genomics Workbench (v.4.1; CLCbio) to obtain a complete reference mitochondrial genome as previously described [START_REF] Hancock-Hanser | Targeted multiplex next-generation sequencing: Advances in techniques of mitochondrial and nuclear DNA sequencing for population genomics[END_REF]. The NGS reads from 12 samples from each ocean basin subspecies were aligned separately to the harbor porpoise mitochondrial genome using BWA mem (v. 0.7.5a; [START_REF] Li | Fast and accurate short read alignment with Burrows-Wheeler transform[END_REF] and non-aligned reads were extracted using samtools (v. 1.2; Li et al., 2009). The extracted (nonmitochondrial) reads from each sample were separately aligned to the killer whale reference genome, combined into one alignment, and the consensus harbor porpoise genome sequence generated. The nuclear DNA reads for each sample were aligned to the new harbor porpoise consensus genome sequence as above, followed by SNP discovery separately from each subspecies using GATK (v. 2.5-2;[START_REF] Depristo | A framework for variation discovery and genotyping using next-generation DNA sequencing data[END_REF][START_REF] Mckenna | The Genome Analysis Toolkit: A MapReduce framework for analyzing next-generation DNA sequencing data[END_REF]. Potential SNPs were filtered to remove SNPs with mapping quality <30, excessive coverage (>2x mean depth of coverage), and estimated minor allele frequency <0.05. To avoid linked loci, SNPs were selected from contigs that were at least 100 kb in length, and SNPs were at least 100 kb apart on the contigs. SNP loci from the two subspecies were compared to identify loci that were polymorphic in both ocean basins to avoid ascertainment bias in application to either subspecies. Finally, candidate SNPs were compared to GenBank using BLAST+ [START_REF] Camacho | BLAST+: architecture and applications[END_REF] and filtered to remove loci that were potentially in repeat regions, gene families, or a close match to non-mammalian species. The resulting set of filtered SNPs are subsequently referred to as "targeted SNPs", as these were the initial targets for our SNP genotyping effort.

Additional SNPs detected in short GT-seq sequences were combined with the targeted SNPs to form microhaplotype genotypes.

Multiplex primer design

Primers were designed from a batch of 500 loci selected randomly from the filtered SNPs using the program FastPCR [START_REF] Kalendar | Java web tools for PCR, in silico PCR, and oligonucleotide assembly and analysis[END_REF]. Parameters for primer selection (from 300 bp sequences with the target SNP at position 151) were: primer length = 15-32, tm = 57 -62, 3' Tm = 25-50, dimer stringency = 5, synchronized Tm for primer pair = 5, Forward primers between position 40 and 150, Reverse primers between position 152 and the 3' end, and addition of 5' GT-seq tails for indexing and library preparation: F-tail = CGACAGGTTCAGAGTTCTACAGTCCGACGATC, R-tail = GTGACTGGAGTTCAGACGTGTGCTCTTCCGATCT. After primer design, all loci with both forward and reverse primers were compared using the "Primers list analysis" function in FastPCR to detect cross-locus primer dimer interactions with Tm > 20°C ("strong" primer dimers). Loci were filtered out of the primer list if either one or both primers had predicted interactions with Tm > 40°C with >2 other primers, or if the predicted primer interactions were >46°C. Primers for 385 loci were synthesized at 100 µM concentration in 96-well plates by Integrated DNA Technologies (Coralville, Iowa, USA).

Multiplex PCR optimization

GT-seq primers were pooled and used for multiplex amplification of one sample initially to optimize the locus set prior to genotyping. Optimization consisted of multiple rounds of GT-seq library preparation as described by [START_REF] Campbell | Genotyping-in-Thousands by sequencing (GT-seq): A cost effective SNP genotyping method based on custom amplicon sequencing[END_REF], and sequencing a small portion of the library (e.g., 1-10 million reads) to determine the relative abundance of reads per locus, and presence of primer artifacts as determined by the published analysis scripts (https://github.com/GT-seq/GT-seq-Pipeline). Loci were removed at each iteration to eliminate loci represented by disproportionately high read depths, evidence of primer artifacts, or low ratios of the probe to primer target sequences. There are no published guidelines for cut-off values, so we removed loci that appeared to be outliers for any of these values, and based on expert advice from experienced users of the GT-seq method (see acknowledgements). Primer sequences for the final set of 292 loci used for genotyping are in supplemental Table S3.

SNP genotyping

Amplicon libraries were prepared following the GT-seq protocol, including the optional Exo-SAP pre-treatment of the samples [START_REF] Campbell | Genotyping-in-Thousands by sequencing (GT-seq): A cost effective SNP genotyping method based on custom amplicon sequencing[END_REF], and pooled libraries were sequenced on an Illumina NextSeq500 sequencer, 1x150 bp reads. Custom scripts for processing GT-seq data [START_REF] Campbell | Genotyping-in-Thousands by sequencing (GT-seq): A cost effective SNP genotyping method based on custom amplicon sequencing[END_REF] were used to demultiplex the sample files and conduct preliminary genotyping. Genotypes were quality checked for duplicate samples, percent missing genotypes per locus and sample, and percent homozygosity using the strataG package in R. Replicate samples were used to estimate genotyping error rates, then fastq files from replicates samples were combined to a single file per sample.

Fastq files were checked for standard quality metrics (e.g., per base quality scores, nucleotide composition, sequence duplication level, overrepresented sequences) with FASTQC v0.11.3 (Babraham Bioinformatics), then trimmed using FASTP [START_REF] Chen | fastp: an ultra-fast all-in-one FASTQ preprocessor[END_REF] to remove adapter sequences and poly-A and poly-G 3' tails that were added during sequencing of amplicons shorter than 150 bp, and to exclude reads shorter than 30bp after trimming. Reads were mapped to the reference locus sequences using the BWA MEM algorithm (v. 0.7.15;[START_REF] Li | Fast and accurate short read alignment with Burrows-Wheeler transform[END_REF], and SNPs detected across all samples using FREEBAYES v1.1.0-54-g49413aa [START_REF] Garrison | Haplotype-based variant detection from short-read sequencing[END_REF] after removing sample files smaller than 1 MB (containing <0.1 M filtered reads).

FREEBAYES was run with minimal filtering (supplemental materials), followed by additional filtering with vcfTools v0.1.12b [START_REF] Danecek | The variant call format and VCFtools[END_REF] to extract the targeted SNP for each locus (for targeted SNP analysis; see below), and to remove sites with low coverage (minimum depth = 10), indels, and loci with less than 30% completed genotypes.

For targeted SNPs (position 151 in all reference sequences), custom scripts in R (supplemental materials) were used to extract the genotype data from the vcf file, generate allelic count plots to visualize the genotype distributions of reads for each allele, and re-call genotypes based on minimum depth and allelic ratios. The minimum depth of 10 reads total and default minor allele read proportion for heterozygotes of >0.3 were adjusted as needed until genotypes clearly fell into separate clusters in the allelic plots.

Loci with poor resolution of plotted genotypes were removed from the data set.

Microhaplotypes (containing the targeted SNPs and/or newly discovered SNPs) were generated for all loci using the R package MicrohaPlot [START_REF] Baetscher | Microhaplotypes provide increased power from short-read DNA sequences for relationship inference[END_REF]. The MicrohaPlot algorithm inserts N's for missing sequence data at SNPs within haplotypes, so we used a custom R-scripts (supplemental materials) to identify SNPs with >10% N's.

The identified SNPs were removed from the original vcf file using vcfTools, and MicrohaPlot was used to generate new microhaplotypes with the remaining variable SNP positions. The unfiltered haplotypes were exported for subsequent filtering with custom scripts to view and call genotypes similar to the methods described above for targeted SNPs (supplemental materials). The few remaining microhaplotypes with N's in them were excluded from genotypes prior to analysis.

A final combined data set for microhaplotypes and targeted SNPs was created by combining the multi-SNP loci with the single-SNP (targeted) loci. Since the targeted SNP loci had been genotyped using two different methods, we selected the genotype data for each locus from the method that provided the higher quality or quantity of genotypes for the targeted SNP. Some microhaplotype loci that were monomorphic in harbor porpoises, or which had similar genotype quality for the targeted SNP but contained other SNPs present in Dall's porpoise samples, were retained (instead of the targeted locus data) to allow genetic identification of intergeneric hybrids. All loci that were found only in the microhaplotype or targeted SNP data set were then added to the filtered loci to generate a final data set.

Quality analysis

Quality analysis and sample and locus filtering were conducted using custom R scripts.

Samples missing >80% of the genotypes, and loci missing more than 45% of the genotypes were also removed. Genetic duplicates (>80% identity) were identified and one from each pair of samples identified as duplicates was removed. Analyses of deviations from Hardy-Weinberg equilibrium (HWE) expectations were conducted across all samples using the R package adegenet, and loci with a difference between observed and expected heterozygosity >0.2 were removed as extreme outliers (>10x the average difference of 0.02), most likely due to non-mendelian loci (e.g., null alleles, duplicated loci, or high error rates). Remaining loci were tested for significant linkage disequilibrium (LD) and deviations from HWE within three discrete, geographicallydefined strata represented by greater than 20 samples (inland waterways (N=88), Neah Bay (N=21), and Northern California/Southern Oregon (N=35)), after correction for multiple tests using a sequential correction [START_REF] Holm | A simple sequentially rejective multiple test procedure[END_REF]. We used the R package Demerelate (v. 0.9.3; Kraemer & Gerlach, 2017) using relatedness estimators "Wang" and "Mxy" to test for inadvertent sampling of close relatives.

Porpoise distribution data

To examine the genetic results in the context of harbor porpoise distribution and relative density along the U.S. West Coast, independent aerial survey data collected during 1991-2017 off California (e.g., [START_REF] Forney | Detecting trends in harbor porpoise abundance from aerial surveys using analysis of covariance[END_REF]), and during 1989[START_REF] Storey | Statistical significance for genomewide studies[END_REF]off Oregon and Washington (e.g., Calambokidis, Laake, & Osmek, 1997) were processed to derive the number of porpoise seen per kilometer surveyed as an index of relative density. Although the transect design differed between these two data sets, the survey protocols, observer team and configuration, and aircraft type were the same. Survey data were truncated spatially to include only the primary porpoise habitat extending from shore out to 90-100 m water depth. Relative densities were calculated for each transect line, assigned to latitude of the transect mid-point, and then smoothed south-to-north using a Loess smoother.

Habitat data

A variety of modeled and measured oceanographic variables (Table 1) were extracted to examine potential environmental correlates of genetic patterns. These predictors included 1) sea surface temperature, sea surface height, mixed layer depth, and the standard deviation of these three variables derived from the Regional Ocean Modeling System (ROMS) outputs [START_REF] Moore | The Regional Ocean Modeling System (ROMS) 4dimensional variational data assimilation systems Part II -Performance and application to the California Current System[END_REF], 2) coastal upwelling indices [START_REF] Jacox | Coastal upwelling revisited: Ekman, Bakun, and improved upwelling indices for the US west coast[END_REF], and 3) multispectral ultra-high resolution sea surface temperature and its standard deviation [START_REF] Chin | A multi-scale highresolution analysis of global sea surface temperature[END_REF].

Analytical methods

Several methods were used to estimate the number of populations and population assignment based on the genetic data. Population structure and individual assignment was examined using STRUCTURE (v. 2.3), which implements a Bayesian clustering method to identify significant genetic clusters based on Hardy-Weinberg equilibrium allele frequency expectations [START_REF] Hubisz | Inferring weak population structure with the assistance of sample group information[END_REF][START_REF] Pritchard | Inference of population structure using multilocus genotype data[END_REF]. We ran ten replicates for each value of k (where k is the number of putative populations), using correlated allele frequencies and an admixture model, with location prior using geographically defined units (supplemental Table S1).

Each analysis consisted of 50,000 burn-in steps followed by 100,000 MCMC steps, and 10 replicates combined using 100 iterations in CLUMPP (v. 1.1.2; [START_REF] Jakobsson | CLUMPP: a cluster matching and permutation program for dealing with label switching and multimodality in analysis of population structure[END_REF]. We also used CLUMPAK [START_REF] Kopelman | Clumpak: a program for identifying clustering modes and packaging population structure inferences across K[END_REF] to assess convergence of the MCMCs and evaluate consistency of replicates across values of K. The ΔK method [START_REF] Evanno | Detecting the number of clusters of individuals using the software STRUCTURE: a simulation study[END_REF] was used to evaluate most likely number of inferred clusters.

STRUCTURE is known to have low power to detect populations (i.e., DIPs) when there is even a very low migration rate (m≥0.005/generation), where demographic independence and biologically meaningful differentiation still exist (e.g., [START_REF] Cullingham | Confidently identifying the correct K value using the DeltaK method: when does K = 2?[END_REF][START_REF] Kalinowski | The computer program STRUCTURE does not reliably identify the main genetic clusters within species: simulations and implications for human population structure[END_REF]. To better differentiate DIPs, we examined spatially explicit principal components with geographic information using sPCA [START_REF] Jombart | Revealing cryptic spatial patterns in genetic variability by a new multivariate method[END_REF]. Spatial distances were based on type 1 (Delaunay triangulation) connection network. We tested for significant evidence of structure in the sPCA using the Eigenvalue test "spca_randtest" [START_REF] Montano | An Eigenvalue test for spatial principal component analysis[END_REF] with 9999 permutations in the R package adegenet (v. 2.1.1; [START_REF] Jombart | adegenet: a R package for the multivariate analysis of genetic markers[END_REF]. Geographical subsets of the data were analyzed hierarchically to evaluate structure at decreasing spatial scales.

Given the nearly continuous distribution of harbor porpoises along a coastline, we tested for genetic isolation by distance using Mantel tests for correlation of both individual and population genetic distances with geographic distances using the adegenet R package.

We used Euclidean distance for individuals (Euclidean distance among vectors of allele frequencies) and pairwise FST distances between population strata, and straight-line geographic distances between individual samples or average latitude/longitude position of samples clustered into a priori geographic populations. We also used Monmonier's algorithm [START_REF] Monmonier | Maximum-difference barriers: an alternative numerical regionalization method[END_REF] as implemented in adegenet. Putative boundaries between populations were inferred based on the default threshold value (third quartile of all distances between neighbors).

For both nuclear and mtDNA data, we tested for a priori population divergence using pairwise estimation of FST with MCMC resampling implemented in the strataG R package. Multiple stratification schemes were tested based on a priori information gained from the Bayesian population genetic and sPCA analyses described above, previously defined management stocks [START_REF] Carretta | U.S. Pacific marine mammal stock assessments[END_REF], and gaps in harbor porpoise distribution [START_REF] Forney | Detecting trends in harbor porpoise abundance from aerial surveys using analysis of covariance[END_REF]. For mtDNA, we generated a median joining network (MJN) using the program POPART [START_REF] Leigh | POPART: fullfeature software for haplotype network construction[END_REF].

We used distance-based redundancy analysis (dbRDA) to investigate genotypeenvironment associations and identify microhaplotypes potentially under selection [START_REF] Forester | Comparing methods for detecting multilocus adaptation with multivariate genotype-environment associations[END_REF]. dbRDA identifies how groups of SNPs or microhaplotypes covary in response to the multivariate environment. It is well-suited to isolation-by-distance demographic scenarios, maintaining both high true positive and low false positive rates [START_REF] Forester | Comparing methods for detecting multilocus adaptation with multivariate genotype-environment associations[END_REF]. Environmental variables (Table 1) were extracted from longitudinal oceanographic data matched to collection date, latitude and longitude of samples. We used Bray-Curtis dissimilarity [START_REF] Bray | An Ordination of the Upland Forest Communities of Southern Wisconsin[END_REF] to calculate the microhaplotype distance matrix. This approach quantifies the dissimilarity among individuals based on their multilocus genotypes, and is equivalent to proportion of shared alleles [START_REF] Shirk | A comparison of individualbased genetic distance metrics for landscape genetics[END_REF]. We performed dbRDA separately for porpoises from the outer coastal and inland waterways regions, because data for all environmental predictors were not available for all locations (Table 1).

For dbRDA, we first removed microhaplotype loci with heterozygosity less than 0.05, then removed individuals with missing data for the retained environmental predictors. We produced three data sets for each region (i.e., outer coastal and inland waterways), representing three thresholds of missing genotype data across individuals: 25%, 20%, and 15%. We imputed missing values for each data set using snmf in the LEA package v.

3.1.2 [START_REF] Frichot | LEA: An R package for landscape and ecological association studies[END_REF], testing values of K from 1-5, and alpha (regularization parameter) values of 10, 100, and 1000. All runs used 25 repetitions, 200 iterations, and a 5% cross entropy withholding. We then performed a dbRDA for each imputed data set, retaining three constrained axes for outlier analysis. We identified candidate microhaplotypes under selection using robust (e.g., not sensitive to outliers) Mahalanobis distance [START_REF] Capblancq | Evaluation of redundancy analysis to identify signatures of local adaptation[END_REF], which identifies outlier microhaplotypes based on their constrained ordination loadings in multidimensional space. We accounted for confounding factors in the dbRDA, such as population structure and isolation-by-distance, using the genomic inflation factor [START_REF] Francois | Controlling false discoveries in genome scans for selection[END_REF], and applied a false discovery rate cutoff of 0.1 [START_REF] Storey | Statistical significance for genomewide studies[END_REF] to identify outlier microhaplotypes showing relationships with environmental variation. Finally, we compared detections across missing data thresholds.

RESULTS

We used a total of 431 mtDNA control region sequences (395bp) from previously published (N=363, including re-sequenced samples) and newly generated (N=68) sequences, and SNP data for 296 porpoises ranging from the southern extent of the harbor porpoise range in Southern California, USA, to British Columbia, Canada (Supplemental Table S1, Figure S1). Resequencing of samples with previously published haplotypes resulted in ten haplotype changes from the [START_REF] Chivers | Small-scale population structure of eastern North Pacific harbour porpoises (Phocoena phocoena) indicated by molecular genetic analyses[END_REF][START_REF] Hanson | Using location data from telemetry tagged marine mammals to improve stock assessments[END_REF] studies, and 53 from [START_REF] Crossman | Population structure and intergeneric hybridization in harbour porpoises Phocoena phocoena in British Columbia, Canada[END_REF], mostly due to resolution of ambiguous positions in the previous electropherograms that resulted in synonymization of multiple haplotypes from each of those studies. We identified 52 harbor porpoise haplotypes, 22 of which were found in only a single individual, and six Dall's porpoise haplotypes (supplemental Tables S4, S5). Our GT-seq locus panel consisted of 340 loci, of which 290 were genotyped in at least 55% of the samples and were polymorphic in harbor porpoises. Two additional loci were polymorphic only in Dall's porpoises and were used to identify hybrids between the 2 species. Of the 290 loci, 151 (52%) contained a single SNP, while the remaining 139 (48%) contained ≥2 SNPs, genotyped as microhaplotypes. None of the loci deviated significantly from HWE, and significant LD was only detected in one locus pair in one of the three tested geographic strata; no loci were removed based on these tests of HWE or LD.

DNA quantity and quality varied substantially among samples, resulting in variable number of completed genotypes, and an inverse correlation of error rate with the number of completed genotypes. Arbitrarily changing the cut-off value for percent completed genotypes can result in slight changes (see error rates below) to overall data quality, but at the cost of reduction in sample sizes in individual strata, reducing statistical power to detect structure. To maximize the sample sizes across strata, we used 20% (≥58 of 292) genotyped loci as the minimum cut-off, resulting in 296 genotyped samples (after removal of unintentional duplicates (see below): 280 harbor porpoise, 11 Dahl's porpoise, 5 hybrids). Of the harbor porpoise samples, 72% were genotyped at >90% of the loci, and 83% were genotyped at >80% of the loci (supplemental Table S1).

Average per-allele error rates for SNPs were calculated from 32 samples genotyped in duplicate from separate GT-seq amplicon libraries, based on single SNPs in 270 loci that were genotyped in >50% of samples. Intentionally replicated samples for which there was sufficient data in both replicates (N=32) had matching genotypes at an average of 96% of the loci (range 87-100%). For samples genotyped at >80% of the loci, the estimated error rate (based on 23 replicated sample pairs) was 0.010/allele. For samples with lower genotype completion rates (46 in the harbor porpoise data set), the mean error rate estimate increased to 0.045/allele (based on N=9 replicate pairs). Nine sample pairs (excluding Dall's porpoise samples, which were too homozygous in this data set) were identified as unintentional duplicates based on at least 85% identical genotypes (range 96-100%), and one from each pair was removed from further analysis. All of the genetic duplicates were inadvertent duplicate samples from the same individual, usually due to sample being stored in different collections with different identification codes. We detected three potential first order relatives (full siblings or parent-offspring pairs) based on the "wang" estimator, and six (including two of the three from the "wang" estimator) based on the "Mxy" estimator (supplemental Table S6) (see [START_REF] Kraemer | Demerelate: calculating interindividual relatedness for kinship analysis based on codominant diploid genetic markers using R[END_REF] for details). All samples were retained for some of the subsequent analyses, but one individual from each putative pair was removed to control for the effects of sampling closely related individuals in STRUCTRURE, sPCA, and FST analyses.

Putative hybrids between parapatric Dall's and harbor porpoise species have been previously identified based on phenotype and genetic profiles [START_REF] Crossman | Population structure and intergeneric hybridization in harbour porpoises Phocoena phocoena in British Columbia, Canada[END_REF].

We genotyped nine known Dall's porpoises and had phenotypic or previous genotypic indication (based on microsatellites; [START_REF] Crossman | Population structure and intergeneric hybridization in harbour porpoises Phocoena phocoena in British Columbia, Canada[END_REF] of 11 putative hybrids in our final data set. Assignment analysis of all samples using STRUCTURE with the number of clusters set to k = 2 correctly assigned all nine Dall's porpoises to one group, and additionally assigned two of the putative hybrids to the same group with 100% probability. Five of the putative hybrids were assigned to the harbor porpoise group with ≥99.9% probability, and four samples were identified as F1 hybrids with assignment probabilities to each species group between 45% and 55%. One additional sample previously identified as a harbor porpoise was also identified as an F1 hybrid (49%/51% assignment to the 2 species groups). The admixture plot from ten combined structure analyses is presented in supplemental Figure S2a. All F1 intergeneric hybrids were from samples collected in the San Juan Islands and Oregon-Washington coast geographic strata, between latitudes 47º and 49º N. All remaining putative harbor porpoise samples were assigned to the second group with >99.6% probability. As the targeted SNP loci were ascertained only from harbor porpoise samples, Dall's porpoise samples had unsurprisingly low diversity (Table 2), but the use of additional SNPs in microhaplotypes provided variable loci useful for species assignment and hybrid identification. Samples identified as hybrids were removed from subsequent harbor porpoise analyses. STRUCTURE is useful for identifying population differentiation at the ESU level, where divergence is sufficient to allow high probability of assignment of samples to populations or clusters [START_REF] Waples | What is a population? An empirical evaluation of some genetic methods for identifying the number of gene pools and their degree of connectivity[END_REF]. We started with STRUCTURE in a hierarchical analysis of the SNP data to identify evolutionarily divergent populations within our harbor porpoise sample distribution. Structure analysis with all samples assigned to geographic strata a priori (supplemental Table S1) did not provide strong evidence of multiple divergent groups along the U.S. west coast, with the number of groups most likely ≤2 (based on ΔK and CLUMPAK similarity score >0.993 for 10/10 replicates for K1-2; supplemental Figure S2b). However, assignment probabilities among three groups differed substantially across geographic region, especially in the three southernmost geographic strata (Figure 2). Subsets of samples representing only high-quality samples (<80% complete) and smaller geographic regions did not result in any additional evidence of population structure (supplemental Figures S2c-e).

Spatial principle component analysis (sPCA; [START_REF] Jombart | Revealing cryptic spatial patterns in genetic variability by a new multivariate method[END_REF] has been shown to be a useful tool for detecting patterns of genetic variability in harbor porpoise [START_REF] Fontaine | Mixing of porpoise ecotypes in southwestern UK waters revealed by genetic profiling[END_REF][START_REF] Lah | Spatially explicit analysis of genome-wide snps detects subtle population structure in a mobile marine mammal, the harbor porpoise[END_REF], and was used to explicitly combine geographic and nuclear genetic data and to investigate spatial patterns of genetic variation without assumptions of Hardy-Weinberg expectations or linkage disequilibrium [START_REF] Montano | An Eigenvalue test for spatial principal component analysis[END_REF].

Plots of spatial principle components across the study range indicated strong evidence of structure (spca-randtest p < 0.001) and show clear separation between samples representing the two southern-most strata, Morro Bay and Monterey Bay, along the first axis (PC1, Figure 3A, 3B), with less separation but a north-south gradient along spatial PC axes two and three (Figure 3A; supplemental Figure S3A). Additional geographic population structure is revealed by analyzing geographic subsets of the data. Hierarchical analysis of sub-regions reveals population structure within the inland waterways (p = 0.004; Figure 3G), and within the Washington inland waters managment stock (excluding the Strait of Georgia strata samples) (p = 0.019; supplemental Figure S3H). sPCA of samples from coastal waters (Figure 3B), excluding inland waterways in Washington and British Columbia, shows strong evidence of structure (p < 0.001), with the first two spatial PC's highlighting unique clusters corresponding to Morro Bay (PC1), and to a lesser extent the remaining strata along PC2. Hierarchical analysis of subsets of neighboring strata indicated significant structure in analyses of all neighboring strata pairs and within the northern California/southern Oregon management stock (p = 0.017; Figure 3D; supplemental Figure S3D). Plots of the first four individual spatial PCs for these stratification sets are in supplemental Figure S3. Replicate analyses based on the smaller number of high-quality samples (>80% complete genotypes) showed similar patterns, but resulted in non-significant p-values in some comparisons (Inland waters, southern OR/OR-WA, and northern California/southern Oregon; Supplemental Figure S3).

A Mantel test supported (p < 0.001) isolation by distance (IBD) for nuclear SNP data along the Pacific coast (excluding inland waterways) for individual distances. A scatterplot of genetic and geographic distances (Supplemental Figure S4) shows discontinuities suggestive of differentiated populations rather than continuous clines of genetic differentiation. The Monmonier's algorithm was used to infer locations of genetic boundaries or discontinuities, and identified the major boundary between coastal Washington and western Vancouver Island at the default threshold level (Figure 4A), with lesser boundaries between Morro Bay and Monterey Bay, and Monterey and SF-RR on the outer coast. Within the inland water ways, the major discontinuities were in the southern Strait of Georgia (Figure 4B). It is not straightforward to interpret the Monmonier plots because the samples are not distributed broadly in two dimensions, but the threshold arrows indicate regions where genetic distances appear to be greater than expected based on geographic distance (see [START_REF] Blair | A simulation-based evaluation of methods for inferring linear barriers to gene flow[END_REF] for discussion of methods of detecting barriers to gene flow).

Genetically similar groups (Figure 3) and IBD discontinuities (Figure 4) coincided closely with several previously defined geographical strata (Figure 5; [START_REF] Chivers | Small-scale population structure of eastern North Pacific harbour porpoises (Phocoena phocoena) indicated by molecular genetic analyses[END_REF][START_REF] Chivers | Additional genetic evidence for population structure of Phocoena phocoena off the coasts of California[END_REF]. Genetic differentiation between adjacent a priori strata was tested using pairwise divergence (FST), for both nuclear data and mtDNA (Table 3).

Results are presented based on inclusion of all samples, but pairwise divergence among adjacent strata after removal of one from each pair of potential first order relatives, and after removal of samples with <50% or <80% completed genotypes, did not qualitatively alter results (supplemental Table S7). As expected when larger numbers of loci are used to calculate FST [START_REF] Willing | Estimates of genetic differentiation measured by F(ST) do not necessarily require large sample sizes when using many SNP markers[END_REF], the FST point estimates did not change substantially (<0.004), but exclusion of samples reduced the sample sizes in some geographically-defined strata, resulting in loss of statistical power to detect structure among some non-adjacent strata (supplemental Table S7D-F). Along the outer coast the only genetic distances between neighboring strata pairs that were significant (p < 0.05) for both mtDNA and nuclear DNA were between Morro Bay at the southern end of the distribution and Monterey Bay (FST = 0.020 (nDNA); 0.310 (mtDNA)). Nuclear locus frequencies were also significantly different between three of the more southern strata pairs (Monterey/SF-RR (FST = 0.005), N. CA/S. OR (FST = 0.008), S. OR/ORWA coast (FST = 0.006), while mtDNA haplotype frequencies were significantly different between only the two northern-most coastal strata pairs (ORWA coast /Spike Rock (FST = 0.049), and Spike Rock/W. Vancouver Is. (FST = 0.128)). Frequency plots of the common mtDNA haplotypes (excluding haplotypes that occurred in fewer than 5 samples) show clear differences among a priori geographical strata (Figure 5, supplemental Figure S5), with near fixation of haplotypes in Morro Bay to the south, and private haplotypes in BC to the north. Relative density data showed low-density regions along the Big Sur coastline (35.5-36.5º N), just south of 38º N, and just north of 39º N at the existing management stock boundary. These low relative density areas correspond to existing boundaries between the four southern-most management stocks. The reason for low density for one region, the Big Sur coastline, is likely due to the very narrow shelf, but reasons for density variation in other regions remain unclear.

Inland waters strata were not differentiated by nuclear marker frequencies except for the Strait of Georgia, which was significantly differentiated from all inland waters and outer U.S. coast strata, but not W. Vancouver Is. or BC (Table 3). The inland waters strata were differentiated from the nearest coastal strata in mtDNA (FST between 0.023 and 0.247, with most point estimates being significant at p < 0.05). Within the inland waters, divergence in mtDNA was significant for Neah Bay, in the outer Strait of Juan de Fuca, versus San Juan Islands (FST = 0.066), and from both neighboring coastal strata (FST = 0.141 (Spike Rock); 0.247 (W. Vancouver Island)). Puget Sound, which has been recolonized in the last 18 years [START_REF] Evenson | Disappearance and return of harbor porpoise to Puget Sound: 20 year pattern revealed from winter aerial surveys[END_REF], was not significantly differentiated from either the neighboring San Juan Islands, or Neah Bay in the Strait of Juan de Fuca.

The genetic discontinuity identified immediately north of the U.S.-Canadian border at approximately 49.1°N. latitude (Figure 4) warranted an adjustment of the boundary between the Strait of Georgia and San Juan Islands strata from 48.8° to 49.1°, resulting in reassignment of ten and five samples (mtDNA, nDNA analyses, respectively) from the Strait of Georgia to the San Juan Islands strata (Figure 4B). Genetic divergence remained significant following the reassignment of these samples and FST increased for both nDNA (0.005 to 0.006) and mtDNA (0.095 to 0.153).

For the genotype-environment association analyses (dbRDA), we analyzed three data sets for each region (i.e., outer coastal and inland waterways), representing three thresholds of missing genotype data across individuals: 25%, 20%, and 15%. We first imputed the data sets with snmf using the optimized settings (identified by minimizing the average cross entropy, with identical optimized parameters for all data sets): alpha=10, and K=1 (supplemental Table S8). We modified the default genomic inflation factors to produce pvalue distributions that better met the uniform distribution assumption (per Francois et al., 2016, p-value histograms provided in supplemental Figures S6-S11), providing a balance between true positive detections and false positives driven by population structure and other potential confounders. We did not include other corrections for confounding factors in the GEA since population structure is low across the study area (e.g., Table 3), and dbRDA is robust to isolation-by-distance demographic scenarios [START_REF] Forester | Comparing methods for detecting multilocus adaptation with multivariate genotype-environment associations[END_REF] such as this one. Environmental predictor variables within each region were not highly correlated (i.e., all pairwise correlations were less than |0.8| and all variance inflation factors, which measure multicollinearity, were less than five, where values >10 can be problematic), and none of the retained predictors were highly correlated with either latitude or longitude (supplemental Tables S9,S10). For the outer strata and inland waterways, we identified as candidates those microhaplotypes that were detected in at least two of the three thresholded data sets. This produced 22 candidates for the outer coastal region and six candidates for the inland waterways region. Candidate microhaplotypes for the outer coastal region showed the strongest relationships with the mean and standard deviation of daily sea surface temperature (sst.mn and sst.SD), and the mean of daily sea surface height (ssh.mn). Triplots of the dbRDA results illustrated relationships among individual's multilocus genotypes, the candidate microhaplotypes, and environmental predictors. For example, the southernmost Morro Bay individuals, located in the ordination space as a function of their genotypes at the candidate loci, showed strong relationships with increasing ssh.mn and sst.mn and decreasing ssh.SD, indicating potential adaptation to local environmental conditions (Figure 6A). The smaller inland waterways data set showed some differentiation between the Neah Bay and San Juan Islands individuals based on sea surface temperature, with the Puget Sound individuals showing no relationships, possibly due to their most recent re-colonization (Figure 6B). Because loci under selection could affect population structure, we confirmed that identification of divergence between adjacent a priori strata were consistent with and without the loci identified by dbRDA analysis (supplemental Table S7). Full dbRDA triplots and microhaplotype biplots are provided in supplemental Figures S12-S15 for the outer coastal data and supplemental Figures S16-S19 for the inner waterways.

Discussion

The North Pacific harbor porpoise is one of several geographically and genetically described subspecies of harbor porpoise, some of which consist of multiple ecotypes [START_REF] Ben Chehida | Evolutionary history of the porpoise family (Phocoenidae) across the speciation continuum: a mitogenome phylogeographic perspective[END_REF][START_REF] Fontaine | Harbour porpoises, Phocoena phocoena, in the Mediterranean Sea and adjacent regions: Biogeographic relicts of the last glacial period[END_REF][START_REF] Fontaine | Postglacial climate changes and rise of three ecotypes of harbour porpoises, Phocoena phocoena, in western Palearctic waters[END_REF][START_REF] Fontaine | History of expansion and anthropogenic collapse in a top marine predator of the Black Sea estimated from genetic data[END_REF].

Studies of the North Atlantic subspecies have suggested both historical biogeographic processes [START_REF] Fontaine | Rise of oceanographic barriers in continuous populations of a cetacean: the genetic structure of harbour porpoises in Old World waters[END_REF][START_REF] Fontaine | Postglacial climate changes and rise of three ecotypes of harbour porpoises, Phocoena phocoena, in western Palearctic waters[END_REF] and ecological processes [START_REF] Fontaine | Mixing of porpoise ecotypes in southwestern UK waters revealed by genetic profiling[END_REF]Fontaine, Tolley, et al., 2007;[START_REF] Lah | Spatially explicit analysis of genome-wide snps detects subtle population structure in a mobile marine mammal, the harbor porpoise[END_REF] resulting in divergent populations and limited dispersal and introgression between adjacent types. In the North Pacific, evidence of structure based on both biogeographic [START_REF] Taguchi | Mitochondrial DNA phylogeography of the harbour porpoise Phocoena phocoena in the North Pacific[END_REF] and ecological/spatial divergence [START_REF] Chivers | Small-scale population structure of eastern North Pacific harbour porpoises (Phocoena phocoena) indicated by molecular genetic analyses[END_REF][START_REF] Chivers | Additional genetic evidence for population structure of Phocoena phocoena off the coasts of California[END_REF][START_REF] Crossman | Population structure and intergeneric hybridization in harbour porpoises Phocoena phocoena in British Columbia, Canada[END_REF] has been more limited, but suggests that similar processes may have acted across the ranges of subspecies in both ocean basins.

Previous genetic studies of North Pacific harbor porpoise have been limited by sample availability, sample quality, and low power to detect population genetic differences within this nearly continuously distributed coastal species. We have developed a set of SNP and microhaplotype loci that provide sufficient genetic power to detect low levels of population structure, while allowing us to make use of poor-quality tissue samples available from beach-cast and fishery bycatch carcasses. We used hierarchical partitioning to explore evidence of population structure across the range, and to infer potential barriers to geneflow. Correlation of remotely-sensed environmental variables with genetic patterns identified factors potentially influencing local adaptation. Our results indicate that North Pacific harbor porpoises exhibit genetic discontinuity and limited dispersal that may be associated with habitat variability and local adaptation.

Previous evidence indicated intrageneric hybridization between the more pelagic Dall's porpoise and the coastal harbor porpoise, based on morphology and genetic analysis of eight microsatellite loci [START_REF] Crossman | Population structure and intergeneric hybridization in harbour porpoises Phocoena phocoena in British Columbia, Canada[END_REF]. While the previous genetic analyses suggested bi-directional hybridization and introgression, the more powerful dataset used here indicates that many of the putative hybrids (based on genetics) could be assigned with >98% probability to one or the other species (N=5 P. phocoena; N=2 P. dalli), or as F1 hybrids (N=5) between male harbor porpoise and female Dall's porpoise. Sightings of putative hybrids in the wild have also all been associated with Dall's porpoises [START_REF] Baird | An intergeneric hybrid in the family Phocoenidae[END_REF]. Harbor porpoise males have disproportionately large testes and have been observed to exhibit forceful and fast ambush mating [START_REF] Keener | The sex life of harbor porpoises (Phocoena phocoena): Lateralized and aerial behavior[END_REF], suggesting a greater reliance on sperm competition rather than mate choice in male harbor porpoise. This mating strategy might result in occasional intergeneric hybrids with Dall's porpoise in areas where their ranges overlap [START_REF] Baird | An intergeneric hybrid in the family Phocoenidae[END_REF]. Although the lack of individuals exhibiting intermediate probabilities of assignment to both species (between 50% and 100%) suggests hybrids are usually infertile, there have been sightings of putative hybrid females with neonatal calves [START_REF] Willis | Natural hybridization between Dall's porpoises (Phocoenoides dalli) and harbour porpoises (Phocoena phocoena)[END_REF].

Within harbor porpoises, our results are consistent with isolation by distance along a mostly linear range from southern California to western Vancouver Island, and through the inland waterways, but also suggest some regions of higher genetic divergence than expected based on geographic distance alone. The two largest population splits are between Morro and Monterey Bays in the south, and between the Strait of Georgia and San Juan Islands in the north, where multiple methods suggest divergence in both mitochondrial DNA and nuclear markers. This break between the Strait of Georgia and San Juan Islands, as well as between the San Juan Islands and the outer Strait of Juan de Fuca, is consistent with the limited movements of tagged harbor porpoises in this area [START_REF] Hanson | Using location data from telemetry tagged marine mammals to improve stock assessments[END_REF]. Another population break is suggested between the Monterey Bay and San Francisco sample sets, based on the Monmonier analysis of the connection network (Figure 4), where FST divergence is significant for nuclear loci but not for mtDNA.

Although the Morro Bay population is designated as a separate management stock based on a hiatus in distribution at the Big Sur coast [START_REF] Carretta | U.S. Pacific marine mammal stock assessments[END_REF], our results provide the first evidence for genetic differences. The most common mtDNA haplotype (CR01) was found in every a priori geographic group except Morro Bay, and the most common haplotype in Morro Bay (CR02) was found in 86% of the Morro Bay samples, but less than 33% of any other population (range 0 -33%). Two (CR30, CR42) of the three other haplotypes found in Morro Bay were unique to that population, and differed from CR02 by only one nucleotide change. The third haplotype (CR03) was distributed across most of the range, but represents a common haplotype that is more similar to CR02 (3bp different) than to the cluster including CR01 (7bp different) and the majority of the other haplotypes distributed throughout the rest of the range (see haplotype network, supplemental figure S5). This suggests long-term isolation of the southernmost population, as well as a persistently small population size, or severe or recurrent bottlenecks. Annual estimates of tens to hundreds of harbor porpoises were killed in gillnet entanglements between approximately 1960 and 2001, resulting in a reduction of the population by 30-97% [START_REF] Barlow | An assessment of the status of the harbour porpoise in central California[END_REF]. The population has since begun recovering, from a low of 560 in 1990 to 4255 in 2012 [START_REF] Carretta | U.S. Pacific marine mammal stock assessments[END_REF]Forney et al., in press). Our results suggest that recovery has been due to internal recruitment rather than migration from the larger Monterey Bay population to the north.

The northernmost strata were the most differentiated, with high mtDNA divergence metrics and the highest frequencies of the most common control region haplotype (CR01; 66% W. Vancouver Is., 63% Strait of Georgia). The frequency of this haplotype ranged from zero to 50% in other strata. The BC stratum (excluding Strait of Georgia samples) spanned a large geographic range with few samples (N=5), and included two private mtDNA haplotypes, both found in two samples. These suggest different haplotype composition in this region, but we cannot rule out sampling of related individuals or isolated regional populations, as the sample size is small and shared haplotypes were found in samples that were collected close in space and time. Combined, these data suggest that there is a haplotype gradient between the northern and southern portions of the range, possibly representing interglacial range expansion [START_REF] Taguchi | Mitochondrial DNA phylogeography of the harbour porpoise Phocoena phocoena in the North Pacific[END_REF] or the presence of a historical phylogeographic separation that has more recently re-connected, as has been suggested for populations in the eastern Atlantic Ocean [START_REF] Fontaine | Rise of oceanographic barriers in continuous populations of a cetacean: the genetic structure of harbour porpoises in Old World waters[END_REF][START_REF] Fontaine | Postglacial climate changes and rise of three ecotypes of harbour porpoises, Phocoena phocoena, in western Palearctic waters[END_REF].

Harbor porpoise habitat in the eastern North Pacific includes the California Current System (CCS), recognized as one of the most productive marine ecosystems on the planet, but spans topographically, oceanographically and temporally complex regions from temperate to Arctic waters. Aside from preference for coastal waters less than 200 m deep, little is known about factors influencing suitable habitat for Pacific harbor porpoises, which vary in population density and abundance along the west coast and inland waters of North America. The remotely sensed and modeled oceanographic variables available throughout most of this range (upwelling index, mixed-layer depth, sea surface temperature and sea surface height) are all proxies for habitat variability with complex relationships to wind, currents, topology, isocline depth, productivity, freshwater input and seasonal and interannual effects [START_REF] Castelao | Upwelling jet separation in the California Current System[END_REF][START_REF] Hickey | Alongcoast structure and interannual variability of seasonal midshelf water properties and velocity in the Northern California Current System[END_REF][START_REF] Venegas | Satellite-derived variability in chlorophyll, wind stress, sea surface height, and temperature in the northern California Current System[END_REF]. Our genotype-environment association analyses provide initial evidence for local adaptation to environmental variability across the study range.

Given the limited scale of genomic sampling in this study these findings should be seen as preliminary, providing a basis for future investigation of local adaptation across the complex seascape environment inhabited by harbor porpoise (see, e.g., [START_REF] Fontaine | Mixing of porpoise ecotypes in southwestern UK waters revealed by genetic profiling[END_REF]Fontaine, Tolley, et al., 2007;[START_REF] Nielsen | Oceanic movements, site fidelity and deep diving in harbour porpoises from Greenland show limited similarities to animals from the North Sea[END_REF].

An overall correlation with latitude, which might be expected across a large latitudinal range, appears to strengthen in the regions south of Oregon. In central Oregon, there is an oceanographic shift around Cape Blanco, where the coastal upwelling jet separates from the coast (~15-30 km offshore), becoming an oceanic jet (>100 km offshore) [START_REF] Castelao | Upwelling jet separation in the California Current System[END_REF]. Populations in the northern half of the coastal range are associated with several environmental variables (e.g., mean and SD of the mixed layer depth, and high variation in both sst and ssh), while those in the central and southern portions are associated with others (e.g., high mean and low variation in ssh) (Figure 6). There is also evidence of population structure across that region, with significant nuclear genetic differentiation among populations from southern Washington to northern California (Table 3). There was a particularly strong correlation between increased mean ssh and decreased ssh variability in the southernmost Morro Bay population, which is also the most genetically distinct, suggesting possible local adaptation linked to oceanographic processes (e.g., thermocline and/or upwelling). While these links are indirect, they provide the first evidence of local adaptation in addition to demographic independence among harbor porpoise populations along the west coast of North America. The Morro Bay stock was previously recognized based only on a distribution hiatus and evidence of historical population decline due to extensive fishery bycatch [START_REF] Barlow | An assessment of the status of the harbour porpoise in central California[END_REF], as well as blubber pollutant ratios that differed from other areas [START_REF] Calambokidis | Chlorinated hydrocarbon concentrations and their use for describing population discreteness in harbor porpoises from Washington, Oregon, and California[END_REF]. It is potentially subject to offshore energy production disturbance [START_REF] Forney | Nowhere to go: noise impact assessments for marine mammal populations with high site fidelity[END_REF], and is also at the southern edge of the species' geographic range, and most likely to be impacted by climate change [START_REF] Learmonth | Potential effects of climate change on marine mammals[END_REF][START_REF] Ruiz-Cooley | Temporal variation in pelagic food chain length in response to environmental change[END_REF]. Our results provide strong support for continued management of this population as a separate management stock, especially in light of potential impacts of coastal development and climate change.

Harbor porpoise populations along the west coast of North America have historically experienced substantial fisheries bycatch in portions of the range [START_REF] Barlow | An assessment of the status of the harbour porpoise in central California[END_REF], and multiple lines of evidence indicate that there is limited dispersal among regions. Within the U.S., this led to identification of five coastal management stocks between southern California and the northern U.S. border (coinciding with four of the geographic strata in this study, plus the combined S. Oregon and N. California strata), and one stock in the inland waterways in Washington [START_REF] Carretta | U.S. Pacific marine mammal stock assessments[END_REF], with stock boundaries placed to coincide with areas of lower density. The genetic results presented here provide additional evidence of limited movement among the currently defined management stocks, as well as between some geographically defined strata within existing stocks, consistent with previous studies of genetics and contaminant levels [START_REF] Calambokidis | Chlorinated hydrocarbon concentrations and their use for describing population discreteness in harbor porpoises from Washington, Oregon, and California[END_REF][START_REF] Chivers | Small-scale population structure of eastern North Pacific harbour porpoises (Phocoena phocoena) indicated by molecular genetic analyses[END_REF], and limited movements of tagged porpoises in the inland waterways of Washington [START_REF] Hanson | Using location data from telemetry tagged marine mammals to improve stock assessments[END_REF]. In particular, sPCA and FST analyses of nuclear data show separation of the N. California and S. Oregon strata (sPCA randtest p = 0.017) with a small but significant divergence (FST=0.008) in nuclear allele frequencies (Table 3, Figure 3D, supplemental Figure S3D). We suggest that the lack of statistically significant difference in mtDNA is likely a result of low statistical power due to small sample size in N. California (n=11) rather than results that contradict the finding of differences between these strata based on nuclear markers.

Placement of a potential stock boundary as indicated in Figure 5 is suggested by the drop in density in a region where there is a gap in our sample distribution, though the location of the density minimum in Figure 5 varies through time (supplemental Figure S20), either due to sampling differences among surveys or small north-south shifts of populations across years.

Genomic methods are evolving rapidly, increasing both the number and variety of genetic markers that can be used to understand evolution, population structure, historical demography, and ecological adaptations (e.g., [START_REF] Tan | Applications of Next-Generation Sequencing Technologies and Computational Tools in Molecular Evolution and Aquatic Animals Conservation Studies: A Short Review[END_REF]. This study of harbor porpoise population structure has built on previous research by expanding the number and geographic distribution of samples and the number and type of genetic markers. We also apply new analytical methods to infer patterns of spatial genetic variation, and investigate their correlation with environmental variables that may drive local adaptation.

Genetic analyses have been identified as potentially high value for stock delineation [START_REF] Martien | The dip delineation handbook: A guide to using multiple lines of evidence to delineate demographically independent populations of marine mammals[END_REF] and these results will be useful for harbor porpoise management stock structure refinement. Nevertheless, inference remains limited by the uneven distribution of opportunistic samples, potential shifts in populations during multiple decades of sample accumulation, and limited availability of environmental variables across the spatial and temporal scale of this study. Additional studies using genome-wide genetic data from across the range are needed to more fully understand habitat use and local adaptation in the eastern North Pacific, and their importance for harbor porpoise management and conservation.

Figure 5 .Figure 6 .

 56 Figure 5. Pie charts represent frequencies of haplotypes represented by ≥ 5 samples in the stratum, plus two haplotypes (CR05, CR35) that were found only in the BC strata, in two samples each (out of 5 total). Relative porpoise density estimates (porpoise per kilometer surveyed) are plotted by latitude, averaged across multiple survey periods (blue line). The current U.S. stocks are named and boundaries are represented by horizontal black lines (as in Figure 1). The proposed boundary location for splitting the Northern California/Southern Oregon stock into two stocks is shown with a dashed red line.

Table 1 .

 1 Environmental predictors for distance-based redundancy analysis. SD = standard deviation.

	Analysis	Environmental predictor	Abbreviation Predictor resolution	Source
	Outer	Mean daily sea surface	sst.mn	Averaged over 0.1 degrees	Regional Ocean Modeling System
	coastal	temperate		around sample location	ROMS/CCSRA 31-year reanalysis
	populations	SD of daily sea surface	sst.SD	Calculated within 0.3 degrees	(1991-2010) or Near Real Time
		temperate		around sample location	(2011+), from the U.C. Santa Cruz
		Mean daily sea surface height	ssh.mn	Averaged over 0.1 degrees	Ocean Modeling and Data Assimilation
				around sample location	group;
		SD of daily sea surface height	ssh.SD	Calculated within 0.3 degrees	http://oceanmodeling.ucsc.edu/; (Moore
				around sample location	et al., 2011)
		Mean daily mixed layer depth	ild.mn	Averaged over 0.1 degrees
		(depth at which the temperature		around sample location
		is -0.5°C from surface
		temperature)
		SD of daily mixed layer depth	ild.SD	Calculated within 0.3 degrees
		(depth at which the temperature		around sample location
		is -0.5°C from surface
		temperature)
		Coastal Upwelling Transport	CUTI	Calculated every 1 degree	http://mjacox.com/upwelling-indices/
		Index		latitude along West Coast	(Jacox et al., 2018)
				from 31N-47N
		Biologically Effective	BEUTI	Calculated every 1 degree
		Upwelling Transport Index		latitude along West Coast
				from 31N-47N
	Inland	Mean daily sea surface	mSSTday.mn0	Averaged over 0.1 degrees	Multispectral Ultra-high Resolution Sea
	waterways	temperate (comparable to sst.mn	.1	around sample location	Surface Temperature (murSST);
	populations	above)			downloaded from NOAA/ERDDAP
		SD of daily sea surface	mSSTday.SD	Calculated within 0.3 degrees	(Simons, R. A. 2015. Environmental
		temperature (comparable to	0.3	around sample location	Research Division Data Access
		sst.SD above)			Program ERDDAP. https://upwell.
		Mean monthly sea surface	mSSTmon.mn	Averaged over 0.1 degrees	pfeg.noaa.gov/erddap.
		temperature	0.1	around sample location	Data set names: jplMURSST41SST
		SD of monthly sea surface	mSSTmon.SD	Calculated within 0.3 degrees	(daily) and jplMURSST41mday
		temperature	0.3	around sample location	(monthly)

Table 2 .

 2 SNP and microhaplotype genotype summary information, based on 292 loci.

	Species	Samples	Hexp.	Hobs.	Monomorphic
	P. phocoena	281	0.366	0.356	2
	P. dalli	11	0.144	0.145	177
	hybrid	5	0.479	0.324	64

Table 3 .

 3 Pairwise FST for nuclear loci (lower left) and mtDNA haplotypes (upper right). Statistically significant values (p<0.05) are in bold. Sample sizes each population is shown in parentheses for nuclear loci on the left, and mtDNA on the top. The p-values for adjacent strata pairwise comparisons are provided in supplemental TableS9.Table4. dbRDA sample information for three data sets filtered for different levels of missing genotypes. Rows in italics report the number of retained individuals in each set of populations. GIF = genomic inflation factor; FDR = false discovery rate.

				Strait of		Puget			ORWA
			W. Vanc.	Georgia	San Juan	Sound	Neah Bay	Spike	Coast	SoOR	NoCal	SF_RR	Monterey	Morro
		BC (5)	Is (29)	(55)	Is (92)	(24)	(20)	Rock (18)	(30)	(29)	(17)	(30)	Bay (37)	Bay (28)
	BC (3)	--	0.223	0.249	0.097	0.147	0.251	0.127	0.151	0.140	0.126	0.108	0.125	0.604
	W. Vancouver Is (15)	-0.003	--	-0.005	0.088	0.076	0.247	0.128	0.002	0.005	0.004	0.046	0.130	0.542
	Strait of Georgia (20)	-0.006	0.005	--	0.095	0.066	0.232	0.152	0.011	0.017	0.012	0.064	0.159	0.526
	San Juan Is (27)	0.005	0.000	0.005	--	0.008	0.066	0.023	0.043	0.037	0.021	0.031	0.059	0.307
	Puget Sound (20)	0.003	0.003	0.009	0.001	--	0.039	0.057	0.029	0.020	0.009	0.046	0.093	0.454
	Neah Bay (21)	0.015	0.000	0.010	0.000	0.002	--	0.141	0.176	0.157	0.160	0.174	0.192	0.531
	Spike Rock (19)	0.002	-0.004	0.007	0.002	0.002	-0.001	--	0.049	0.083	0.033	0.005	0.024	0.233
	ORWA Coast (28)	0.015	-0.001	0.014	0.004	0.005	0.004	-0.001	--	-0.006	-0.023	-0.005	0.067	0.440
	SoOR (24)	0.009	0.002	0.011	0.002	-0.001	0.004	0.001	0.006	--	-0.027	0.019	0.090	0.480
	NoCal (11)	0.014	0.003	0.015	0.001	0.005	0.003	0.004	0.002	0.008	--	-0.012	0.059	0.448
	SF_RR(26)	0.012	0.004	0.015	0.005	0.004	0.005	-0.002	0.004	0.002	0.004	--	0.023	0.347
	Monterey Bay (29)	0.021	0.006	0.016	0.007	0.011	0.008	0.001	0.006	0.012	0.007	0.005	--	0.310
	Morro Bay (28)	0.027	0.028	0.036	0.022	0.024	0.028	0.023	0.026	0.025	0.024	0.020	0.020	--

Acknowledgements

We are grateful Andy Foote and Marie Louise for help with SNP discovery methods; Diana Baetscher, Cassie Columbus, Carlos Garza and the SWFSC Santa Cruz genetics lab for help with GT-seq optimizations, and Nathan Campbell and Ninh Vu for help with GT-seq method implementation and optimization; Eric Anderson, Thomas Ng, Anthony Clemento and Neil Thompson for help with MicrohaPlot software and methods; Philipp Kraemer for help with the Demerelate R package; Steven Head and John Shimashita, TSRI Sequencing core for sequencing assistance; David Anderson, Cascadia Research Collective for help with Washington Inland Waterways data and observations; Samuel Woodman for help with analysis of porpoise relative density data. Samples were generously provided by Lisa Ballance, Barbara Taylor, Jay Barlow at the SWFSC; Jan Straley, University of Alaska Southeast; Marilyn Dahlheim, NOAA National Marine Mammal Lab; Jorge Velez-Juarbe, Los Angeles County Museum of Natural History; Merrill Gosho, Bill Walker, Jeff Foster, Stephen Clausen, The Whale Museum, Olympic Coast National Marine Sanctuary and many contributors through the Marine Mammal Stranding Networks. Samples were collected under various Marine Mammal Protection Act permits, and stored in the SWFSC Marine Mammal and Sea Turtle Research (MMASTR) collection. Funding support for this study was provided by NOAA Office of Protected Resources. We are grateful to Aimee Lang and 3 anonymous reviewers for helpful comments on previous versions of this manuscript.

Data accessibility -Mitochondrial DNA haplotype sequences used in this study are complete or partial sequences of sequences previously available in Genbank (Accession numbers provided in Supplemental Table S2). Nuclear genotypes and environmental data associated with the samples: Dryad doi: https://doi.org/10.5061/dryad.4tmpg4f6v

Figure 4. A) Gabriel connection network from a Monmonier analysis of samples (adjusted to prevent sample location overlap) with local nuclear genetic differences above the threshold value indicated by thicker blue lines and arrows between geographic distance edges for A) the outer coastal samples, with subsets of samples shown for localized thresholds, and B) inland waterways of Washington and British Columbia. Inset maps are for orientation of sample networks. Thickness of the lines and arrows is proportional to the level of genetic difference above the threshold level. Colors correspond to a priori strata as in Figure 1.

San Juan Is.

Neah Bay