extfieldsofstudy:
Materials Science
Biology
Computer Science

text:
Multichannel Silicon Probes for Laminar Cortical and Hippocampal Recordings in Large Animals

1) Laminar information across different deep brain structures and cortical regions need to be decoded in order to understand the spatiotemporal ensembles that represent cognition and memory. Pig models are widely used in translational research, and may provide an excellent bridge between rodents and non-human primates for CNS disease models due to their gyrencephalic neuroanatomy and significant white matter composition. One of the major obstacles to addressing the research questions and approaches currently accessible to smaller-brained animals, are technical limitations in microelectrodes. 2) We tested various designs of multichannel silicon probes developed for large animal electrophysiology acutely in Yucatan pigs by recording neurophysiological signals from laminar structures under anesthesia. 3) Electrophysiological parameters of single units and local field potentials were analyzed to evaluate performance of given silicon probes acutely. 4) Cross-sectional area of multichannel silicon probes was found to be a crucial determinant of silicon probes’ performance. EDGE style probes have the highest yields during cortical and intra-hippocampal recordings in pigs, with a potential for use in chronic implantations and awake behavior. The new CAMB 64-channel probes with EDGE style and a poly-2 site arrangement have an even better single unit separation and a denser sampling of the laminar structure, identifying them as potential candidates for chronic implantations. This study provides a comprehensive analysis of acute silicon probes designed for large animal laminar electrophysiology.

Introduction
Oscillatory activity and corresponding power bands in the brain are different within layers of laminar structures as different inputs are driven at different frequencies depending on their area of origin [1][2][3][4][5]. Laminar analyses of oscillatory activity have proven valuable for analyzing circuit dynamics as well as changes between states in both cortex and hippocampus [5,6]. Traditionally, single neurons have been analyzed independently on the basis of their tuning to sensory stimuli or movement. Although tuning curve approaches are unaffected by growing numbers of simultaneously recorded neurons, newly developed techniques that analyze interactions between neurons become more accurate and more complex as the number of recorded neurons increases [7]. For spike -field entrainment analyses, it is also important to know where inputs that drive the cells are arriving from and how they interact locally. Cell types, dendritic arborization, local and long-range axonal projections all distribute unevenly across layers, shaping the integration and segregation of neural signals. Ultimately, we must decode laminar information across different structures and cortical regions to understand the spatiotemporal ensembles that represent cognition and memory [6].
The hippocampus is an example of a deep laminar structure, and is highly involved in encoding episodic memory formation presumably via spatiotemporal ensembles of neurons. Once multichannel silicon probes were developed, they quickly became a standard tool for hippocampal neurophysiological recordings in rodents. Increasing the density of laminar contacts in deep structures allows researchers to gain more information about these structures and more recently, study the potential effects of neurological disease states on single unit activity [8][9][10]. Although rodent cortical and hippocampal laminar structures were studied in great detail, the transition to large animal models has been slower. One of the major obstacles to addressing the research questions and approaches currently accessible to smaller-brained animals, are technical limitations in laminar silicon microelectrodes.
Much of the large animal brain is still inaccessible due to insufficient length of available silicon probes. In large animals and humans, the hippocampus is located deep inside the brain and is difficult to study electrophysiologically. Even many cortical targets cannot be reached with the current technology such as those typically utilized for brain machine interface (BMIs) such at the Utah electrode. In addition, laminar structure cannot be detected with a single wire or even tetrode recordings. Existing non-laminar solutions can scale to high density, but have single contacts at the tip, limiting the ability to differentiate single units as well as localize them within the laminar circuitry. These are therefore not suitable for understanding the relationships of neurons within laminar structures, as no information can be gathered about the relationships between neurons and their laminar inputs.
In order to reach deep regions of cortex and subcortical structures, silicon probes need to be long enough while still remaining narrow and balance flexibility to better match brain tissue properties to reduce gliosis and neuronal death. Technology necessary to produce finer features at longer length has only recently became available, including stitching across the die in the photomask process to utilize smaller features. New multichannel silicon probes designed for large animal electrophysiology allow for simultaneous recordings of fields and spikes from multiple layers of laminar structures such as cortex and hippocampus. In this study, we tested laminar silicon probes of various designs suitable for large animals and compared variables such as electrode size and coating material, probe thickness and width, and site arrangement and their effects on insertion ability, neurophysiological characteristics and biomechanical compatibility across probe designs.

Animals
Male Yucatan miniature pigs, which were chosen for their docile temperament, ease of handling, and slow growth rate for future chronic studies were used in this study. Pigs were purchased from Sinclair and underwent the current studies at the approximate age of 5 -6 months at a mean weight of 38 ± 3 kg (n = 17, mean ± SEM). We chose this age as the pigs are post-adolescent with near-fully developed brains, but young enough to be of a manageable weight for procedures and behavior [11][12][13]. All pigs were pair housed when possible, and were always in a shared room with other pigs. All animal procedures were performed in accordance with the University of Pennsylvania animal care committee's regulations, an AALAC accredited institution.
Surgical Procedure Yucatan miniature pigs were fasted for 16 hours then induced with 20 mg/kg of ketamine and 0.5 mg/kg of midazolam. Animals were intubated with an endotracheal tube and anesthesia was maintained with 2-2.5% isoflurane per 2 liters O2. Each animal was placed on a ventilator and supplied oxygen at a tidal volume of 10 mL/kg. A catheter was placed in an auricular vein to deliver 0.9% normal saline at 200 mL per hour. Additionally, heart rate, respiratory rate, arterial oxygen saturation, end tidal CO2 and rectal temperature were continuously monitored, while pain response to pinch was periodically assessed. All of these measures were used to titrate ventilation settings and isoflurane percentage to maintain an adequate level of anesthesia. A forced air warming system was used to maintain normothermia throughout the procedure.
All animals in this study underwent electrophysiological recordings and analysis as described previously [14]. Briefly, pigs were placed in a stereotactic frame, with the surgical field prepped and draped, and a linear incision was made along the midline. A 13-mm diameter burr hole was made centered at 7 mm lateral to the midline and 4.5 mm posterior to bregma, and the bony opening was subsequently expanded using Kerrison punches. Skull screws were placed over the contralateral cortex as a ground and non-probe reference signal. The dura was opened in a cruciate manner and the dorsal hippocampus was first mapped with high-impedance Tungsten electrode (AE = 125 µm, length = 60 mm, impedance = 0.5 MΩ; FHC, Cat# UEWSEGSEBNNM) in the sagittal plane, utilizing observed spiking activity to generate a twodimensional map. We then chose coordinates for insertion of multichannel silicon probes based on the map such that the spread of electrode contacts would span the laminar structure of the dorsal hippocampus at its maximal thickness in the dorsal-ventral plane, perpendicular to the individual hippocampal layers. The cortical recordings were done in the same track on the way to dorsal hippocampus. At the end of the procedure the anesthetic level was deepened (to 5% isoflurane) and pigs were sacrificed for histological analysis to locate tracks generated by probe insertion.

Neural Data Collection and Analysis
All electrophysiological recordings were made under isoflurane anesthesia (2 -2.5 %). Neural signals were amplified and acquired continuously at 32 kHz on a 64-channel Digital Lynx 4SX acquisition system with Cheetah recording and acquisition software (Neuralynx, Inc.).
Spike Detection and Analysis: Signals acquired from the tungsten monopolar electrode during the mapping procedure were bandpass filtered (600 Hz to 6 kHz) and thresholded for spike detection in real time. Thresholds were chosen based on observed signal to noise ratios during the session. Recorded spike trains and waveforms then underwent off-line automated spike sorting using KlustaKwik software (http://klusta-team.github.io/klustakwik/, RRID: SCR_014480) and further manual refinements using SpikeSort3D (Neuralynx, Inc.). Neural signals acquired from the 31 channels on the silicon probe were bandpass filtered (0.1 Hz to 9 kHz) in real time prior to sampling. Off-line spike detection and sorting was performed on the wideband signals using the Klusta package (http://klusta-team.github.io/klustakwik/, RRID: SCR_014480), which was developed for higher density electrodes, and manually refined with KlustaViewa or Phy (https://github.com/klusta-team/klustaviewa). The Klusta packages are designed to construct putative clusters from all probe channels simultaneously by taking advantage of spoke timing and the special arrangement of channels [15]. After manual refinement, resulting single-unit clusters were then imported into Matlab software, version R2017a for visualization and further analysis using custom and built-in routines (https://www.mathworks.com/products/matlab.html, RRID: SCR_001622). In order to minimize waveform shape distortion for visualization, the wideband signal was high pass filtered using a wavelet multi-level decomposition and reconstruction filter (level 6, Daubechies 4 wavelet) [16]. Waveforms were then extracted from this filtered signal and averaged for display. Spike trains were imported into NeuroExplorer software, version 6 for further analysis and to generate autocorrelograms and interspike intervals histograms (http://www.neuroexplorer.com/, RRID:SCR_001818).

Tissue Handling and Histological Examinations
Histological analyses were performed on brain tissue from male Yucatan miniature pigs in order to locate tracks generated by electrodes' insertion. While under 5 % isoflurane anesthesia, all animals underwent transcardial perfusion with 0.9% heparinized saline followed by 10% neutral buffered formalin (NBF). After post-fixation for 7 days in 10% NBF at 4°C, each brain was dissected into 5 mm blocks in the coronal plane and processed to paraffin using standard techniques. 8 µm sections were obtained at the level of the hippocampus in either the coronal or sagittal planes. Standard hematoxylin & eosin (H&E) staining was performed on all animals to visualize electrode tracks. For Luxol Fast Blue / Cresyl Violet (LFB/CV) Staining, tissue sections were dewaxed in xylenes and rehydrated to water via graded ethanols before being immersed in 1% LFB solution (Sigma, S3382) at 60°C for 4 hours. Excess stain was then removed by immersion of sections in 95% ethanol. Differentiation was performed via immersion in 0.035% lithium carbonate for 10 seconds followed by multiple immersions in 70% ethanol until the gray and white matter could be clearly distinguished. Slides were rinsed and counterstained via immersion in preheated 0.1% CV solution (Sigma, C5042) for 5 minutes at 60°C. After further rinsing, slides were differentiated in 95% ethanol with 0.001% acetic acid, followed by dehydration, clearing in xylenes and cover slipping using cytoseal-60.

Statistical Analysis
The data was analyzed using Graphpad Prism software, version 7 (http://www.graphpad.com/, PRID: SCR_002798). Single unit waveforms recorded with silicon laminar probes are presented as mean ± SD. Impedance of silicon probes are shown as mean ± SEM.

Results
Performance of multichannel silicon probes was tested acutely in Yucatan pigs during cortical and hippocampal recordings under anesthesia as described previously [14]. We have acutely tested silicon probes from various manufacturers and evaluated these probes based on their ability to record laminar oscillatory and single unit activities in cortical and hippocampal deep structures. Silicon probes were also evaluated based on their ability to cleanly separate distinct single-units.

Silicon Probe Design
Specific parameters of multichannel probes tested in this study are summarized in Table 1. Later version of silicon probes used in the study (P-trode w/ref, Vector, EDGE, and CAMB probes) were designed to have one channel substituted for a low-impedance reference site, which is placed 1-2 mm above the most-proximal probe site (depending on the probe design). The individual probe sites consisted of iridium oxide (IrOx) on NeuroNexus probes, platinum (Pt) on ATLAS probes, and gold (Au) on CAMB probes. 1 the most proximal site was substituted for a reference site with surface area of 4200 µm 2 ; 2 sites with alternating diameters, Ø = 35 µm and Ø = 25 µm; 3 an arrangement of four electrode sites placed close together; 4 spacing between individual tetrodes (at first site) as well as between linear sites; 5 similar to the linear layout, but electrode sites are strategically positioned at the edge of the substrate; 6 similar to the linear layout, but electrode sites are off-set by 21 µm relative to each other.
Multichannel silicon probes used in this study had a total of either 32 or 64 electrode sites and can be divided into sub-categories based on the site arrangement (Scheme 1). Most of the multichannel silicon probes used in pig hippocampal recordings were of linear design to resolve laminar structure (P-trode, ATLAS, EDGE and CAMB probes). For dorsal hippocampal targeting, this design provides 31 (or 63) channels for laminar recordings and results in the reference channel being positioned within the temporal horn of the lateral ventricle sitting just above the hippocampus (not shown).

Laminar Structure is Revealed by Site Spacing
The laminar structure of pig dorsal hippocampsus was examined electrophysiologically with various designs of silicon probes (Table 1). Laminar structure of the dorsal hippocampus shown in Figure 1 [14]). Since the hippocampal region of interest (at our standard recording coordinates in medial-lateral (ML) and anterior-posterior (AP) planes) is about 1600 µm, a silicon probe with 100 -200 µm site spacing will provide good coverage for the laminar structure with enough resolution to identify every layer (Figure 1). While the P-trode type silicon probe with 300 µm could be potentially used to examine the whole dorsal hippocampus (total length = 9600 µm), some of the layers could be easily missed (such as L-M). Similarly, while the ATLAS type silicon probe with 35 µm site spacing could be potentially used to analyze some hippocampal layers in more detail, it does not provide enough coverage to cover all hippocampal layers (total length = 1120 µm). Representative histopathological section is overlaid on laminar oscillatory activity to match representative layers as described previously [14].

Internal Reference Reduces Noise During Laminar Fields Recordings
Designing the probes with an internal reference for hippocampal recordings may help to reduce noise conferred from the electrical environment by using an electrically inactive but conductive structure such as the CSF of the ventricle in contrast to ECoG style skull screws. To estimate how an introduction of the internal reference affects the noise during hippocampal recordings, local field potetnial (LFPs) were recorded with EDGE style probe similar to the experiment displayed in Figure 1. Power specrum density (PSD) was then calculated at the lavel of pyramidal CA1 layer (Figure 2). During the same recording, the signals were first referenced to the skull screw similar to standard referencing protocol used in rodent and human recordings (Figure 2, red line). Afterwards, the internal reference on Edge probe was used, which eliminated most of the slow "drift" oscillations as well as 60 Hz frequencies presumably from AC noise (Figure 2, blue line). Frequencies in 120-300 Hz band ( Figure 2) and above were not affected. Since single unit activity usually lies within 600 -6,000 Hz frequency band, the choice of skull vs. internal reference will not affect the single unit recordings.

Figure 2.
Noise reduction with the internal reference on silicon probe. During the same recording, the signals were first referenced to the skull screw similar to standard referencing protocol used in rodent and human recordings. Power Spectral Density (PSD) of the channel presenting stratum radiatum was calculated in 2-300 Hz frequency band and then normalized by the its peak value and displayed in logarithmic scale ( Figure  2, blue line). Afterwards, the reference on EDGE style probe was switched to the internal reference. PSD was again calculated and normalized by the its peak value (Figure 2, red line). The internal reference eliminated most of the slow "drift" oscillations as well as most of the noise produced in 60 Hz-cycle frequencies (the peak at 60 Hz).

Sites Arrangement for Single Unit Detection and Sorting
In planar silicon probes, the individual electrode sites are placed on a face of the designs (P-trode, ATLAS, and Vector probes). In an attempt to decrease the damage near electrode sites and increase the exposure of the electrode sites to the parenchyma, we implemented a probe with EDGE design (Scheme 1a, panel 4). In contrast to P-trode and Vector style probes, the individual electrode sites on Edge probe are strategically positioned on the edge of the substrate, potentially reducing the interference of the insulator with the surrounding signals.
To characterize the laminar structure of the pig hippocampus, silicon multichannel probes with a linear design were used throughout the study. While the linear design is fit for a wide range of applications, isolation of single unit activity may be difficult with only a linear site arrangement. To utilize commercial spike sorting software such as SpikeSort 3D (Neuralynx), individual sites on silicone probes with a linear design were artificially grouped into sets of four (tetrodes) and then sorted. Single units recorded with 32channel probes of the earlier designs (P-trode, Scheme 1a, panels 1; and ATLAS, Scheme 1a, panels 2) could not be properly isolated with spike sorting software available at the time (Figure 3a), partially due to potential overlap between putative units between the artificial "tetrodes" (Figure 3a). A single unit (pink) was recorded on two adjacent sites, which got placed into separate tetrodes (the last site of top tetrode and the first site of bottom tetrode). As a result, this single unit was counted twice, creating a duplicate unit. To address this issue, we designed the Vector probe with custom arrangement of electrode sites (Scheme 1, panel 3). Most of the individual sites were arranged in four sites placed close together (tetrodes), allowing for high-quality cell discrimination in hippocampal recordings, with extra sites in between groups of tetrodes in order to maintain the laminar analyses. Four tetrodes were placed in pyramidal CA1 hippocampal layer (top part of the probe), while three tetrodes were placed in granular cell layer (bottom part of the probe). Single units recorded with Vector style probes and sorted with Spike Sort 3D software are shown in Figure 3b. The same unit could be registered on multiple sites of the same tetrode (blue color, first and last site) to facilitate spike sorting. As Vector style probes were created with the individual sites arranged closely to form tetrodes, more single units were isolated but cross-over of the units onto neighboring tetrodes were still observed occasionally in the CA1 layer due to the large size and dendritic arbor of these neurons (not shown). Since the advent of laminar spike-sorting software [15] the necessity no longer exists for specific tetrodes to be recreated in the silicon probes, however probe geometry is still an important part of the spike-sorting process (see section 3.2.3 below).
(a) (b) Figure 3. Sites Arrangement Defines Single Unit Separation. Single units were sorted with Spike Sort 3D software available from Neuralynx at the time. (a) Single units recorded with linear P-trode type probe (300 mm site spacing). A single unit (pink) was recorded on two adjacent sites, which got placed into separate tetrodes (the last site of top tetrode and the first site of bottom tetrode). In result, this single unit was counted twice, creating a duplicate unit. (b) Single units recorded with tetrode type Vector probe. The same unit is registered on multiple sites of the same tetrode (blue color, first and last site) which facilitates spike sorting by increasing the probability of isolating clusters.

Electrode Material Affects Spike Shape
To compare signal-to-noise quality as well as single unit separation, we compared standard electrode material irridium oxide (IrOx) used in rodent silicon probes with platinum (Pt) electrode material used in ATLAS style probes. Vector and ATLAS probes were tested acutely in the same animal by consequitive insertion in the cortical column at the same distance from the brain surface (dorsal-ventral (DV) coordinates). Single units were recorded and then sorted with Spike Sort 3D software for either IrOx Vector probe (Figure 4a) or Pt ATLAS probe (Figure 4b). While single units of various distinct shapes were sorted for Vector probe (Figure 4), all single units sorted for ATLAS probe did not have the standard action potential shape, suggesting a filtering effect for higher frequencies (Figure 3b). There was also a lower number of units picked up during recording with ATLAS probe, potentially due to the lower signal to noise ratio. The total of 45 single units were recorded and sorted in cortex and hippocampus (in the same one animal), with 18 single units on big electrode sites (Ø = 35 µm) and 27 small electrode sites (Ø = 25 µm) picking up unit activity. The size of an individual site doesn't appear to affect recordings of single unit activity as the ratio of units registered with both big and small sites in the same recording was about the same (40% and 60 % respectively).

Width and Thickness of Silicon Probes Defined by Production Methodology
The available silicon probe technology for large animals suffers from the limitations of low-resolution lithography, which determines the size and spacing of the lines (approx. 4µm features) and thereby limits the number of sites and increases the size of the probes (typical thicknesses of 50-100 µm and widths of 225-250 µm; Scheme 2). NeuroNexus large animal depth probes are 50µM thick, stiff and brittle, thus increasing insertion force but making them prone to fracturing during handling and decreasing yields beyond 10-15mm. This necessitated a novel interface cable solution whereby a cable is threaded through a 250µm stainless steel tube and bonded to the contacts at the end of the probe, allowing the probes to be 200 μA 100 ms 200 μA produced in lengths great enough to reach any part of the large animal brain. In addition, 32 individual sites are the current upper channel limit with these probes using contact lithography while maintaining reasonable widths. In order to produce longer silicon probes without the need for the extension tube, novel CAMB 64-channel flexible silicon probes are produced with newer lithography techniques with 1-4µm resolution features, and utilizing a new process that allows for stitching across the stepper lithography dies. This technology allows for a doubling of the number of sites on the probe, leading to the first 64-channel probes for large animals that are made exclusively from silicon wafer (Scheme 2, orange panel). Additional testing revealed that mock CAMB style probes with 35 µm thickness and 80 µm width (not shown) could be inserted into large animal brain with an 80µm width, giving a minimum dimension for the probes as stiffness increases linearly with width. Increasing a probe's width imparts greater mechanical stiffness of the silicon probe, which benefits insertion, but causes a further compliance mismatch between probe and brain tissue, potentially contributing to gliosis and local cell death for chronic implementations. By reducing the tissue impact during surgery, smaller, more flexible silicon probes may increase the survival rate of neurons in close proximity to the travel path.

CAMB 64-channel Silicon Probes
We also compared single-unit separation between a linear 64-channel CAMB probe and a 64-channel CAMB probe with a poly-2 design, which is similar to the linear layout, but electrode sites are off-set by 21 µm relative to each other (Scheme 1b). For the single units on the laminar probe, modern laminar spike sorting methods that classify units based on waveforms across all channels on the probe simultaneously (i.e. Klusta) lead to clusters with poorly defined classic separation measurement thresholds (i.e. L-ratio, Mahalanobis distance) for higher dimensions than tetrodes [15,20] To answer an open question whether non-linear probe designs are necessary for true resolution of single units using modern software such as Klusta/PHY, CAMB silicon 64-channel type probes were inserted in the pig cortex. Single units were recorded with CAMB 64-channel probes of either linear (Scheme 1b, panel 5) or poly-2 (Scheme 1b, panel 6) design, 30mm in length. Single units sorted with PHY software for both designs are shown in Figure 5.  Many separated clusters from the CAMB linear 64-channel probe recording contained multiple units, which were not well isolated from each other with current spike sorting methods. While some clusters from CAMB 64-channel poly-2 design probe contained multiple units, the proportion of multi-unit clusters was less when compared to the linear design (33% and 54% respectively) as shown in Figure 6. # sorted cells times more single units (green, n = 48) were identified. Also, the proportion of multi-unit clusters for poly-2 design was less when compared to the linear design (33% and 54% respectively).

Discussion
In this study, we implemented acute in vivo electrophysiology to characterize multichannel silicon probes from various manufacturers, and compared designs of a new higher density research probe. Silicon multichannel probes designed and used for large animal cortical and hippocampal electrophysiology allow for a longer area of coverage in contrast to single site recordings with tungsten electrodes, and also allowing for a slight drift in electrode placement in case of chronic implantations. They may also replace the need for multiple-insertion experiments if the animal is chronically implanted and/or a drive is utilized. Most of our silicon multichannel probes were designed with an internal reference to improve signal-to-noise ratio while leaving physiological slow oscillation like alpha and theta rhythms intact. In addition to having an extra reference signal to utilize or switch to in case of reference failure, it can also be used to compare the signals to the signals recorded with standard skull reference in rodent and human studies.
We tested various designs of silicon probes available from different manufacturers for a quality of single unit isolation. Linear, tetrode plus linear and poly-2 designs arrangement design were evaluated for single unit vs. multi-unit cluster isolation. Although most of the probes with a linear design site arrangement (P-trode, ATLAS, and EDGE style probes) were able to record neuronal activity of individual neurons, silicon probes with a tetrode and a poly-2 design site arrangement had better cluster separation with currently available sorting algorithms, as the proportion of single units to multi-units recorded with a given probe increased. While electrode material should be carefully considered as some materials affect shape of single units suggesting a filtering effect for higher frequencies, the site size did not appear to influence the number of single units recorded. Previous acute examinations in large animals have also noted the usefulness of the parallel geometry in isolating units, even prior to the advent of new sorting techniques [21].
We have also identified the advantages of minimizing probe dimensions such as shank width and thickness while increasing channel counts. Since the current commercially available probe technology employs contact lithography, it sets the resolution of patterned features to 4 µm and limits the number of electrode sites on the probe without increasing widths to those that will create impediments to long-term chronic recordings. Research probes described here (CAMB 64-channel) have been designed with maximal taper at the active electrode sites as well as placing at least one row at the edge of the probe is has been accomplished with the NeuroNexus EDGE probe. Application of newer photo-lithography techniques decreases the width/electrode ratio, minimizes local tissue damage, and increases the length that silicon probes can be produced for large animal work. Neuropathological examination of CAMB mock probes demonstrated that thinner probes produced less neuropathology, which is important point to consider for chronic probe implantations. Although a key benefit of long thin shanks is that they offer mechanical flexibility, thereby minimizing the mechanical compliance mismatch with the brain, long probes must remain stiff enough to penetrate the brain. By testing CAMB 25mm-long probes, we showed that 35µm shank thickness for 80µm wide probes strikes an ideal compromise between ease of tissue penetration and appropriate targeting vs. the need to maintain small device dimensions. This allowed for production and testing of the 30mm long 64-channel probes described above, which yielded substantially more units in the poly-2 style probe.

Conclusions
This study is a comprehensive survey of multichannel silicon probes designed for large animal electrophysiology available for use, as well as comparing them to a new design and process. NeuroNexus EDGE style probes were determined to be ideal commercially available probes for acute recordings from laminar structures due a linear site arrangement, 6mm of coverage, and the highest single-unit yields, with a potential for use in chronic implantations and awake behavior. In addition, cross-sectional area was found to be a crucial determinant of silicon probes' performance. Novel CAMB 64-channel probes with a poly-2 design were found to have an even better single unit separation and a denser sampling of the laminar structure than linear probes. By increasing channel density, we were able to better visualize laminar structure and to reduce tissue damage from probe insertion, identifying them as potential candidates for chronic implantations. Channel density, site arrangement, and physical profile of the silicon probe are all important factors to consider when designing probes for acute and chronic implantations to study laminar structures over time in future awake behaving animals. These results should be confirmed and expanded in chronic preparations, as it is likely that acute unit visibility immediately (<12 hours) post insertion is not predictive of chronic neural interface response and recovery post implantation. We hope these results will lower the threshold for adoption in acute and future chronic implantations, by demonstrating consistent yields in laminar electrophysiological recordings in large animals using commercially available probes, leading to new discoveries in both hippocampal and neocortical neurophysiology. In addition, better detection and understanding of laminar circuitry and changes in human disease (i.e. epilepsy, traumatic brain injury) are needed, and therefore viable electrodes need to be tested first in translational large animal models.


extfieldsofstudy: Biology Medicine

text:
Contemporary ancestor? Adaptive divergence from standing genetic variation in Pacific marine threespine stickleback

Background Populations that have repeatedly colonized novel environments are useful for studying the role of ecology in adaptive divergence – particularly if some individuals persist in the ancestral habitat. Such “contemporary ancestors” can be used to demonstrate the effects of selection by comparing phenotypic and genetic divergence between the derived population and their extant ancestors. However, evolution and demography in these “contemporary ancestors” can complicate inferences about the source (standing genetic variation, de novo mutation) and pace of adaptive divergence. Marine threespine stickleback (Gasterosteus aculeatus) have colonized freshwater environments along the Pacific coast of North America, but have also persisted in the marine environment. To what extent are marine stickleback good proxies of the ancestral condition? Results We sequenced > 5800 variant loci in over 250 marine stickleback from eight locations extending from Alaska to California, and phenotyped them for platedness and body shape. Pairwise FST varied from 0.02 to 0.18. Stickleback were divided into five genetic clusters, with a single cluster comprising stickleback from Washington to Alaska. Plate number, Eda, body shape, and candidate loci showed evidence of being under selection in the marine environment. Comparisons to a freshwater population demonstrated that candidate loci for freshwater adaptation varied depending on the choice of marine populations. Conclusions Marine stickleback are structured into phenotypically and genetically distinct populations that have been evolving as freshwater stickleback evolved. This variation complicates their usefulness as proxies of the ancestors of freshwater populations. Lessons from stickleback may be applied to other “contemporary ancestor”-derived population studies. Electronic supplementary material The online version of this article (10.1186/s12862-018-1228-8) contains supplementary material, which is available to authorized users.

Background
The ecological theory of adaptive divergence predicts that populations diverge phenotypically and genetically if they reside in distinct environments [1,2], potentially resulting in speciation. Some of the most striking examples of adaptive divergence come from species in which contemporary populations persist in environments putatively occupied by ancestral populations (e.g. [3][4][5][6][7][8][9][10][11]). To the extent that such populations have not undergone evolution, contemporary descendants of populations in ancestral environments can be used as a proxy for the ancestor, with phenotypic and genetic differences between these "contemporary ancestors" and "derived" populations being used to infer the direction, source, and pace of adaptation to derived environments. However, recent changes in the ancestral environment may stimulate evolutionary responses in the contemporary populations that inhabit it, complicating their utility as a proxy.
Standing genetic variation (SGV), defined as the variety of alleles segregating in a population [12,13], is expected to play an important role in parallel evolution. In particular, SGV permits rapid adaptation compared to de novo mutation, and increases the likelihood that the same beneficial allele will be present in different derived populations [14][15][16]. The role of SGV in adaptive divergence is readily measurable: if an allele fixed in the derived population is present in the contemporary ancestor at low frequencies, it likely contributed to adaptation [13]. However, the inference that an allele present in the contemporary ancestral population resulted in adaptation via SGV requires three assumptions. (1) The subset of individuals that originally colonized the derived environment must have contained the rare adaptive allele at some frequency; otherwise it arose from de novo mutation or subsequent gene flow. (2) The contemporary ancestor has undergone little evolution, including gene flow from the derived population (e.g., [17]). (3) The ancestral population has to have been properly characterized (i.e., population structure and allele frequencies associated with SGV). If there are multiple potential ancestral populations, each with a different pool of SGV, inference about the source and pace of evolution in the derived population will vary depending on which putative ancestral population is investigated (e.g. [12]). These assumptions must be verified to characterize accurately the role of SGV during population divergence.
Threespine stickleback (Gasterosteus aculeatus) provide perhaps the best documented examples of adaptation from SGV. Marine threespine stickleback occur widely in the northern hemisphere, including along the Pacific coast of North America from Alaska south to southcentral California. Across the north Pacific coast, much freshwater habitat formed recently (~10,000-20,000 years ago) in association with isostatic rebound following glacial retreat. The subsequent colonization of this habitat by stickleback allows tests of the significance of de novo mutation and SGV for adaptation (e.g. [17][18][19]). For instance, marine stickleback bodies are often covered by > 29 bony lateral plates, but fewer plates (0-10) have evolved in parallel in freshwater populations through selection on a rare marine allele [18]. Despite numerous studies indicating the role of SGV at either a single locus for platedness (Ectodysplasinhereafter Eda) or for multiple loci with unknown phenotypic effects [20][21][22], assumptions about the appropriateness of considering extant marine sticklebacks as representative of the ancestors of freshwater populations remains untested. Despite evident genetic variation in threespine stickleback among geographic clades [23][24][25][26], marine stickleback on the eastern Pacific are largely assumed to constitute a single population (e.g. [27][28][29][30][31]). This assumption is justified by the absence of barriers to gene flow in the marine environment [27], the migratory capacity of marine stickleback [32], the relative "evolutionary stasis" of marine stickleback inferred from the fossil record [28,29], and the low marine population structure reported from several local studies ( [20,33], but see [34]). Nevertheless, substantial evidence indicates local adaptation even in highly migratory marine fishes [35][36][37], and indeed in Baltic Sea threespine stickleback [38][39][40]. If Pacific threespine stickleback constitute a single population, their large population size should limit the effects of genetic drift and high gene flow should offset local adaptation [41]. Given these conditions, marine stickleback should not exhibit local differentiation that would generate regional differences in the initial colonists of freshwater lakes and streams. SGV should be the same along the Pacific coastand all freshwater populations could have evolved from the same initial pool of marine SGV facilitating parallel adaptation. These assumptions require formal testing to elucidate the role of SGV during adaptive divergence.
In this study, we consider phenotypic and genotypic variation of > 200 marine threespine stickleback from eight locations from California to Alaska to test hypotheses about the genetic structure of marine stickleback and its evolutionary consequences. Based on variation in plate phenotypes and genotypes associated with SGV at Eda, three-dimensional body morphology from micro-computed tomography (μCT) scans quantified using geometric morphometrics, and Genotype-by-Sequencing [42], we assess whether marine stickleback constitute a single population. By doing so, we test assumptions about the distribution of SGV in "contemporary ancestors" of freshwater stickleback. Additionally, we test predictions regarding the influence of SGV on the source of adaptation based on genomic sequences for a freshwater population from British Columbia. We specifically test the following null predictions: (1) Marine stickleback populations will not vary in the frequency and content (e.g. private alleles) of SGV. (2) Similarly, marine stickleback will not exhibit genetic population structure, which would otherwise influence the SGV regionally available for selection. (3) Marine populations will exhibit no phenotypic divergence in body shape or platedness [e.g. 18]. (4) If differences in SGV among marine populations occurs, genetic variation will show no evidence of having been shaped by natural selection. (5) If population structure in marine stickleback occurs, geographic proximity to a freshwater population will determine the extent of genetic divergence between marine and freshwater stickleback. (6) Differences in SGV among marine populations, if they occur, will not affect the candidate loci identified as contributing to adaptation in freshwater stickleback.

Methods
Threespine stickleback (n = 383, Table 1) were collected with minnow traps or seines during the summers of 2010 (Brannen Lake, British Columbia, hereafter BCFW), 2012 (Alaska) and 2013 (all other localities). Sampling locations extended along a 21.8 degree latitudinal spread (Table 1, Fig. 1), from California (south to north, CA01, CA02, CA03), through Oregon (OR01, OR02), the Puget Sound area of Washington (WA01), Vancouver Island (BC01) and Alaska (AK01). Locations varied in terms of benthos, freshwater input, and protectionfor instance, CA01 fish were sampled in a slough with freshwater input determined by precipitation, while OR02 were collected near a tidal gate close to the mouth of a river. Other marine species were collected alongside stickleback, such as bay pipefish (Syngnathus leptorhyncus) or smelt (Atherinops/ Atherinopsis sp.). Adults were captured in all localities with the exception of WA01, while OR01 contained a range of age classes. Stickleback were euthanized using buffered tricaine methanesulfonate (MS-222) or Eugenol (clove oil) and preserved in 70% ethanol. Fin clips were preserved in 95% ethanol for later sequencing. All collections were conducted in accordance with CCAC guidelines (AUP AC13-0040) and state/provincial/ national collection and import permits.
Sex was identified using primers developed by [43] that amplify sex-specific alleles at the idh locus. Alleles were visualized in a 2% agarose gel for 367 individuals.

Library preparation and analysis
Reduced representation DNA sequencing was used to generate Single Nucleotide Polymorphisms (SNPs) in order to assess population structure and adaptive divergence. Two hundred nanograms total genomic DNA was extracted per fish in January 2016 using Qiagen  [44]). Thirty to thirty-five individuals were included from each marine location, and 15 from BCFW. After digestion-ligation, fish were pooled into groups of nine. Cleanup and size selection were performed simultaneously using SPRI beads (Beckman Coulter), at a bead ratio of 0.8× and 0.61× for left and right-side cleanup, respectively. This left a fragment range between 250 and 600 bp. Pooled samples were divided into three technical replicates to ameliorate stochastic differences during PCR, and were amplified.
Replicates were pooled and left-side cleaned using SPRI beads. Pooled samples were quantified using a 2200 TapeStation (Agilent) and Qubit (Thermofisher) dsDNA high sensitivity assay. Equal volumes of each 2 nM pooled sample were pooled to make the final library. Library preparation followed the Illumina protocols for the Illumina NextSeq 500 Mid-Output kits with version 2 chemistry. Two sequencing runs were completed on the Illumina Next Seq 500 using 150 cycles and different final library concentrations -the first at 1.8 pM final concentration, the second at 1.1 pM final concentration. A 20% PhiX spike-in was used for both to compensate for the low diversity nature of the library. Results from the two sequencing runs were merged. Sequenced reads were cleaned and processed using Stacks v.1.35 [45,46]. Reads were de-multiplexed and cleaned using process_radtags, rescuing barcodes if the correction of a single sequencing error made them identifiable. GSnap [47] was used to align reads to the stickleback reference genome (Ensembl release 72 [48]), allowing for five mismatches with soft masking disabled. SNP calls were corrected in rxstacks using a bounded SNP model with an upper error rate of 0.1. Stringent filtering criteria were applied to the data, with slightly different filtering criteria used to address different questions. In order to determine population structure, each sampling site was assumed to constitute a distinct population. The filtering criteria for this "marine site" data set included: i) log likelihood threshold > − 60; ii) sequenced in more than 75% of individuals; iii) in 6 of 8 populations; iv) with minimum 4× coverage; v) minimum minor allele frequency of 2%, and vi) F IS > − 0.3. Individuals were included if they retained > 10,000 RADloci after cleaning (n = 239, Table 1). After population structure was determined, Adegent-recognized clusters rather than sampling location were used in the second Stacks run. This "marine cluster" data set used the same filtering criteria as above, with the exception that the variant needed to be present in all clusters. Finally, a "marine-freshwater" data set was run, which treated each sampling site as a distinct population and additionally included a freshwater population (see below). In this case the variant needed to be present in all eight marine populations and in the freshwater population. For each data set, all population statistics except for F ST were calculated using the populations module of Stacks.

Population genetic structure
Pairwise global and per-locus F ST were calculated using the Weir and Cockerham [49] adaptation implemented in hierfstat v0.04-22 [50] in R [51]; pairwise global F ST were tested for significance (> 0) using 999 permutations in GenoDive [52]. Discriminant Analysis of Principal Components (DAPC) [53] was used on the "marine site" data to assess population structure in the marine environment using Adegenet v2.0.1 [54], as it has shown to perform better than Structure under a stepping-stone model of dispersal [53]. The optimal number of Principal Components (PCs) to retain was calculated using both xvalDAPC and a-scores, which gave similar answers. The optimal number of clusters was assigned based on the lowest Bayesian Information Criteria (BIC) score using k-means clustering. As several possible k clusters had similarly low BIC scores, analyses were run and compared using 3 to 8 clusters.
An analysis of molecular variance (AMOVA), implemented in poppr v2.3.0 [55,56] using the Ade4 package [57], was used to determine the proportions of genetic variance among versus within sampling sites or Adegenetrecognized clusters. Missing values were replaced with the average frequency for a locus; ignoring missing values did not alter overall patterns. To explore the possibility of cryptic population structure, each sampling site was further analysed individually using Stacks and Adegenet.
The distance between each sampling site was measured as distance along the coast (km) using Google Maps. Distances were measured to or from the mouth of each bay. Neighbouring localities were separated by 242-479 km, except for BC01-AK01, which were separated by approximately 2500 km of coastline. The location of WA01 in Puget Sound resulted in all locations south of Washington being closer to BC01 than they were to WA01. Genetic distance was calculated using the pairwise global Weir and Cockerham F ST measures from hierfstat. Geographic and genetic distance matrices were compared using a Mantel test from the Adegenet package with 999 replications to determine Isolation-by-Distance (IBD).
Population statistics were also estimated for the "marine cluster" data set in Stacks using the optimal Adegenet-identified clusters.
A phylogenetic network was calculated using SNPhylo [58] and visualized using FigTree v.1.4.3 [59]. The "marine sites" data set was used, but SNPhylo additionally filtered loci based on linkage disequilibrium. Individuals were colour-coded according to their recognized genetic cluster. A hierarchical clustering tree was additionally constructed using BayPass v.2.1 [60].

Platedness
Plate variation among populations was assessed using plate number and Eda genotype. Adult stickleback (i.e. fish > 30 mm standard length (SL)) (n = 281, Table 1, Additional file 1: Table S2) were stained in Alizarin red. Plate number, including keel, was counted on both sides of the body and summed. Low-plated individuals without keel (LPNK) were defined as individuals with < 20 anterior plates. Partially-plated keeled (PPK) fish had 21-59 plates, including at least one plate at the caudal keel. Fully-plated keeled (FPK) fish had ≥ 60 plates. Additionally, some low-plated fish had a keel (LPK) and were defined as having < 20 anterior plates plus additional plates at the caudal keel. Partially plated stickleback that lacked a keel (PPNK) had > 20 anterior plates but had no plates at the caudal keel. Individuals were also genotyped at the Stn382 locus [18] as this microsatellite is linked to an indel in intron 1 of the Eda gene, yielding a 218 bp "fully-plated" allele (C) or a 158 bp "low-plated" allele (L) [61]. Genotyping followed the protocol of [43]. Individuals were genotyped as LL (homozygous for the low-plated allele), CL (heterozygous), or CC (homozygous for the fully-plated allele). This approach allowed juveniles (< 30 mm SL) to be included in the analysis (total n = 361), and provided genetic information at a locus with known adaptive significance that was not recovered from sequencing. Hardy-Weinberg equilibrium (HWE) was assessed for Stn382 for each marine site using a goodness-of-fit Chi-squared test.

Morphometrics
Phenotypic variation among populations was further assessed using morphometric analysis. Stickleback > 30 mm SL that had been preserved with relatively little bending (n = 272) were straightened, and spines and fins held flat against the body using plastic wrap. μCT scanning at a resolution of 20 μm was conducted in a standardized fashion for all individuals using a Scanco μCT35 instrument (Scanco AG). Three-dimensional images were generated from the anterior point of the premaxilla to the posterior tip of the pelvic spine, using standardized isosurface thresholds in Amira 5.4 (FEI Visualization Sciences Group). Fifty-five landmarks were plotted on the left side of each fish (Additional file 1: Table S1, Fig. 2) and raw landmark scores were exported to MorphoJ v1.06a [62] for further analyses. A prior study had removed the operculum on the left side of all AK01 stickleback, so landmarks were plotted on their right sides. Data were first transformed to remove differences associated with isometric scaling, rotation and translation using Procrustes superimposition. Residuals from a within-marine site multivariate regression on centroid size were estimated and used in all subsequent analyses. Principal Components Analysis (PCA) determined the major axes of phenotypic variation. Canonical Variate Analysis (CVA) was used to determine Procrustes distances using sex and marine site as categorical variables, although BC01 had only one female and OR01 included one individual of unknown sex. The significance of Procrustes distances among pairwise comparisons of marine site-sex combinations was determined based on 10,000 permutations and a corrected α of 0.0005. Discriminant Function Analysis (DFA) was used to determine the likelihood that individuals could be reassigned to their site of origin, given their phenotypes. For this analysis the effect of sex on reassignment success was not assessed.

Selection on phenotypic variation
Patterns of phenotypic variation among populations may be the result of genetic drift, natural selection, or phenotypic plasticity. One way to rule out neutral evolutionary processes is to compare estimates of phenotypic divergence (P ST ) with neutral expectations based in part on observed neutral genetic differentiation. Observed phenotypic divergence was estimated as: where σ 2 B and σ 2 W were the between-and withinpopulation components of variance, respectively, for plate count and the first four PCs from the morphometric analysis (as per [63,64]). Variance components were estimated for all marine sampling sites together (global P ST ) and pairwise using lme4 [65], with population as a random effect. Genetic divergence at the Stn382 locus for Eda (F STQ ) was estimated using the Weir and Cockerham method in Genepop V4 [66]. Neutral genetic divergence (F ST ) was estimated in hierfstat using non-genic SNPs identified from our data set using Biomart [67]. Non-genic SNPs may still be linked to loci under selection, so this approach provides a conservative estimate of neutrality.
Selection was inferred based on two methods. The first assessed the association between P ST -F ST and F STQ -F ST using Mantel tests. This measure is based on the expectation that phenotypic or QTL divergence will be uncorrelated with neutral genetic divergenceby extension implicating selection to explain such patterns. The second test involved Whitlock and Guillaume's [68] method using the R-code from [69]. In brief, the expected between-population variance component for a neutral phenotype was estimated by using observed non-genic F ST and observed within-population variance component for the phenotype: As per [69], the distribution of neutral σ 2 B was estimated by generating a χ 2 distribution with six degrees of freedom (one less the number of sampling sites excluding Washington), and multiplying a randomly drawn value from this distribution by σ 2 B . From this new distribution expected neutral σ 2 B were drawn 10,000 times and used to create a distribution of neutral P ST -F ST . The observed P ST -F ST was then compared to this distribution and the quantile of the neutral distribution that lay beyond the observed value was used as the probability of the observed outcome in the absence of selection, p. Under the expectation of no selection, p is > 0; selection is evidenced if p = 0. F STQ -F ST values was also compared to the neutral P ST -F ST , as per [40].
To quantify the relation between phenotype and Stn382 genotype, a generalized linear model (GLM) was fit to the data in R, using the glm routine, as per [40], with plate number as the dependent variable, genotype as a fixed effect, and using a log-link function with a quasi-Poisson error distribution. Furthermore, a Mantel test was used to estimate the correlation between pairwise F STQ and P ST measures.

Selection on genetic variation in the ocean
Under the assumption that marine stickleback populations have a shared history, the covariance matrix of population allele frequencies (Ω) was estimated in BayPass [60] using the "marine site" data. From this a hierarchical clustering tree [60] was generated, assuming no gene flow. A covariate-free genome scan was then performed to identify outlier loci putatively under selection, using per-locus measures of differentiation (XtX). The simulate.baypass function was used to estimate the posterior predictive distribution of XtX using a pseudo-observed data set (POD) [60]. Any loci in the "marine site" data set with XtX values above the POD-estimated threshold were scored as outlier loci potentially under selection. Genic outliers were identified using BioMart [67].

Marine-freshwater genetic divergence
To assess the extent to which the choice of putative "contemporary ancestor" affected inference about adaptation to fresh water, pairwise F ST was estimated for each marine-freshwater pair using hierfstat. In BayPass, covariance matrices, PODs and XtX thresholds were estimated for each marine-freshwater comparison. Outlier loci were examined to determine if the same outliers were being consistently recovered irrespective of the origin of the marine fish.

Sequencing results
Over 192 million reads passed initial filters (Additional file 1: Table S3) in the "marine sites" data set. Two hundred eighty-two RAD-loci were excluded due to excess heterozygosity. Filtering minor alleles at a threshold of 2% reduced the number of retained loci by 32%. After filtering, between 230,010 (AK01) and 426,018 (OR01)  Table S1 for identity of landmarks loci were retained for each site sampled, generating between 1877 (AK01) and 5204 (OR01) SNPs (Table 2), for a total of 6655 variant loci.

Standing genetic variation
All marine samples exhibited SGV, ranging from an average of 0.82% (AK01) to 1.22% (OR01) of the total SNPs genotyped in a given population; however, the pool of SGV varied from California to Alaska (Table 2). Stickleback from each marine location contained multiple private alleles (alleles found only at that location) (Fig. 3) and were polymorphic for a portion of the variant loci (loci that were polymorphic in at least one marine site). Polymorphism among variant loci varied from 57% (AK01) to 80% (OR01). For variant loci, the average frequency of the major allele (present in > 50% of all sequenced stickleback) ranged from 88% (OR02) to 93% (AK01) (Additional file 1: Figure S1), suggesting that the frequencies of SGV also differed among locations. Heterozygosity ranged from 0.10 (AK01) to 0.17 (OR02) for variant loci. Population-level average F IS over all variant loci varied from 0.032 (CA01) to 0.064 (OR01) ( Table 2). Although F IS was close to 0 for most loci, it approached 1 for a few loci (Additional file 1: Figure S2).
Significant population genetic structure was detected. The best supported number of clusters from the eight marine locations sampled was five (BIC = 1379, Fig. 4, Additional file 1: Table S4). The five clusters were, from south to north, CA01, CA02, CA03-OR01, OR02, and WA01-BC01-AK01. The CA03-OR01 cluster also contained seven individuals from OR02 and a single individual from AK01; otherwise individuals clustered with others from their sampling locality. The possibility of a single genetic cluster was as well-supported as ten genetic clusters (BIC = 1392). Three to six clusters had BIC values that differed little from the best-supported model. Altering the number of putative clusters revealed different population structures (Additional file 1: Table S4), with WA01 and AK01 continuing to cluster until k = 13. Cryptic population structure was not evidenced for any sampling locality (Additional file 1:   Table S5). Even when WA01, BC01, and AK01 were included in a single analysis, k = 1 was the best supported cluster (BIC = 436.9). However, k = 2 (BIC = 438.2) and k = 3 (BIC = 440.2) still separated individuals by locality. The genetic variance partitioned between clusters by AMOVA was low (19%) compared to within clusters (81%). When clusters were not considered, 33% of variation occurred between marine sites. A Mantel test of pairwise geographic distances and pairwise Weir and Cockerham F ST was non-significant when all sites were included (r = − 0.2, p = 0.8). This was largely driven by the extreme distance between the WA01, BC01, and AK01 sites. If AK01 was excluded from the analysis the association between geographic and genetic distance was weakly significant (r = 0.5, p = 0.02) (Additional file 1: Figure S3).
The phylogenetic network largely agreed with Adegenet assignment (Fig. 5). The network revealed greater intermixing of groups than did Adegenet, but WA01, BC01, and AK01 still largely clustered together and comprised a separate lineage from most southern  Table 1 for label meanings  Table S2).

Morphometrics
Phenotypes varied extensively among sites. The first eight Principal Components explained 71% of all phenotypic variance (Additional file 1: Figure S5). The first two Canonical Variates (CVs) explained 57% of the variation among combinations of site and sex, and the first four CVs explained 77% (Additional file 1: Table S8). CV1, after accounting for differences in centroid size, revealed that BC01 fish had narrow, streamlined bodies with dorsal and ventral landmarks both shifting inward relative to the consensus fish (Figs. 7 and 8). Californian fish were grouped close together on CV1 and had squatter, less streamlined bodies with dorsal and ventral landmarks shifted away from one another relative to the consensus. AK01, OR01 and OR02 had intermediate phenotypes between BC01 and California. CV2 showed a gradual transition from CA01 to BC01, but AK01 was clearly distinct from all other sites along this axis. AK01 showed substantial dorsolateral and anterior-posterior constriction of the body relative to all other sites (Figs. 7 and 8).
The sexes differed morphologically at all sites except AK01 (Additional file 1: Table S9), but the sexes still largely grouped together according to sampling location. Fig. 4 Adegenet-identified clusters for k = 5. Inset shows hypothetical range of each cluster. Note that the cluster identified as CA03, OR01 contains one AK01 and seven OR02 individuals. See Table 1 for label meanings CA02 and OR02 were exceptions, with males from both sites clustering with OR01 males. Similarly, CA01 and CA02 females had morphologies that were not significantly distinct.

P ST -F ST and F STQ -F ST comparisons
P ST was estimated as 0.46 for platedness, 0.60 for PC1, 0.29 for PC2, 0.23 for PC3, and 0.09 for PC4. Fig. 5 a Result of the phylogenetic analysis using SNPhylo. Individuals are colour-coded according to their five Adegenet-recognized clusters (Fig. 4). b Hierarchical clustering tree using BayPass. See Table 1 for label meanings Fig. 6 The frequency of different Eda genotypes using the Stn382 marker, for each sampling site. "The North" refers to samples from Washington, British Columbia, and Alaska. CC = homozygous for the fully-plated allele. LL = homozygous for the low-plated allele. CL = heterozygous F STQ was 0.60. Plate P ST and F STQ greatly exceeded the range of the neutral P ST -F ST distribution (p = 0 for both), as did P ST for PC1 (p = 0). P ST for PC2 was marginally significant but within the tail of the neutral distribution (p = 0.002), while P ST for PC3 (p = 0.02) and PC4 (p = 0.7) were well within the neutral distributions (Fig. 9).
Of all of the Mantel tests between P ST , F ST , F STQ , and distance, only one was significant and two were marginally significant with corrected α = 0.05/18 = 0.0027: PC1 P ST correlated positively with geographic distance, plate count P ST correlated positively with F STQ , and, surprisingly, PC1 P ST correlated positively with F STQ (Table 4, Additional file 1: Figure S6 and S7). The relation between plate P ST and F STQ was substantiated with a generalized linear model which showed a decrease in plate number with number of L alleles (null deviance = 4813 with 277 d.f., residual deviance = 349 with 275 d.f., p < 0.0001, Hosmer and Lemeshow goodness of fit text:

Selection on genetic variation in the ocean
102 of 6655 loci from the "marine sites" data set were flagged as outliers, using an XtX threshold of 15.52 (Additional file 1: Figure S8). Although variant loci were sequenced across all 21 chromosomes and an additional 65 scaffolds, outliers were only detected on 16 chromosomes and two scaffolds. Of these, 8 of 36 (22%) variant loci on scaffold 37 were outliers, followed by 13 of 178 (7%) on linkage group (LG) XXI, 11 of 184 (6%) on LGXIX, and 18 of 498 (4%) on LGIV. Of these 102 outliers, 20 were located within 16 genes (Table 5), although none of these genes have been previously studied in stickleback.

Genetic differentiation and outlier analysis for marinefreshwater comparisons
The "marine-freshwater" data set identified 132,415 loci, of which 1912 were variant. As expected, marine-freshwater divergence was high and in all but one instance was > 0.25 ("very great" differentiation).  Table S11). The maximum difference in pairwise per-locus F ST estimates for a single locus was 0.991, and the average difference was 0.19. 314 (17%) SNPs had a minimum F ST estimate of little genetic differentiation (< 0.05) in at least one marine-freshwater contrast, but great genetic differentiation (> 0.25) in another.
97 of 1912 loci were identified as outliers in at least one marine-freshwater comparison, using XtX thresholds of 4.9 (CA03 -FW) to 6.0 (CA01 -FW). Seventeen outlier loci resided in 13 genes (Table 6). Only 3 of 97 loci were flagged as outliers in all eight marine-freshwater comparisons, while nearly a third (31%) were flagged in only a single comparison. Five of the genic loci were outliers in seven comparisons (tub, S100P, and three novel genes), and five were outliers in only one (including a different S100P locus).
Surprisingly, 12 of the loci flagged as outliers in the "marine sites" data set were also flagged as outliers in at least one marine-freshwater comparison (Table 6), including a novel gene on LGII (n = 7 comparisons), and csnk1g2a (n = 3 comparisons). The gene OVGP1 was also flagged as containing outlier loci in both analyses, but different loci were flagged in each case.  Table 1 for label meanings

Marine stickleback exhibit between-population genetic variation
The significance of SGV for parallel evolution depends on its occurrence in the ancestral populations. Marine stickleback harbour SGV, and it differs between populations. Differences in gene expression between two marine BC populations suggested this possibility [10], but here it has been quantified across an extensive latitudinal range. The extent of SGV, 0.8-1.2% of all sequenced loci, is intermediate to that reported from other studies [34,70]. Nucleotide diversity varied from 0.0016 to 0.0027, consistent with results from Alaska (0.0022 and 0.0025 [20]) and slightly lower than that reported from Oregon (0.003-0.0036 [34]). All marine locations harboured some degree of private alleles, even after ignoring minor alleles at < 2% frequency, suggesting that not just frequencies of SGV but also content of SGV can differ from site to site. Furthermore, the superiorly. a CV1 for a BC01-type body shape; b CV1 for a CA01-type body shape; c CV2 for an AK01-type body shape; d CV2 for a BC01-type body shape. Light blue wireframe shows the consensus morphology, while dark blue shows the conformational change best-known example of SGV, Eda, was present at varying frequencies between populations and is likely under selection in the marine environment. Such variation in the content and frequency of SGV, in turn, led to compelling evidence for population genetic structure. Marine threespine stickleback showed substantial population genetic structure along the Pacific coast of North America. Although F ST values (average F ST = 0.088) were generally lower than those reported for marine-freshwater divergence (e.g. [20,34,71]), they were higher than those reported for other marine stickleback populations along the North American Pacific coast (two Alaskan populations: F ST = 0.0076 [20]; three Oregonian populations: F ST = 0.007 [34]). However, they align with studies from Europe [39,40,72].
Five genetic clusters were identified for the eight sampled localities, although structuring was hierarchical, with some clusters more genetically diverged than others. The most widespread cluster occupied > 2700 km of coastline from Washington to Alaska. Marine stickleback from this genetic cluster have been well characterized, with genetic divergence reported to be low between populations separated by up to 1000 km [20,33,73]. Such low structuring between proximate marine populations has led to a basic assumption in stickleback literature that marine stickleback exhibit little genetic or phenotypic diversity globally (e.g. [27-29, 64, 74])an assumption supported by the iconic image of distinct freshwater stickleback forms radiating from a single fully-plated marine stickleback type (e.g. [75]). In contrast, results from a broad range suggest that such generalizations should be restricted to the northern genetic clusterand even it contains morphological and genetic differentiation that could be adaptively significant.
The southern genetic clusters were sequentially separated by a few hundred kilometres, well within the migratory ability of marine stickleback [76][77][78]. IBD was evident only after removing AK01 from the dataset, suggesting that limited migration could explain patterns of divergence between the southern genetic clusters. However, IBD needs to be interpreted with caution, as geographic distance was correlated with latitude, and latitudinal variation can be associated with environmental clines [79]. Whatever the causes that shape genetic variation between stickleback populations, the distribution of SGV among different marine populations affects inference about the source and pace of selection in the freshwater environment, and complicates attempts to uncover loci that are under selection in derived populations.

Marine stickleback exhibit between-population phenotypic variation
Marine stickleback are generally considered to be fullyplated (e.g. [27]), yet the Eda genotype for low-platedness has an ancient marine origin [18]. The low-plated allele has been hypothesized to exist in the marine environment as SGV only when transported from the freshwater environment [17] or when masked by marine modifying alleles [18]. If the low-plated allele exists at low frequencies in the ocean, behaviours that facilitate the movement of low-plated marine stickleback into fresh water could also account for the consistent colonization of rare low-plated stickleback in lakes and streams [80]. We found substantial variation in the frequency of the low-plated allele, to the point that it was the major allele in some Californian and Oregonian populations, but was absent from BC01. This finding is consistent with other records of low-plated marine stickleback in California ( [81], but see [18,82]) and high frequencies of low-platedness in European marine stickleback [61,74,[83][84][85]. However, it contrasts with other Pacific and Atlantic North American studies that focussed on northern sites [86][87][88][89][90][91].
The observed SGV at Eda could affect the rate at which adaptation to lakes occurred in the past. Indeed, it may explain why reduced plate size has evolved in some freshwater populations, rather than reduced plate number, as an alternative strategy that may have been required in the absence of SGV at Eda [92,93]. Such variation at Eda is a particularly striking reminder that the function of full-platedness in marine stickleback remains unknown (for an example of its possible use, demonstrated in freshwater populations, see [94]).
Marine stickleback body shape also varied extensively along the coast. Californian populations tended to have squatter body shapes that appeared to be less streamlined than their northern counterparts. The functional Table 4 The observed correlation and p-values (p) for Mantel tests between geographic distance, neutral genetic distance (F ST ), phenotypic distance (P STfor plates or the first four Principal Components (PCs) of body shape), or genetic distance at Eda (F STQ )

Comparison
Observed  LG Linkage Group. XtX provides an average if the number of SNPs is > 1. The XtX threshold, above which a SNP was considered an outlier, was 15.52  Table 6 Pairwise XtX values and outlier locus information for each of the eight marine-freshwater comparisons.
LG is linkage group, although scaffolds are also included. Position refers to the nucleotide position along the linkage group (Continued) LG significance of these differences requires testingbut it is interesting that the streamlined fish were from a single genetic cluster, while the squat Californian fish exhibited significant population structuring between neighbouring localities. This potentially indicates extensive migration along the northern coast that is not mirrored in the south. Jamniczky et al. [95] reported considerable morphological divergence between neighbouring sampling sites in British Columbiagroups presumably with little to no genetic divergence, implicating plasticity as a driver of morphological variation. Morris et al. [79] similarly reported variation in vertebral number and standard length with latitude. However, two related analyses suggest that selection may also play a role in shaping phenotypic diversity. Based on DeFaveri and Merilä's [40] method, Pacific coast stickleback exhibited selection for platedness similar to that observed for Baltic Sea stickleback. There was also suggestive evidence for selection on PC1 of body shape, which largely corresponded to CV1more streamlined northern fish, more squat southern fish. The role of plasticity in affecting these results remains to be determined. PC1 was marginally associated with Eda genotype. To our knowledge, this is the first study to demonstrate pleiotropic or linked effects of Eda on body shape in the marine environment (for marine-freshwater or freshwater-only evidence for pleiotropy, see [96]; for other forms of phenotypes associated with Eda see [80,[97][98][99][100][101]), which could result in low-platedness being an indirect target of selection in some marine habitats. It is possible, for instance, that a pleiotropic relationship exists between body shape, Eda, and thermal tolerancea possibility suggested by the relationship between low plate frequency and latitude in anadromous populations of Europe [83]. This is one of several possible explanations that requires formal testing.
Although Eda is the best-characterized example of SGV, the outlier analysis revealed other potential candidate genes for selection in the marine environment. Extensive differentiation was evidenced at some loci. Local adaptation despite gene flow has been found in other marine fish populations [35,[102][103][104], including European marine stickleback [72]; but the extent to which stickleback south of Washington exhibited population structure was unanticipated.
Smaller bodies and distinct body shapes tend to evolve in freshwater stickleback populations [105,106], often with significant correlations between morphology and freshwater biotic and abiotic factors [107]. Given the morphological variation among marine stickleback, elucidating whether freshwater morphology is the result of plasticity, SGV, or de novo mutation requires informed decisions about what constitutes the marine ancestor.

Inferring the source and pace of adaptation in freshwater stickleback
The choice of marine stickleback affected inference concerning the source and pace of adaptation in one freshwater population from Vancouver Island. Few outlier loci were consistently recovered in all marine-freshwater CA01 CA02 etc. identifies the marine group that was compared to the Brannen Lake British Columbia freshwater population. The number in brackets () after each population identifier is the XtX threshold value used for that marine-freshwater comparison above which a locus was flagged as an outlier. If the locus was an outlier in that comparison its XtX value is given; otherwise the cell is blank. Total refers to the number of comparisons for which that locus was flagged as an outlier. If Marine = YES that same locus was also flagged as an outlier in the marine environment comparisons; had only a single marine population been used in this study, the outliers reported would differ depending on which marine population had been chosen. Many studies involve comparisons between geographically proximate marine and freshwater stickleback pairs (e.g. [21,22]), presumably to account for the possibility of population structure in the marine environment. Yet F ST was the lowest when the freshwater population was paired with a geographically distant population from northern Oregon, a finding that is difficult to reconcile with the assumption that the nearest marine population is the most suitable ancestral type. The occurrence of Japanese mtDNA haplotypes in Haida Gwaii lake populations that are not present in Haida Gwaii marine populations [108,109] suggests that this may not be an isolated incident. Similarly, several studies have noted, but not explained, the fact that northern marine stickleback are genetically more similar to southern than northern freshwater stickleback [17,21,22]. Nearly one third of outlier loci were only outliers for a single marine-freshwater comparison. Thus, one's choice of marine population could produce spurious inferences about the loci under selection, or miss true candidate genes. Clearly more information is needed going forward about the relationship between marine and freshwater stickleback.
Inferring the role of SGV during adaptive divergence in ancestral-derived comparisons requires three conditions to be met, for which we have provided varying evidence. The first, that colonists had to contain the same variants as present in the ancestral population, could not be tested here. Second, the contemporary ancestor has to have been relatively evolutionarily static. Both phenotypic and genetic evidence demonstrates that this condition does not hold for marine stickleback. Contemporary marine populations have diverged phenotypically in a manner beyond neutral evolutionary expectations, and genetically in a way partially explicable by selection. This means that contemporary marine threespine stickleback populations are genetically and phenotypically distinct from their own ancestorsand it was these ancestors that also originally colonized lakes and streams along the coast. Thus the term "contemporary ancestor" is a misnomer, as contemporary marine threespine stickleback populations do not reflect the ancestral condition. Interpretations of freshwater stickleback evolution need to be tempered by marine stickleback evolutionary history.
Third, the ancestral population has to have been properly characterized. Eastern Pacific stickleback exhibit some population genetic structure, although consistent with other reports [20,33] stickleback north of Oregon constitute a single genetic population. This means that along much of the coast freshwater environments were likely colonized by distinct marine stickleback populations, which differed in SGV frequency and content. Furthermore, it is likely that marine stickleback have exhibited range contractions and expansions along the southern and northern coasts throughout their evolutionary history, most recently in the north after the last glacial retreat [110]. This means that there is no a priori reason to expect that a marine population currently proximate to freshwater populations are descendants of the ancestors of those freshwater populations. Population structuring and evolutionary history thus changes our understanding of the ancestral condition of marine stickleback, and requires that we carefully consider the use of contemporary marine populations when addressing evolutionary questions.

Conclusions
Studies that compare marine and freshwater stickleback may need to adjust their methodologies in light of marine stickleback variation. The typical image in textbooks is of a single marine stickleback form from which numerous freshwater forms radiate [75]. Our data suggests that there is phenotypic and genetic variation in marine threespine stickleback which likely impacts freshwater stickleback diversificationto assume a single marine form is no longer tenable. Furthermore, assuming ancestral status for the marine population most geographically proximate to the freshwater population of interest is problematic, unless it can be directly demonstrated (e.g. [111]). So where does this leave the comparative method? One possibility would be to reconstruct the genotype of the ancestor to all eastern Pacific marine sticklebackbut this would ignore the important role that local variation has played in freshwater stickleback evolution. Another possibility could be to conduct larger-scale geographic sampling than has heretofore been done, of both marine and freshwater forms, in order to determine a more thorough evolutionary history of this species. Then, having taken evolutionary relationships into account, comparisons can be made using better-justified "contemporary ancestors".
"Contemporary ancestors" are used in a number of systems [3][4][5][6][7][8][9][10][11] for addressing evolutionary questions. They are particularly useful for determining the role of SGV during evolution, and for identifying the alleles involved in adaptation to new environments. Clearly, care must be exercised in characterizing these proxies of the ancestral form, as unaccounted population structure and current evolution can lead to spurious interpretations of adaptation. Whether the lessons from stickleback apply to other species with smaller geographic distributions or more limited opportunities for gene flow waits to be seen.

Additional file
Additional file 1: Additional tables and figures. This document contains information supporting the main text, including: the list of landmarks used for 3D morphometrics, population-specific details on sex and platedness, additional information pertaining to methodology, the distributions of major allele frequencies and FIS per population, isolation-by-distance analyses, details regarding Adegenet-recognized clusters, CVA and DFA morphometrics results, additional Mantel tests, global pairwise FST values for each marinefreshwater comparison, and results pertaining to the outlier analysis.