s2fieldsofstudy: Biology Medicine

text:
The metabolomic plasma profile of myeloma patients is considerably different from healthy subjects and reveals potential new therapeutic targets

Introduction Multiple myeloma (MM), a malignant plasma cell disorder, is still an incurable disease. Thus, the identification of novel therapeutic targets is of utmost importance. Here, we evaluated the peripheral blood-based metabolic profile of patients with MM. Material & methods Peripheral blood plasma levels of 188 endogenous metabolites, including amino acids, biogenic amines, acylcarnitines, glycerophospholipids, sphingomyelins, and hexoses were determined in patients with plasma cell dyscrasias: monoclonal gammopathy of undetermined significance, a precursor stage of MM (MGUS, n = 15), newly diagnosed MM, (NDMM, n = 32), relapsed/refractory MM (RRMM, n = 19) and in 25 healthy controls by mass spectrometry. Results Patients with NDMM, RRMM and MGUS have a substantially different metabolomic profile than healthy controls. The amount of eight plasma metabolites significantly differs between the NDMM and MGUS group: free carnitine, acetylcarnitine, glutamate, asymmetric dimethylarginine (ADMA) and four phosphatidylcholine (PC) species. In addition, the levels of octadecanoylcarnitine, ADMA and six PCs were significantly different between RRMM and MGUS patients. 13 different concentrations of metabolites were found between RRMM and NDMM patients (free carnitine, acetylcarnitine, creatinine, five LysoPCs and PCs). Pathway analyses revealed a distinct metabolic profile with significant alterations in amino acid, lipid, and energy metabolism in healthy volunteers compared to MGUS/MM patients. Conclusion We identified different metabolic profiles in MGUS und MM patients in comparison to healthy controls. Thus, different metabolic processes, potentially the immunoregulation by indoleamine 2,3 dioxygenase-1 (IDO), which is involved in cancer development and progression supporting inflammatory processes in the tumor microenvironment and glutaminolysis, can serve as novel promising therapeutic targets in MM.

Introduction
Multiple myeloma (MM) still remains an incurable malignancy of antibody-producing plasma cells. Monoclonal gammopathy of undetermined significance (MGUS) is a common asymptomatic precursor state of MM. Process of myeloma initiation is mediated by interaction of inherited genetic and environmental factors [1]. Clinical manifestations of MM occur due to proliferation of malignant plasma cells or effects of aberrant proteins released by myeloma cells. Most common signs and symptoms of MM include hypercalcemia, renal function impairment, anemia and lytic bone disease [2]. Major progress in treatment of multiple myeloma has been seen mainly due to development of novel agents, as proteasome inhibitors. Although survival of MM has increased, most of the patients ultimately relapse and become refractory [3]. Therefore, novel therapeutic options are needed. Dysregulated metabolism has been considered as one of the hallmarks of cancer [4]. The "Warburg effect" gave the first insights into the altered metabolism of cancer cells reporting that tumor cells prefer aerobic glycolysis for nutrient production and proliferation [5]. Recently, the field of metabolomics is attracting considerable interest due to novel methods and technical progress, which enable high quality measurements and contributed significantly to a better understanding of metabolic rewiring in pathological diseases as cancer [4]. Metabolic profiling in cancer emerges as a tool for the diagnosis, classification, treatment decisions and assessment of treatment efficacy, and identification of novel therapeutic targets. Metabolic profiling can be performed in different tissue samples and biofluids, such as serum, plasma, saliva, urine and in a high throughput fashion [6]. Nevertheless, metabolomics reporting and data analysis is not yet standardized and often more technical oriented and less biology driven [7].
The first metabolic analysis of MM cells identified that MM depend on glucose and glutamine metabolism [8]. Higher levels of isoleucine, arginine, acetate, phenylalanine, and tyrosine, and decreased levels of 3-hydroxybutyrate, lysine, glutamine, and some lipids were observed in myeloma patients at diagnosis, but not after achieving remission [9]. Bajpai et al. showed that targeting glutamine metabolism sensitizes MM cells to the bcl-2 inhibitor venetoclax. Cellular metabolic investigations revealed LDHA and HIF1α as novel targets for drug resistance in MM under bone marrow hypoxic conditions. Inhibition of LDHA and HIF1A can restore sensitivity to therapeutic drugs such as bortezomib [10]. Ludwig et al. highlights alterations in bone marrow metabolism as an early feature of the development of MGUS and MM [11]. Despite the interest in metabolomics of MM, the role and potential application in diagnostics, classification and prediction of therapy response remains unclear. Therefore, this study focuses to determine differences in the metabolic phenotype between healthy subjects, MGUS-, NDMM-and RRMM patients.

Ethics statement
Investigations have been conducted in accordance with the Declaration of Helsinki and the national and international guidelines and have been approved by the authors' institutional review board (Ethics committee of the Medical University of Innsbruck, Austria, number: 1085/2017 and Ethics committee of the Faculty of Medicine, Masaryk University, Brno, Czech Republic, number: 20/1/2011). Patients provided written informed consent.

Patients and sample collection
Patients with monoclonal gammopathy of undetermined significance (MGUS, n = 15), newly diagnosed MM (NDMM, n = 32) and relapsed/refractory MM (RRMM, n = 19) according to the International Myeloma Working Group (IMWG) criteria [12], were included in the study. All relevant clinical data and disease characteristics are shown in Table 1. Peripheral blood samples underwent centrifugation for 10 min at 1000-x g and obtained peripheral blood plasma was collected and stored at -80˚C.

Metabolite analysis
A targeted metabolomics approach with the AbsoluteIDQ™ p180 kit (BIOCRATES Life Sciences AG, Innsbruck, Austria) was used for quantification, based on electrospray ionization liquid (ESI-LC-MS/MS) and flow-injection analysis mass spectrometry (FIA/MS) measurements. The assay allows simultaneous quantification of in total 188 metabolites out of 10 μL peripheral blood plasma, including amino acids, biogenic amines, acylcarnitines, glycerophospholipids, sphingomyelins, and the sum of hexoses. Analytic measurements were carried out on an API 4000 and API 5500 LC-MS/MS System (AB Sciex Deutschland GmbH, Darmstadt, Germany) controlled by the Analyst 1.5.1 and Biocrates MetIDQ software. For the calculation of metabolite concentrations, external standards served as a reference. Concentrations of all metabolites were calculated in μM and normalized with respect to internal quality control samples.

Metabolomic data pre-processing and statistical analysis
All multivariate (Principal Component Analyses, PCA and Partial-Least-Squares Discrimination Analysis, PLS-DA and Hierarchical Cluster Analysis, HCA) and univariate statistical analyses were performed with the statistic platform R (Version 3.2.4). The raw data was cleaned applying a modified 80% rule (at least 80% valid values above the limit of detection (LOD) need to be available per analyte in the samples for each group), which resulted in exclusion of 48 metabolites (S1 Table). Remaining values below LOD in the dataset were replaced applying a logspline imputation method with values between LOD and LOD/2 (Parts A and B in S2 Table) [13]. To meet assumptions of statistical tests, data were additionally log2 transformed. Analysis of variance (ANOVA) and post-hoc t-tests were performed to identify significant metabolite alterations between the different patient groups and healthy controls. Data were corrected for multiple testing using the Benjamini-Hochberg (BH) procedure. Metabolites with p-values below a significance level of α = 0.05 were considered as statistically significant.

Supervised metabolic analysis of MM patients versus healthy controls
Multivariate supervised PLS-DA, that relies on the class membership of each observation, was used on the transformed metabolomic data set (in total 188 metabolites) to identify overall differences. First PLS-DA was applied to compare peripheral blood plasma samples of healthy controls with those of the MM group (MGUS+NDMM+RRMM). As shown in the score-plot (Fig 1), the samples from controls were largely separated from the MM group. This finding indicates clearly a distinct metabolic composition of these two populations. Second, PLS-DA was used to detect separations between MGUS-, NDMM-and RRMM patients (Fig 2). Because of the high complexity and in-group variation in the data, differences based on the largest principle components PC1 and PC2, were not as prominent as compared to controls. The substantial overlap remaining between the three patient groups indicates an overall similar metabolic composition. The difficulty in separating these three individual groups by the measured metabolites is likely due to the fact that their differences that segregate the samples are reduced by the intrinsic variations found from patient to patient.
In addition, based on the known percentage of plasma cells in the bone marrow of MM patients, we evaluated differences in the metabolic profile with respect to this parameter. In total, 18 metabolites from various classes such as the branched-chain amino acids (BCAAs; isoleucine, leucine and valine), various medium-and long-chain acylcarnitines (C10, C16 and C18) were significantly altered between the plasma cell ranges (S3 Table).

Metabolite alterations among the healthy controls and MGUS/MM patients
Student's t-tests, fold change calculations and biochemical pathway mapping were performed in the following data analysis step to gain more precise insights about metabolic alterations between the four different groups. Prominently, three branched-chain amino acids (leucine, isoleucine, and valine) displayed low concentrations in NDMM patients, indicating a possible higher consumption of these metabolites (uptake by the MM cells in the bone marrow). In addition, glutamate, one of the most abundant amino acids in the human body, was significantly increased in NDMM patients compared to healthy controls, suggesting that MM cells in the bone marrow secrete or detoxify glutamate. In line, a recent meta-analysis of clinical metabolic profiling in blood and tissue of cancer patients identified glutamate as the most increased metabolite in blood as well as in tissue [7].
The activation of the kynurenine pathway occurs when the activity of indoleamine 2,3 dioxygenase-1 (IDO1), which catalyzes the conversion of tryptophan to kynurenine, is increased [14,15]. This process reduces blood tryptophan levels, which in turn decrease the availability of a higher kynurenine/tryptophan ratio, as detected in the NDMM patients. In addition, eight lysophosphatidylcholines (LysoPCs) were significantly lower in NDMM patients. Out of the 25 significantly altered phosphatidylcholines (PCs), 23 were increased in NDMM samples. Metabolite differences in acylcarnitines, amino acids and biogenic amines between healthy controls and RRMM patients were comparable to those observed between controls and NDMM patients. 32 metabolites were significantly altered and long-chain acylcarnitines (C16, C18, C18:1) were elevated in RRMM patients in comparison to the control group ( Fig  3C and Part C in S1 Fig). All three BCAAs and ADMA were increased in healthy subjects. The kynurenine/tryptophan ratio, reflecting possible increased IDO activity, was higher in RRMM patients.
Metabolic alterations between healthy controls and MGUS patients (in total 36 significant metabolites) were similar to the changes in the RRMM and NDMM groups, with significantly enriched plasma acylcarnitines in MGUS samples (C2, C18, C18:1 and C18:2). Increased kynurenine pathway activation and lower ADMA levels could be again detected in MGUS patients.
Differences between NDMM-, RRMM-and MGUS patients were not as prominent as compared to the healthy control group with fewer significantly affected metabolites. In total, eight metabolites were significantly altered between the NDMM and the MGUS group, including free carnitine, acetylcarnitine, glutamate, ADMA, and four PCs species (Fig 3D and Part D in S1 Fig). Eight metabolites (ADMA, octadecanoylcarnitine and 6 PCs) were significantly altered between RRMM-and MGUS patients (Fig 3E and Part E in S1 Fig). In RRMM-and NDMM patients the 13 significantly altered metabolites comprised free carnitine, acetylcarnitine, creatinine, several LysoPCs and PCs (Fig 3F and Part F in S1 Fig). Mean concentration and standard deviation (SD) of significantly altered metabolites are shown in Parts A-F in S4 Table. Discussion MM is a heterogenic disease with a dynamic metabolic and proteomic process in bone marrow and its microenvironment [16]. In the literature only marginally reports exist about a specific metabolic profile in MM [10]. In this study we determined differences in metabolic phenotypes between healthy subjects, MGUS-, NDMM-and RRMM patients.
Using multivariate PLS-DA of peripheral blood plasma samples, we could clearly separate healthy controls from MGUS-, NDMM-, and RRMM patients. Thus, it appears that MM patients have a different metabolic profile than healthy subjects. We found 57 metabolites to be significantly different between healthy controls and NDMM patients. Free carnitine (C0), acetylcarnitine (C2) and five long-chain acylcarnitines (C14:1, C16, C18, C18:1 and C18:2) were all increased in the NDMM population. Carnitine is essential for the transport of longchain acyl groups from fatty acids into the mitochondrial matrix. Inside the mitochondrial matrix fatty acids can be broken down to acetyl-CoA through a process called ß-oxidation, so they can enter the tricarboxylic acid (TCA) cycle for energy production [17]. Thus, acylcarnitines represent the carrier form of activated fatty acids for the transport across the inner mitochondrial membrane and alterations potentially point toward changes in lipid breakdown.
Moreover, 36 amino acids and biogenic amines were significantly increased or decreased between controls and MGUS patients. The three BCAAs (leucine, isoleucine, valine) had a lower concentration in NDMM patients. Notably, BCAAs are not only an essential nutrient source, but can also function as potent signaling molecules and changes in their levels have been linked to cancer progression [18]. BCAAs can be efficiently used for protein synthesis or oxidized for energy purposes by cancer cells. As BCAAs are essential, cancer fully rely on dietary BCAA intake or their release from protein degradation [19]. In addition, it has been shown that the enzymes catalyzing the first step in BCAA degradation are overexpressed in many cancer types [20,21]. The cytosolic branched chain-amino transferase 1 (BCAT1), which converts BCAAs to their corresponding branched-chain a-keto acids is involved in cancer proliferation and has been proposed as a prognostic cancer marker [18,[22][23][24]. Moreover, a study investigating the metabolic profile of bortezomib resistance in cell culture and primary cells showed that the serine synthesis pathway (SSP) has significantly increased activity in bortezomib resistant MM and plasmacytoid lymphoma cell lines. Importantly, the study also observed a strong correlation between SSP activity and the ability of cells to withstand increasing bortezomib concentrations in all bortezomib resistant cell lines tested [25]. In this concern, it may be that serine and glycine were significantly increased in our NDMM-and RRMM patients.
Furthermore, kynurenine, which is converted from the essential aromatic amino acid tryptophan by the enzyme IDO1 plays an important role in several cellular functions. The activity of IDO and the activation of the kynurenine pathway, which is indirectly reflected by the kynurenine to tryptophan ratio, was elevated in NDMM patients. The IDO enzyme has been reported to be a key player in both cancer development and progression as it supports inflammatory processes in the tumor microenvironment, mediates immune tolerance to tumor antigens and suppresses T-cells and natural killer cells [26,27]. In addition, the extent of the abnormality in tryptophan metabolism has been shown to directly correlate with the aggressiveness of the malignancy [28,29]. Therefore, the results of this study corroborate that the aberrant kynurenine pathway activation by IDO may serve as a potential cellular mechanism, which is present in NDMM patients to induce a systemic deregulation of immune responses [30,31]. Based on the evidence for immune tolerance, there has been increasing interest in IDO as a novel therapeutic target for the development of new anti-cancer and anti-myeloma drugs. Thus, promising IDO-inhibiting drugs for use in multiple myeloma are now the focus of research.
Glutamate, which was elevated in NDMM patients, is a nonessential amino acid and serves as the major bioenergetic substrate for proliferating normal but also neoplastic cells. Glutamate functions as an excitatory neurotransmitter that plays essential roles in metabolic, and oncogenic signaling pathways [32]. As an intermediate in glutaminolysis, glutamine is converted into TCA cycle metabolites through the activity of multiple enzymes. Aberrant glutaminolysis has been described as a hallmark of cancer, as many cancer cells undergo metabolic reprogramming that makes them highly glutamine dependent for their survival and proliferation [33,34]. Especially, targeting the glutaminolysis pathway by inhibition of glutamine synthase (GS), cellular glutamine transporter (GLUTs) or oncogenes involved in the regulation of this metabolic pathway is a promising approach for heavily pretreated multiple myeloma patients [35].
Asymmetric dimethylarginine (ADMA), which was enriched in healthy controls, is an endogenous inhibitor of nitric oxide synthase (NOS) and is derived from methylation of arginine residues in proteins. Increased blood levels of ADMA have been linked to cardiovascular diseases, renal failure and hypertension, in which ADMA was described as an independent risk factor [36,37].
Moreover, in our study we observed that eight lysophosphatidylcholines (LysoPCs) were significantly lower in NDMM patients. LysoPCs result from the partial hydrolysis of PCs, at which one fatty acid group is enzymatically removed by Phospholipase A2 [38]. LysoPCs are emerging as a novel class of inflammatory lipids, joining thromboxanes, leukotrienes and prostaglandins, with which they share common metabolic pathways and regulatory mechanisms.
Additionally, we found that 32 metabolites were significantly changed and again several long-chain acylcarnitines (C16, C18, C18:1) were elevated in the RRMM group, indicating again changes in overall lipid energy metabolism. In line with the NDMM group, all three BCAAs as well as ADMA were increased in the control samples, whereas a higher kynurenine/ tryptophan ratio was present in RRMM patients, pointing towards an elevated kynurenine pathway activation. Compared to the NDMM group, a smaller number lysoPCs and PCs were significantly altered between controls and RRMM patients.
The metabolic alterations between healthy controls and MGUS patients were again similar to the changes in the RRMM-and NDMM groups, with significantly enriched peripheral blood plasma acylcarnitines in MGUS samples (C2, C18, C18:1 and C18:2). An increased kynurenine/tryptophan ratio and lower ADMA levels again could be detected in MGUS patients. Metabolic alterations of several different biochemical classes in the bone marrow environment from MGUS patients were recently reported [11]. These changes, especially for the lipid metabolism, could also be observed in our study with peripheral blood.
In our study we found significantly different changed metabolites in MGUS, NDMM and RRMM patients. Moreover, we identified a distinct metabolic profile with significant alterations in amino acid, lipid and energy metabolism in healthy volunteers compared to MGUS/MM patients. We assume that several cellular metabolic processes, most likely immunoregulation by IDO or glutaminolysis may serve as novel promising therapeutic targets in MM. Promising future treatments may combine approaches to target different metabolic pathways such as glutamine metabolism, glycolysis or the TCA cycle. Further investigations are needed to examine the association of certain metabolites with disease progression and drug resistance.
Supporting information S1


extfieldsofstudy: Chemistry Medicine

text:
In Vitro-Generated Tc17 Cells Present a Memory Phenotype and Serve As a Reservoir of Tc1 Cells In Vivo

Memory CD8+ T cells are ideal candidates for cancer immunotherapy because they can mediate long-term protection against tumors. However, the therapeutic potential of different in vitro-generated CD8+ T cell effector subsets to persist and become memory cells has not been fully characterized. Type 1 CD8+ T (Tc1) cells produce interferon-γ and are endowed with high cytotoxic capacity, whereas IL-17-producing CD8+ T (Tc17) cells are less cytotoxic but display enhanced self-renewal capacity. We sought to evaluate the functional properties of in vitro-generated Tc17 cells and elucidate their potential to become long lasting memory cells. Our results show that in vitro-generated Tc17 cells display a greater in vivo persistence and expansion in response to secondary antigen stimulation compared to Tc1 cells. When transferred into recipient mice, Tc17 cells persist in secondary lymphoid organs, present a recirculation behavior consistent with central memory T cells, and can shift to a Tc1 phenotype. Accordingly, Tc17 cells are endowed with a higher mitochondrial spare respiratory capacity than Tc1 cells and express higher levels of memory-related molecules than Tc1 cells. Together, these results demonstrate that in vitro-generated Tc17 cells acquire a central memory program and provide a lasting reservoir of Tc1 cells in vivo, thus supporting the use of Tc17 lymphocytes in the design of novel and more effective therapies.

Memory CD8 + T cells are ideal candidates for cancer immunotherapy because they can mediate long-term protection against tumors. However, the therapeutic potential of different in vitro-generated CD8 + T cell effector subsets to persist and become memory cells has not been fully characterized. Type 1 CD8 + T (Tc1) cells produce interferon-γ and are endowed with high cytotoxic capacity, whereas IL-17-producing CD8 + T (Tc17) cells are less cytotoxic but display enhanced self-renewal capacity. We sought to evaluate the functional properties of in vitro-generated Tc17 cells and elucidate their potential to become long lasting memory cells. Our results show that in vitro-generated Tc17 cells display a greater in vivo persistence and expansion in response to secondary antigen stimulation compared to Tc1 cells. When transferred into recipient mice, Tc17 cells persist in secondary lymphoid organs, present a recirculation behavior consistent with central memory T cells, and can shift to a Tc1 phenotype. Accordingly, Tc17 cells are endowed with a higher mitochondrial spare respiratory capacity than Tc1 cells and express higher levels of memory-related molecules than Tc1 cells. Together, these results demonstrate that in vitro-generated Tc17 cells acquire a central memory program and provide a lasting reservoir of Tc1 cells in vivo, thus supporting the use of Tc17 lymphocytes in the design of novel and more effective therapies.
Keywords: Tc17 cells, cD8 + T cell memory, secondary expansion, persistence, homing, oxidative metabolism inTrODUcTiOn Effector CD8 + T cells can be classified into at least two functional subsets. Type 1 CD8 + T (Tc1) cells produce interferon (IFN)-γ and are endowed with high cytotoxic capacity. Type 17 CD8 + T (Tc17) cells produce IL-17 and display enhanced self-renewal potential. The presence of Tc17 cells has been described in several infectious and autoimmune diseases in both human and mice (1,2), thus confirming the role of Tc17 cells in inflammatory processes. The accumulation of Tc17 lymphocytes has also been observed in several types of human tumors including gastric (3), uterus (4), head and neck (5), liver (6), and colon cancer (7) and also in various cancer models such as B16 melanoma (8). In spite of their relatively low abundance in several tumors, we and others have reported that the adoptive transfer of in vitro-generated Tc17 cells promotes the elimination of established tumors in C57BL/6 mice, supporting evidence of antitumor properties of Tc17 lymphocytes (9)(10)(11)(12).
The fact that in vitro-generated Tc17 lymphocytes do not express cytotoxic molecules such as granzymes and perforins (13) and thus are not able to induce apoptosis of target cells (14) opens questions regarding the mechanisms that could explain the antitumor potential of Tc17 lymphocytes in vivo. In search of evidence which may account for this dichotomy, we demonstrated that Tc17 cells are highly plastic and some cells differentiate into Tc1 IFN-γ + /granzyme B + cells upon transfer into tumor-bearing mice (12). Moreover, Restifo and collaborators further suggested that the effectiveness of Tc17 lymphocyte adoptive therapy depends on the ability of these cells to differentiate into IFN-γ-producing cells (9). Accordingly, Dutton's group reported that after adoptive transfer, Tc17 cells promote the recruitment of macrophages and neutrophils to the tumor microenvironment through the production of IFN-γ, TNF-α, and IL-17 (10).
Recently, it has been suggested that Tc17 cells may display phenotypic features of memory T cells. The group of Takiguchi described that human Tc17 cells are found within the CD28 + CD45RA − memory subset in healthy donors (15). Also, other reports have shown that both endogenous and in vitro-generated Tc17 cells express high levels of the TCF-7 transcription factor (16,17), suggesting an active Wnt/β-catenin signaling, which constitutes a major pathway in maintaining the molecular program of memory T lymphocytes. In the context of cancer, it has been reported that adoptively transferred Tc17 lymphocytes have greater persistence in vivo compared to Tc1 lymphocytes (9) and that Tc17 lymphocytes exhibit high IL-7 receptor alpha chain (CD127) expression (11), which may explain their increased survival capacity. Finally, it has been shown that Tc17 cells can differentiate into tissue-resident memory T cells after the resolution of psoriasis (18). However, it is unclear whether Tc17 lymphocytes possess other functional features of memory cells, such as the ability to mount a secondary response and an OXPHOS-based metabolism.
Based on these previous findings, we sought to evaluate whether Tc17 cells present functional properties of memory CD8 + T cells. For this, we generated OVA-specific Tc17 cells and transferred the cells into wild-type mice to test their capacity to expand upon a secondary challenge. We show that Tc17 cells persist in the long term, expand upon a secondary stimulation with OVA, shift to a Tc1 phenotype and recirculate within secondary lymphoid organs. In agreement, Tc17 cells are endowed with a higher mitochondrial spare respiratory capacity (SRC) compared to conventional (Tc1) cytotoxic CD8 + T cells. Finally, Tc17 cells express higher levels of memory-related markers and lower levels of molecules required for the effector program compared to Tc1 cells. Together, these results suggest that Tc17 cells present functional properties of central memory CD8 + T cells.

MaTerials anD MeThODs
Mice C57BL/6 (CD45.2 + ), B6SJL-PTPRC (CD45.1 + ), OT-I, and Rag1 −/− mice were purchased from The Jackson Laboratory. All mice were kept in the animal facility at Fundación Ciencia y Vida and maintained according to the "Guide to Care and Use of Experimental Animals, Canadian Council on Animal Care. " Animal work was carried out under institutional regulations of Fundacion Ciencia & Vida and Facultad de Ciencias, Universidad de Chile and was approved by the local ethics review committees.

intracellular staining and Flow cytometry
To assess cytokine production by adoptively transferred Tc1 and Tc17 cells, cells obtained from lymph nodes, spleen, and peritoneal cavity of recepient mice were re-stimulated with 0.25 µM PMA (Sigma-Aldrich) and 1 µg/ml ionomycin (Sigma-Aldrich) in the presence of GolgiPlug (BD Biosciences) for 4 h. Cells were stained with antibodies against the cell surface markers (congenic markers) and then resuspended in a fixation/ permeabilization solution (Cytofix/Cytoperm; BD Pharmingen). Following fixation and permeabilization, the cells were incubated with antibodies against IFN-γ, IL-2, IL-17, and TNF-α for 30 min at 4°C. The cells were then washed with permeabilization buffer and resuspended in PBS supplemented with 2% FCS for FACS analysis (FACSCanto II; BD Bioscience). In some cases, Fixable Viability Dye (eBioscience) was used to discard dead cells from the analysis.
Intracellular staining for Ki67 was conducted by incubating surface stained cells in a fixation/permeabilization solution (Foxp3/transcription factor permeabilization and fixation buffer, eBioscience). Following fixation and permeabilization, the cells were incubated with the Ki-67 antibody (clone 11F6, Biolegend) in a permeabilization buffer (Foxp3/transcription factor permeabilization buffer, eBioscience). After staining, cells were washed and analyzed by FACS. FACS data analysis was performed using the FLOWJO software (Tree Star Inc., Ashland, OR, USA). qPcr Tc1 and Tc17 cells were differentiated in vitro, and total RNA was extracted. Naive CD8 + T cells were isolated by cell sorting (CD8 + /CD44lo/CD62Lhi/CD25 − ) before RNA extraction. Total RNA was obtained using EZNA Total RNA Kit I (Ω Bio-Tek). RNA (1 µg) was reverse-transcribed using M-MLV reverse transcriptase (Promega). The PCR reaction was performed using Brilliant II SYBR Green QPCR Master Mix (Agilent Technologies) in a Stratagene Mx3000P real-time PCR machine. For relative quantitation, the amplified fragments were normalized according to constitutive transcription of the housekeeping gene HPRT. The sequences of the primers used for quantification of each measured transcript were the following:

Western Blot
Cell lysis was performed using Cytobuster protein extraction reagent (Novagen) supplemented with protease and phosphatase inhibitors (complete PhosSTOP, Roche). Protein extracts were separated in 15% SDS-polyacrylamide gels and transferred to PDVF membranes (Millipore). Blocking was at room temperature for 1 h in 5% fat-free milk, and membranes were incubated overnight at 4°C with 1:2,000 dilution of Total OXPHOS Human WB Antibody Cocktail (ab110411), and then for 1 h at room temperature with a secondary antibody conjugated to horseradish peroxidase. Chemiluminescence detection used SuperSignal West Pico Chemiluminescent Substrate (Thermo Scientific) and series of timed exposure images were acquired with a Chemidocs imaging system (BioRad) to ensure densitometric analyses were performed at exposures within the linear range. To ensure equal protein loading across gels, membranes were submerged in stripping buffer (Restore Western blot stripping buffer; Pierce), incubated at 37°C for 30 min, and re-probed with a loading control antibody. Image J was used for densitometric analysis.
adoptive Transfer of Tc1 and Tc17 cells To assess secondary expansion of transferred cells, Tc1 and Tc17 lymphocytes were generated in vitro from OT-I/CD45.2 + and OT-I CD45.1 + /CD45.2 + mice, respectively. On the day of the transfer, the Tc1 and Tc17 lymphocytes were counted and mixed at a 1:1 ratio. The initial ratio (input) on the cytometer was then evaluated by staining against CD45.1 and CD45.2. A total of 2 × 10 6 cells from this mixture were transferred intravenously into B6.SJL (CD45.1 + ) congenic mice. At day 23 post-adoptive transfer, the mice were immunized with 500 µg of OVA protein (Sigma-Aldrich) and 100 µl of Freund's Complete Adjuvant (CFA, Sigma) intraperitoneally. On days 10, 20, and 27 post-adoptive transfer blood samples were obtained and analyzed by flow cytometry. Finally, on day 27 after the adoptive transfer, the mice were sacrificed, and cells obtained from the peritoneal cavity, secondary lymphoid organs, and spleen were analyzed by FACS.
To examine the ability of Tc17 cells to be retained in effector tissues, we performed intradermal injections of in vitrodifferentiated CD45.2 + Tc1 and CD45.1 + /CD45.2 + Tc17 cells (1.5 × 10 6 cells per mice, at a 1:1 ratio) into CD45.1 + congenic mice, and 30 days later, we analyzed skin, secondary lymphoid organs, spleen, bone marrow, lung, and liver for the transferred cells. The intradermal transfer was performed by injecting 100 µl of cellular suspension in PBS divided in five spots of 20 µl each in the lower back skin of each mouse.

homeostatic survival
For in vitro homeostatic survival analysis, Tc1 and Tc17 cells were cultured in the presence of IL-15 (25 ng/ml, eBioscience) or IL-7 (10 ng/ml, eBioscience). Three days later, cells were harvested, counted, and analyzed by FACS.

seahorse
Oxygen consumption rates (OCRs) and extracellular acidification rates (ECARs) from in vitro generated Tc1 and Tc17 cells (0.5 × 10 6 cells/well) were measured in non-buffered DMEM without phenol red (containing 25 mM glucose, 2 mM lglutamine, and 1 mM sodium pyruvate) under basal conditions and in response to 1 µM oligomycin, 0.5 µM fluoro-carbonyl cyanide phenylhydrazone, and 1 µM rotenone + 1 µM antimycin A (Sigma). The cells were analyzed with the XFe-96 Extracellular Flux Analyzer (Seahorse Bioscience).  Data are presented as mean ± SEM. Differences between two groups were determined using two-tailed Mann-Whitney test. Differences between more than two groups were determined using one-way ANOVA with Tukey or Bonferroni post-test. A two-way ANOVA with Bonferroni post-test was used to compare secondary expansion rates. Statistical analysis and graphs were obtained with GRAPHPAD PRISM (GraphPad Software Inc).

Tc17 cells expand rapidly following a secondary challenge
In a previous study, we have demonstrated that Tc17 cells are endowed with high antitumor capacity and persistence, despite displaying a reduced cytotoxic potential in vitro (12). These observations lead us to hypothesize that Tc17 cells present a memory phenotype. To test this hypothesis, we co-transferred in vitrogenerated Tc1 and Tc17 cells obtained from OT-I mice into wild type mice and then 3 weeks later, the mice were immunized i.p. with OVA-CFA. We analyzed the persistence and expansion of Tc1 and Tc17 cells before and after OVA immunization. As shown in Figure 1A, 10 and 20 days after the co-transfer, we could only detect the presence of Tc17 cells in blood, and these cells suffered a two-fold expansion following immunization with OVA, whereas Tc1 cells remained undetected. Accordingly, Tc17 cells were mainly found in the peritoneal cavity while Tc1 cells were present in low frequencies, ruling out the possibility that Tc1 cells were absent in blood because they were more efficiently recruited to the peritoneal cavity (Figures 1B,C). Interestingly, besides being detected in the peritoneal cavity, Tc17 cells were also present in spleen and secondary lymph nodes such as peripheral lymph nodes (PLN), mesenteric lymph node (MLN), and parathymic lymph nodes (pTLN) (Figures 1B,C). This observation suggests that Tc17 cells may be recruited to effector tissues but also, that a portion of Tc17 cells may remain circulating within secondary lymph nodes as a reservoir for future encounter with their cognate antigen. Next, we tested whether Tc17 cells continued to produce IL-17 or could give rise to IFN-γ-producing cells. As shown in Figure 1D, a high percentage of transferred Tc17 cells differentiated into IFN-γ-producing cells in spleen and peritoneal cavity, whereas a small fraction of Tc17 cells continued to produce IL-17 in spleen, confirming their high plasticity in this setting. Tc17 cells also expanded more rapidly and outnumbered Tc1 cells following in vitro polyclonal activation with anti-CD3 and anti-CD28 antibodies (data not shown). Taken together, these results demonstrate the ability of Tc17 cells to persist in the long term and rapidly expand upon a challenge. Next, we evaluated whether Tc17 cells could differentiate to memory cells in vivo which could expand rapidly following a rechallenge with their cognate antigen. For this, we co-transferred in vitro-generated Tc1 and Tc17 cells obtained from OT-I mice into wild-type mice, and 24 h later, the mice were challenged with OVA. Four weeks later, mice were re-challenged with OVA and the frequencies of transferred cells in blood, PLN and spleen were analyzed. In this setting, the composition of Tc1 and Tc17 cells remained similar to the input before the re-challenge with OVA (dotted line, Figure S1C in Supplementary Material). However, following re-challenge with OVA protein (day 34), we found that Tc17 cells rapidly expanded ( Figures S1A,B in Supplementary Material) and became almost 90% of transferred cells ( Figure  S1C in Supplementary Material), confirming that Tc17 cells are superior than Tc1 cells in proliferating in response to their cognate antigen. Similar results were obtained when analyzing the composition of transferred cells in spleen and PLN ( Figure  S1C in Supplementary Material). Interestingly, 1 week after the re-challenge, Tc17 cells found in the spleen and PLN gave rise to IFN-γ-producing cells (54.6 ± 7.3%) that do not express KLRG1 but express high levels of CD44 and CD127 ( Figure S1D in Supplementary Material). Moreover, we found that a proportion of these Tc17 cells (29.6 ± 3.7%) express CD62L, which explains their potential to migrate within secondary lymph nodes ( Figure  S1D in Supplementary Material).

Tc17 cells Present a Migration Pattern of central Memory T cells
In agreement with our previous observations in tumor bearing mice (12), we demonstrated that transferred Tc17 cells were always found in the spleen and secondary lymphoid organs (Figure 1B; Figure S1C in Supplementary Material), suggesting that they have a circulation pattern similar to central memory (TCM) cells. To study the migration pattern of Tc17 cells, in vitro-generated Tc17 and Tc1 cells were stained and transferred into congenic mice. 24 h later, we analyzed spleen, PLN, MLN, bone marrow, liver, lung, and blood for the presence of the transferred cells. As shown in Figure 2A, Tc17 cells migrated preferentially to secondary lymphoid organs such as the spleen, PLN, and MLN compared to Tc1 cells. In contrast, Tc1 cells preferentially migrated to the bone-marrow compared to Tc17 cells (Figure 2A).
Next, we tested whether Tc17 cells injected in peripheral tissues can become established as memory cells in secondary lymphoid organs. For this, we intradermally co-transferred Tc1 and Tc17 cells in a 1:1 ratio into CD45.1 + mice and 30 days later, the presence of memory cells was analyzed in skin, spleen, PLN, MLN, bone marrow, liver, lung, and blood. As shown in Figure 2B, we did not find Tc1 or Tc17 cells in the skin, however, Tc17 cells were found in blood and in all the organs analyzed. Interestingly, when analyzing the absolute numbers, we observed that Tc17 cells were preferentially found in secondary lymphoid organs ( Figure 2B). These results demonstrate that Tc17 cells are not retained in the skin but migrate into secondary lymphoid organs.

Tc17 cells Proliferate More robustly than Tc1 cells in response to homeostatic stimuli
To analyze the homeostatic response of Tc17 cells, we cultured Tc17 and Tc1 cells with the homeostatic cytokines IL-15 and IL-7, and the proliferation and absolute numbers of cells were analyzed 3 days later. As shown in Figure 3A, Tc17 cells responded to both homeostatic cytokines while Tc1 cells proliferated only in response to IL-15, which is in agreement with the fact that only Tc17 cells express CD127. In spite of their rapid proliferation, Tc17 cells retained the ability to produce IL-17 and only a small fraction of cells produced IFN-γ (data not shown).
To further analyze the in vivo homeostatic proliferation of Tc17 cells, Tc1 and Tc17 cells were transferred into Rag1 −/− mice at a 1:1 ratio and several days later, the absolute number of transferred cells was analyzed. As shown in Figures 3B-D, both Tc1 and Tc17 cells were present in Rag1 −/− mice 11 and 25 days following the transfer, but the absolute number of Tc17 cells in blood always doubled the absolute number of Tc1 cells. Moreover, a higher percentage of Tc17 cells were positive for Ki67 staining at day 25 compared to Tc1 cells ( Figure 3E). Thus, Tc17 cells are much more responsive to homeostatic stimulation than Tc1 cells.

Tc17 cells Present a Transcriptional and Metabolic Profile of Memory T cells
We have previously shown that Tc17 cells express higher levels of tcf-7 and lef1 and lower levels of GzmB compared to Tc1 cells (12). To further characterize the phenotype of Tc17 cells, in vitrogenerated Tc1 and Tc17 cells were analyzed by real-time PCR and by FACS in search of additional markers related to memory cells and effector cells. As shown in Figure 4A, we observed that Tc17 cells express markers related to memory T cells such as bcl-6 and lack markers of terminal differentiation such as tbx-21, eomes, and prdm1. On the other hand, Tc1 and Tc17 cells did not express Klf-2 transcriptional factor, which commands CD62L expression. Interestingly, in contrast to Tc1 cells, Tc17 cells present low levels of CD25, granzyme B, and TNF-α, while they express RORγt and produce high levels of IL-2 ( Figure 4B).
Since Tc17 cells express high levels of bcl-6 and low levels of prdm1 which suggests that Tc17 cells present a metabolic signature related to memory T cells, we sought to determine the metabolic profile of Tc17 cells. A distinctive characteristic of memory CD8 + T cells is that they rely on an oxidative metabolism and present a high SRC (19). Thus, we determined the OCR and ECAR of in vitro-generated Tc1 and Tc17 cells and found that Tc17 cells present a higher OCR and SRC than Tc1 cells (Figures 5A,B). Moreover, Tc17 cells exhibit higher levels of complex IV proteins compared to Tc1 cells (Figures 5C,D and Figure S2 in Supplementary Material), suggesting that they present an increased efficiency in respiration. All these data demonstrate  that Tc17 cells rely mainly on an oxidative metabolism which has been already related to a memory CD8 + T cell phenotype.

DiscUssiOn
Two main features characterize memory T cells; first, they are long-lived and maintained independently of stimulation, and second, they are functionally improved and respond more quickly and efficiently to repeated encounters with a specific pathogen than naive T cells. Also, depending on their migratory preferences, they can be classified into at least three subtypes: central memory T cells (TCM), effector memory T cells (TEM), and tissueresident memory T cells (TRM) (20)(21)(22). Although the potential of naive CD8 + T cells to generate memory cells has been extensively studied, the potential of different CD8 + T cell effector subsets, especially Tc17 cells, to persist and become memory cells has not been fully characterized. Here, we show that Tc17 cells present functional features of TCM cells such as the long-term maintenance, secondary expansion capacity, an oxidative metabolic profile, and a preferential trafficking within peripheral lymphoid tissues. An essential property required for an efficient recall response is the ability of memory T cells to proliferate in response to antigenic stimulation (secondary expansion). In this work, we observed that unlike Tc1 cells, Tc17 lymphocytes were able to rapidly expand in vivo after secondary exposure to the antigen. This quality has been linked to the production and autocrine signaling of IL-2 (23), and consequently, it has been reported that T cells with a lower degree of differentiation such as TCM produce higher levels of IL-2 in comparison to more differentiated T cells as TEM (20,24,25). Moreover, it has been suggested that IL-2 production by TCM cells is important in their enhanced proliferative capacity (26). Consistently, we observed that a high percentage of Tc17 lymphocytes generated in vitro produce IL-2, and even a fraction of the secondary progeny of Tc17 lymphocytes continues to produce IL-2 in vivo (12).
Recent studies have shown that CD25 (IL-2 receptor alpha chain) expression may determine the fate between memory versus effector cells (27,28). It has been suggested that during the asymmetric distribution of CD25 during the first cell division of naive T cells, the population of CD25 low cells give rise  to memory T lymphocytes (28). We have shown that Tc17 cells express low levels of CD25, which is consistent with the phenotype of memory precursors as previously reported by Arsenio and collaborators. The ability of memory T lymphocytes to persist in vivo depends primarily on exposure to cytokines such as IL-7 and IL-15 (29), which promote cell survival by inducing the expression of antiapoptotic proteins such as Bcl-2 (30). In this work we observed that Tc17 lymphocytes show a high expression of the IL-7 receptor alpha chain (CD127) in vivo, similar to that reported in the literature (11), suggesting that Tc17 lymphocytes may have advantages in long-term survival. In agreement with this idea, we observed that Tc17 lymphocytes expand more rapidly than Tc1 cells in vitro when cultured in the presence of IL-7 and also when transferred into Rag1-deficient mice. In support of this idea, the group of Klein has recently reported that vaccine-induced memory Tc17 cells express CD127, are durable and stable, and can persist for a year as functional IL-17A + memory cells (31).
Interestingly, only Tc17 lymphocytes were found in lymphoid organs (spleen and lymph nodes) of mice after antigenic stimulation. This pattern of migration through the lymphoid organs is a hallmark of central memory T lymphocytes and depends on L-selectin (CD62L) and the chemokine receptor CCR7 expression (20,21). Paradoxically, we observed that Tc17 lymphocytes generated in vitro do not express CD62L at the time of adoptive transfer (data not shown), nor the Klf-2 transcriptional factor, which commands CD62L expression. However, a fraction of the Tc17 lymphocytes (almost 30%) can re-express CD62L in vivo and maintain its expression even after secondary antigenic stimulation. In agreement, evidence from Arsenio and collaborators has shown that following antigenic stimulation, CD62L is initially down regulated and in time, it is subsequently upregulated in a population of memory T cells (28). This property of migrating through the lymphoid organs could allow the contact of the Tc17 lymphocytes with antigen-presenting cells and in this way favor the secondary response. In this line, previous studies have positively related the ability of the T lymphocytes to migrate to the lymph nodes with greater effectiveness in schemes of antitumor cellular immunotherapy (24,32).
Evidence from the literature suggests that the Wnt-β-catenin pathway promotes the acquisition of the memory T cell phenotype while restricting the differentiation of CD8 + T lymphocytes into terminal effector T cells (32)(33)(34). Interestingly, in this work, we observed that Tc17 lymphocytes generated in vitro show a high expression of the transducers of the Wnt-β-catenin pathway, namely Tcf-7 and Lef-1, similar to the phenotype already reported in memory CD8 + T lymphocytes (32). This is especially relevant since studies conducted by the Restifo group have shown that the pharmacological activation of the Wnt pathway is sufficient to restrict terminal differentiation in CD8 + T lymphocytes, favoring the generation of highly undifferentiated memory T lymphocytes with a potent therapeutic potential (32,35).
Bcl-6 is another important factor in maintaining the memory phenotype and has been reported to accumulate as CD8 + T lymphocytes transit from the effector into the memory phenotype (36). In agreement, overexpression of Bcl-6 in CD8 + T lymphocytes has shown to promote the generation of TCM lymphocytes in vivo (37). The fact that Tc17 lymphocytes express higher levels of Bcl-6, Tcf-7, and Lef-1 could be a consequence of stimulation of the IL-6/STAT-3 pathway. Indeed, several cytokines related to the Tc17 program such as IL-6, IL-21, and IL-23 stimulate STAT-3 phosphorylation, which induces the expression of RORγt and IL-17, and represses the expression of granzyme B and IFN-γ (13,14). Moreover, it has been reported that STAT-3 controls the expression of Bcl-6 (36) and Tcf-7 (38) in CD8 + and CD4 + T lymphocytes, respectively. In fact, STAT3-deficient CD8 + T lymphocytes are short-lived, fail to generate TCM lymphocytes in vivo, and have a lower Bcl-6 expression compared to normal cells (36). Moreover, recent evidence from the group of Bosselut (39) has clearly demonstrated that STAT3 inhibits cytotoxic effector molecules such as granzymes and cytotoxic regulators such as T-bet and Eomes. Thus, signaling through the IL-6/STAT-3 pathway could contribute to the acquisition of the memory program observed in Tc17 lymphocytes.
Besides being involved in inhibiting Blimp-1 expression (40), Bcl-6 participates in the repression of genes that are critical to the glycolysis pathway, favoring the catalytic pathway (41). Memory CD8 + T cells are characterized by being highly dependent on a mitochondrial metabolism (19, 42), and the inhibition of glycolysis has been shown to improve the generation of memory T cells (43). It has been shown that memory CD8 + T cells present a high SRC, high OCRs, and mitochondrial mass compared to naive T cells, which allows explaining their long term survival and function (19, 42). The SRC constitutes the extra mitochondrial capacity available under stress conditions and increased work. In this line, our results demonstrate that Tc17 cells present higher OCR and SRC compared to Tc1 cells, thus suggesting that they exhibit a metabolic signature related to memory CD8 + T cells.
The results presented here demonstrate that following adoptive transfer, IL-17 production by Tc17 cells is lost when these cells enter effector tissues such as the peritoneal cavity; however, we detected a small fraction of Tc17 cells that still produce IL-17 in the spleen. Some reports from the literature demonstrate that IL-17 production by Tc17 and Th17 cells favors the recruitment of neutrophils into different tissues such as the lung in a murine model of lethal influenza infection (44), the tumor in a model of murine melanoma (10) and tumor draining lymph nodes in a murine colon carcinoma model (45). In all these reports, neutrophils have demonstrated to be pivotal in the induction of a protective immunity, thus IL-17 production by Tc17 cells in lymph nodes may be part of the effector mechanisms used by Tc17 cells to promote effective immune responses. Further studies are needed in order to demonstrate a role of IL-17 production by Tc17 cells in the recruitment of neutrophils to secondary lymph nodes.
Both phenotypic and functional properties observed in Tc17 lymphocytes (persistence, secondary expansion, lymph node migration, and self-renewal) strongly suggest that Tc17 lymphocytes possess qualities suitable for use in antitumor adoptive therapies. In agreement with this idea, our group together with other authors has shown that the adoptive therapy with Tc17 lymphocytes generated in vitro promotes the elimination of established melanoma tumors in mice (9, 10, 12) through an IFN-γ-dependent mechanism (10). Interestingly, following adoptive transfer, Tc17 cells gain the ability to produce IFN-γ. In light of these results, we can conclude that the role of Tc17 cells in this context may be to serve as a self-renewing reservoir of cells that will continuously give rise to Tc1-like cells. In spite of this, the therapeutic effectiveness of Tc17 lymphocytes has been questioned by studies showing that Tc1 lymphocytes are more efficient than Tc17 lymphocytes in the eradication of established tumors based on the number of transferred cells (10,12). It is possible to reconcile these results by arguing that the high effector capacity of Tc1 lymphocytes is crucial to mediate a rapid elimination of tumor cells in a first stage, whereas the memory qualities of Tc17 lymphocytes would support a continuous attack in time, conferring immunity for longer periods.

eThics sTaTeMenT
All mice were kept in the animal facility at Fundación Ciencia y Vida and maintained according to the "Guide to Care and Use of Experimental Animals, Canadian Council on Animal Care. " Animal work was carried out under institutional regulations of Fundacion Ciencia & Vida and Facultad de Ciencias, Universidad de Chile and was approved by the local ethics review committees.
aUThOr cOnTriBUTiOns FF-S designed the study, did experiments, analyzed the data, and wrote the manuscript; BC did experiments, analyzed the data, and wrote the manuscript; DF did experiments and analyzed the


extfieldsofstudy: Biology Medicine

s2fieldsofstudy:
Biology
Environmental Science

text:
Contrasting life histories contribute to divergent patterns of genetic diversity and population connectivity in freshwater sculpin fishes

Background Life history characteristics are considered important factors influencing the evolutionary processes of natural populations, including the patterns of population genetic structure of a species. The sister species Cottus hangiongensis and C. koreanus are small bottom-dwelling freshwater sculpin fishes from South Korea that display marked life history divergence but are morphologically nearly indistinguishable. Cottus hangiongensis evolved an ‘amphidromous’ life history with a post-hatching pelagic larval phase. They spawn many small eggs in the low reaches of rivers, and hatched larvae migrate to the sea before returning to grow to maturity in the river mouth. In contrast, C. koreanus evolved a ‘fluvial’ landlocked type with benthic larvae. They release a smaller number of larger eggs, and the larvae undergo direct development, remaining benthic in the upstream rivers throughout their entire lives. We tested whether there were differences in patterns and levels of within-population genetic diversities and spatial population structure between the two closely related Korean sculpins using mitochondrial DNA control region sequences and seven nuclear microsatellite loci. Results The combined analyses of both marker sets revealed that C. hangiongensis harboured considerably higher levels of within-population genetic diversities (e.g. haplotype/allelic richness, heterozygosities) than C. koreanus. In contrast, the fluvial sculpin exhibited noticeably more spatial population structure than did the amphidromous sculpin, as suggested by pairwise FST statistics. The finding that C. hangiongensis individuals comprised a single random mating population across the east-flowing river basins in the Korean Peninsula, whereas C. koreanus individuals comprised genetically discrete individual populations, was further supported by an individual-based Bayesian population assignment and also factorial correspondence analyses. Conclusions The higher genetic diversity, but lower population structure, of the amphidromous sculpin relative to the fluvial sculpin may have resulted from its greater larval dispersal and also possibly, higher fecundity accompanied by an amphidromous life history. Hence, we conclude that contrasting early life histories – including the presence or absence of the pelagic larval phase – may have led to divergent patterns of within-population genetic diversities and spatial population structure between the sister Cottus species following speciation from a common ancestor of marine sculpin. Electronic supplementary material The online version of this article (10.1186/s12862-018-1171-8) contains supplementary material, which is available to authorized users.

Background
Life history traits are considered important biological factors affecting the evolutionary processes, as they shape the patterns of population genetic structure of a species [1][2][3][4][5][6][7]. An organism's life history can be defined as a set of age-or stage-specific traits that contribute to survival and/or reproduction (i.e. fitness), upon which selection acts [8]. Life history attributes can influence population connectivity as they constrain the movement of individuals (i.e. dispersal), which affects the distribution of within-and among-population genetic variation and contributes to the present-day population genetic structure of a species [3][4][5]. Dispersal is, therefore, thought to play a pivotal role in population dynamics, colonization of novel habitats and geographic distribution of freshwater and marine fishes [9]. In this regard, the presence or absence of a post-hatching planktonic larval period at an early life history stage has been suggested to be a key component in characterizing the degree and pattern of population connectivity or geographic population structure in both marine [10,11] and freshwater species [3,5,6,12].
A number of previous studies suggest that life history characteristics are significantly associated with the degree of population connectivity of freshwater fishes [6,13]. In salmonid fishes that have evolved and diverged into migratory and resident populations within species, the greater dispersal capacity of migratory populations has generally resulted in lower degrees of spatial population structure relative to resident populations [14]. Although several studies have been undertaken to understand the impacts of adult movement (i.e. migratory behaviour) on levels of genetic divergence among freshwater fish populations [4,7], the effects of life history divergence with a focus on the presence or absence of the post-hatching larval phase on the genetic structure of freshwater fishes have not been widely tested (but see [3,5]). Therefore, studying evolutionarily closely related species with contrasting life history types over small geographic scales will provide a good opportunity for directly testing the influences of the larval phase on population connectivity. Moreover, it will allow us to determine the role of life history in shaping the patterns of intraspecific genetic variation in light of the evolution of geographic population structure [6].
Freshwater fishes have evolved an extraordinarily diverse array of life history characteristics. The best reported example of the divergent life histories between populations or closely related species of freshwater fishes comes from the salmonid fishes [13]. For instance, Chinook salmon (Oncorhynchus tshawytscha) is characterized by two distinct life history types with different patterns of geographic distribution as 'ocean-type' and 'stream-type'. These life history types are distinguished based on the duration that juvenile salmon remain in freshwater habitats before migrating to the sea [15,16]. Another case of the diverse life history types among closely related freshwater fish species includes the genus Rhinogobius, in which species with or without a pelagic larval phase exist [17,18].
The freshwater sculpins of the genus Cottus are bottom-dwelling fishes that comprise approximately 68 species from the subarctic areas to the temperate northern hemisphere, including Europe, Siberia, Central and East Asia and North America [19][20][21][22]. These typically cold-water fishes are believed to have originated several times in parallel from different ancestral species of marine sculpins [20]. Interestingly, within the genus Cottus, four distinct life history styles have evolved, including marine, catadromous, amphidromous, and fluvial (landlocked/lacustrine) types [20,23].
The sister species C. hangiongensis and C. koreanus [23] have contrasting life histories and differ in reproductive behaviour and geographic distribution, but are morphologically virtually indistinguishable as adults [24]. Cottus hangiongensis shows an amphidromous life history, where they spawn large numbers of small eggs (752-1376 per clutch) in the lower courses of rivers during early spring [24]. The hatched pelagic larvae then migrate to the sea by river flows. After the pelagic larvae spend approximately 4-6 weeks in marine environments, the juveniles return to rivers to grow and then maintain a benthic life in downstream rivers [25,26]. Whether the juveniles of C. hangiongensis return to the same rivers where they hatched, however, remains unknown. On the other hand, C. koreanus shows a fluvial life history, where they produce relatively fewer, but larger eggs (538-880 per clutch), and the larvae undergo direct development, remaining benthic in the upstream rivers immediately after hatching [24,27]. For spawning, although mature individuals of C. hangiongensis migrate downstream to the spawning ground [27], C. koreanus individuals do not move and stay at their natal habitats [28]. These reproductive behaviours, however, imply that both species remain in the same rivers throughout their entire lives as adults. Males of these species usually mate continuously with several females and remain in the nest after spawning (i.e. polygyny). The migration ranges of settled adult individuals are estimated to be around 10 m for C. koreanus [28] and about 100 m for C. hangiongensis [29], which suggests restricted dispersal capacities for both species at an adult life history phase. Cottus hangiongensis is distributed across the eastern Korean Peninsula, northern Japan and Russia, whereas C. koreanus is endemic to the Korean Peninsula [30] (Fig. 1a). Within South Korea, C. hangiongensis occurs in rivers that flow into the East Sea from central to northern parts of eastern Korea, and C. koreanus is distributed in the Han River, Imjin River, and Nakdong River in inland regions [31] (Fig. 1b). In these Korean sculpins, speciation from a common ancestor of marine sculpin was accompanied by ecological divergence into amphidromous and fluvial life histories, as seen for another group of Japanese sculpins, C. amblystomopsis and C. nozawae [32,33]. The divergence time between C. hangiongensis and C. koreanus has been inferred from a mitochondrial DNA (mtDNA) phylogeny to be approximately 1.7-4.0 Mya [23]. Cottus hangiongensis populations in South Korea have experienced recent sharp declines because they are susceptible to changes in water environments caused by anthropogenic pressure, such as agricultural run-off, water pollution by sewage disposal, and impoundments for water retention [34,35]. Therefore, this species has been protected by the Korean government as a 'legally protected species II' since 2012 (endangered wild species class II), but C. koreanus has recently been removed from the Red List of Korean endangered species.
Several studies on Korean populations of C. hangiongensis and/or C. koreanus have been performed to understand their morphology [24], reproductive and spawning behaviours [27,36], feeding ecology [37,38], and capacity for natural hybridization [31,39]. Previous morphological surveys suggest natural hybridization is possible between these two species [31,39]. However, little attention has been paid to the population genetic structure of Korean populations of both species (but see [39]). In particular, the potential role of life history in shaping the population genetic structure of this species pair remains largely unexplored.
In this study, we examined whether there were differences in patterns and levels of genetic diversity and population connectivity, as well as geographic population structure between the amphidromous species C. hangiongensis, which has a pelagic larval stage, and the fluvial species C. koreanus, which has direct development. To this end, we used mtDNA control region sequences and seven nuclear microsatellite loci to analyse and compare the levels of within-population genetic diversity among five populations of C. hangiongensis and 10 populations of C. Fig. 1 The geographic distribution and sampling localities of C. hangiongensis and C. koreanus. a The entire geographical ranges of both species (C. hangiongensis in blue and C. koreanus in yellow). b Within South Korea, the geographic distribution of C. hangiongensis is limited to the east-flowing rivers (shaded areas in blue), but that of C. koreanus ranges from the Han River further south to the Nakdong River in inland areas (in yellow). Circles in blue denote sampling sites for C. hangiongensis (N = 6) and those in yellow indicate sampling localities for C. koreanus (N = 13). The latitude and longitude of each sampling location is shown in Additional file 1: Table S1. Population abbreviations as in text and Additional file 1: Table S1 koreanus in South Korea. Due to its greater capability to disperse, C. hangiongensis is expected to show much lower levels of geographic population structure, but a higher extent of within-population genetic diversity than C. koreanus. The results of this study will advance our understanding of how life history features influence the population structure and also provide conservation implications for these endangered species.

MtDNA control region sequencing
Genomic DNA was isolated using a DNeasy Blood and Tissue Kit (Qiagen, Hilden, Germany). A 465-and a 467-bp sequence of the mtDNA control regions (CR) of C. hangiongensis and C. koreanus, respectively, were amplified using the published forward and reverse primers L-Thr, H12Sr5 and internal primer LCCR [40]. Polymerase chain reaction (PCR) amplification was performed in a reaction volume of 15 μl comprising 1× PCR buffer, 25 μM of each dNTP (Bio Basic Inc., Markham, ON, Canada), 0.6 μM of each of the forward and reverse primers and 0.2 U of Taq polymerase (Thermo Fisher Scientific, Waltham, MA, USA). The following thermal cycling conditions were used: initial denaturation at 94°C for 1 min followed by 35 cycles of denaturation at 94°C for 1 min, annealing at 58°C for 1 min and extension at 72°C for 1 min, followed by a final extension at 72°C for 20 min in a 2720 thermal cycler (Applied Biosystems, Foster City, CA, USA). PCR products were checked on 2% agarose gels stained with RedSafe (iNtRon Biotechnology, Daejeon, Korea). The amplified PCR products were purified enzymatically with Exonuclease I (New England BioLabs, Ipswich, MA, USA) and Shrimp Alkaline Phosphatase (New England BioLabs). The purified mtDNA fragments were subjected to direct sequencing using the internal forward LCCR and the same reverse H12Sr5 primers as in the PCR, and the BigDye Terminator 3.1 Cycle Sequencing Ready Reaction Kit in an ABI 3730xl automated DNA sequencer (Applied Biosystems). The DNA sequences were edited using CHROMAS v2.01 computer software and aligned with Clustal Omega [41], implemented in BioEdit v7.2.5 [42], and finally verified by visual inspection.

Microsatellite genotyping
The seven published primer pairs, which were originally developed for the European bullhead species Cottus gobio [43,44], were used for microsatellite genotyping of our study species. Seven polymorphic nuclear microsatellite loci were chosen, including Cgo56, Cgo05, Cgo18, Cgo22 [43], Cott138, Cott207 and Cott112 [44]. The forward primers were labelled with fluorescent dye (FAM, NED, VIC and PET). PCR reactions were carried out as described for amplifying the mtDNA control region. PCR cycling conditions comprised an initial denaturation phase at 94°C for 1 min, 35 cycles of 94°C for 1 min, 57°C-64°C for 1 min, 72°C for 1 min and a final extension step at 72°C for 20 min. The PCR products were electrophoresed on an ABI 3730xl automated DNA sequencer (Applied Biosystems). Fragment sizes were determined to a ROX 500-bp size standard (ABI), as detected using GEN-EMAPPER software v5.0 (Applied Biosystems).

Statistical analyses Genetic diversity
To determine levels of mitochondrial diversity in both C. hangiongensis and C. koreanus, the numbers of polymorphic sites and haplotypes, haplotype diversity (h) and nucleotide diversity (π) were estimated for each population, as well as for the entire species using ARLE-QUIN v3.5 [45]. A rarefaction method was applied using CONTRIB v1.02 [46] to estimate haplotype richness (HR), corrected for unequal sample sizes among the samples. The HR estimates were calculated after excluding the samples of OK (N = 1) for C. hangiongensis, and WJ 1 (N = 9), GS (N = 7) and YG (N = 2) for C. koreanus, because of insufficient sample sizes (Table 1). These four samples were omitted from downstream population genetic analyses, with the exception of the mtDNA haplotype network and microsatellite-based STRUCTURE and factorial correspondence analyses (FCA). The OK sample was omitted from STRUCTURE and FCA analyses because it was comprised of only a single individual. To investigate relationships among the mtDNA haplotypes of each species, a haplotype network was constructed with HAPSTAR v0.7 [47], based on the minimum spanning tree generated from ARLEQUIN. To assess microsatellite diversity in each species, mean number of alleles per locus (N A ), expected (H E ) and observed (H O ) heterozygosity, inbreeding coefficient (F IS ), and allelic richness (AR) corrected for unequal sample sizes were estimated using GENEPOP v4.0 [48] and FSTAT v2.9.3.2 [49]. The presence of null alleles was assessed using MICROCHECKER v2.2.3 [50] with 1000 randomizations at the 95% confidence level. Genotypes at the seven microsatellite loci were tested for linkage disequilibrium (LD) separately for each species, and multi-locus tests for Hardy-Weinberg equilibrium (HWE) were undertaken for every sample of each species using GENEPOP. The 95% significance levels for every exact test for both LD and HWE were adjusted using a Bonferroni correction.
A recent reduction in effective population sizes (N e ) was tested using BOTTLENECK v1.2.02 [51], by applying 10,000 permutations for the two-phase mutation (TPM) model of microsatellite evolution. The TPM model was used because it has been shown to be the best fitting for a microsatellite dataset [51] and allows for multi-step mutations. A Wilcoxon sign-rank test was used to determine whether the populations of both species showed a significant heterozygosity excess, which would be expected with founder effects or recent population bottlenecks [52]. Contemporary N e was also calculated for each of the samples in both species according to the LD method implemented in NEESTIMATOR v2.01 [53].
Two independent Mann-Whitney U tests were conducted to test whether there were significant differences in the levels of mitochondrial (HR) and microsatellite (AR) diversities between populations of C. hangiongensis (N = 5) and C. koreanus (N = 10). The significant differences in levels of H E and H O between the two species were also statistically analysed, as done for HR and AR.

Spatial population genetic structure
To examine the spatial genetic differentiation between populations within each species, exact tests for population differentiation [54], as well as calculation of pairwise estimates of F ST [55] at both markers, were carried out using ARLEQUIN and GENEPOP, respectively. The 95% significance levels for the pairwise population comparisons were adjusted using a Bonferroni correction. In addition, we tested for isolation by distance (IBD) among samples separately for each species at both markers using the Mantel test in GENALEX v6.502 [56]. The geographic surface distance in kilometres between two sampling sites was obtained from the website (http://www.movable-type.co.uk/ scripts/latlong.html), based on the coordinate information (latitude/longitude) for each location.
We further analysed the spatial population genetic structure of each species using an individual-based Bayesian population assignment test in STRUCTURE v2. 3.1 [57], under an admixture model, and correlated allele frequencies with no a priori information on the geographic origins of the samples. We tested 10 iterations at each K = 1-18 (C. hangiongensis: number of samples [N] = 5; C. koreanus: N = 13), with 50,000 burn-in steps followed by 500,000 Markov chain Monte Carlo generations. The most probable number of clusters (K value) was estimated using the ΔK method implemented in the web-based tool Structure Harvester (http://taylor0.biol ogy.ucla.edu/structureHarvester/) [58], based on the rate of change in the log probability of data between successive K values. STRUCTURE analyses were also performed separately for populations of C. hangiongensis and C. koreanus. Finally, FCA, based on genetic relationships among individuals with multi-locus genotypes, was performed for both species using GENETIX v4.05.2 [59].
The levels of microsatellite diversity (AR) were significantly higher in amphidromous C. hangiongensis than in fluvial C. koreanus (C. hangiongensis: mean AR = 13.23; C. koreanus: mean = 5.02; Mann-Whitney U, P < 0.001; Fig. 2a; Table 1). The mean number of alleles across the seven loci (N A ) per population for C. hangiongensis was 14.92, ranging from 13.29 (UJ) to 16.57 (SCH), whereas the mean in C. koreanus was 4.83, ranging from 1.71 (YG) to 11.14 (PC 2) (Table 1) Table 1). However, the number of private alleles (PA) detected was approximately 3.2 times higher for C. koreanus (N = 74) than for C. hangiongensis (N = 23), despite more populations having been analysed for C. koreanus. The highest number of PA was identified in SCH (N = 7) for C. hangiongensis and in GS (N = 22) for C. koreanus. The F IS values ranged from − 0.004 (MC) to 0.058 (YYH) for C. hangiongensis and from − 0.291 (YY) to 0. 060 (JS) for C. koreanus. Based on our multi-locus tests for HWE expectations, only one population in C. hangiongensis (SCH) and five populations in C. koreanus (PC 1, JS, IJ, GP and GG) may be experiencing non-random mating (inbreeding or outbreeding) at the seven loci analysed ( Table 1). The estimated frequencies of null alleles at the seven loci were close to zero, ranging from 0.004 (Cgo18) to 0.094 (Cott112) for C. hangiongensis and from 0.048 (Cgo22) to 0.314 (Cott138) for C. koreanus, indicating a low probability of null alleles. Tests of LD across the seven loci were not statistically significant (P > 0.05) for either species after a Bonferroni correction, suggesting that all loci were independent markers. BOTTELNECK analysis revealed that C. hangiongensis had no sign of population bottlenecks in the five populations analysed (Additional file 1: Table S2). Similarly, only one population (CA) in C. koreanus had allelic distribution shifts (mode-shift), which are typically considered to be evidence for population bottlenecks or founder effects (Table S2). However, CA showed no significant heterozygosity excess using a Wilcoxon sign-rank test (P = 0.469). The LD method gave median estimates of effective population size (N e ) of 161.2 (95% confidence interval [CI]: 74.3-∞) and an infinite N e (95% CL: 686.2-∞) for the GN and SCH populations of C. hangiongensis, respectively (Table 2). However, the estimates of N e for the populations in C. koreanus were generally much lower, ranging from 3.8 (95% CI: 1.3-23.7) for GG to an infinite N e (95% CI: 46.6-∞) for GP ( Table 2).
Degree of spatial genetic structure between C. hangiongensis and C. koreanus The degree of spatial genetic differentiation (F ST ) between populations at both mtDNA and microsatellite markers was much lower for the amphidromous species C. hangiongensis than for the fluvial species C. koreanus ( Table 3). The pairwise estimates of F ST between five samples of C. hangiongensis at mtDNA and microsatellites ranged from − 0.026 to 0.019 and from − 0.002 to 0.050, respectively. Unexpectedly, we observed weak but statistically significant genetic differentiation between samples of C. hangiongensis at only microsatellites (except YYH vs MC) ( Table 3). The F ST values estimated between 10 samples of C. koreanus were all highly significant (P < 0.01) and considerably higher than those of C. hangiongensis, ranging from 0.285 to 1.000 at mtDNA and from 0.065 to 0.761 at microsatellites. The Mantel tests of both species showed a lack of a significant correlation between geographic (km) and genetic (F ST values) distances among the populations (mtDNA: C. hangiongensis, r = 0.01, P = 0.50; C. koreanus, r = 0.31, P = 0.45; microsatellites: C. hangiongensis, r = 0.50, P = 0. 07; C. koreanus, r = 0.15, P = 0.20). STRUCTURE analysis showed that the 18 populations analysed across both species are most likely to form seven genetically distinct clusters (K = 7) (Fig. 4a). This number of genetic clusters was determined by: LnP (D) = 148.671 and ΔK = 10.952. Among the seven groups identified, two groups comprised populations solely of C. hangiongensis with a signature of high levels of genetic admixture and five clusters constituted the13 populations of C. koreanus (Fig. 4a). Within C. koreanus, only the PC 2 population showed moderate levels of genetic admixture, and the remaining populations were barely admixed (except the CA and SC populations, which consisted of a single genetic integrity). Some individuals of C. koreanus showed a signal of genetic introgression from C. hangiongensis into these individuals, although population-level admixture between the species was very low (Fig. 4a).
When STRUCTURE analyses were performed separately for the two sculpin species, C. hangiongensis populations showed homogeneous distributions of individual genotypes, in which all individuals were assigned to approximately equal proportions of the inferred three genetic clusters (Fig. 4b), suggesting those populations are genetically indistinguishable. However, C. koreanus' eight genetically distinct clusters (K = 8) best fitted the data (Fig. 4b). In this reduced analysis, C. koreanus populations clustered similarly to the full analysis, except that the YG, GS and SC populations, which flow into the East Sea, formed two new genetic clusters (YG and GS, and SC; Fig. 4b). The eight genetic clusters observed were determined by: LnP (D) = 220.347 and ΔK = 14.934. In addition, the second most likely K for C. koreanus was equal to 10 [LnP (D) = 941.295 and ΔK = 3.381]. Similar to the results of STRUCTURE, the FCA of seven microsatellites also revealed a single genetically indistinguishable group within C. hangiongensis, but every population is its own group within C. koreanus (Fig. 5).

Differences in levels of genetic diversity between amphidromous and fluvial sculpins
We find that an amphidromous sculpin species C. hangiongensis, with a post-hatching larval phase has considerably Fig. 3 Haplotype networks of C. hangiongensis and C. koreanus based on mtDNA control region (465, 467 bp each). A total of eight and 24 haplotypes were detected in C. hangiongensis and C. koreanus, respectively. The size of the circle is proportional to number of individuals that belong to the respective haplotypes and small white circles denote intermediate haplotypes that are not present in our samples, but are necessary to connect all of the observed haplotypes to the network. Each node in the network represents a single mutational step between haplotypes irrespective of its length. Population abbreviations as in text and Additional file 1: Table S1 higher levels of within-population genetic diversities than a fluvial landlocked species C. koreanus with direct development. This was true for both mtDNA control region sequences and seven microsatellite loci, although only microsatellites were statistically significant (Fig. 2). The observed higher genetic diversity in the amphidromous sculpin species, relative to the fluvial sculpin, agrees with a previous study of allozyme markers that found that populations of amphidromous C. amblystomopsis had greater levels of heterozygosity than those of fluvial C. nozawae on Table 2 Estimates of contemporary effective population sizes (N e ) for five and 10 populations of C. hangiongensis and C. koreanus, respectively in South Korea based on linkage disequilibrium (LD) method using NEESTIMATOR v2.01 [53]  Population abbreviations as in text and Additional file 1: Table S1   Table 3 Spatial genetic differentiation (as indicated by pairwise estimates of F ST ) between five and 10 populations of C. hangiongensis and C. koreanus, respectively based on mtDNA control region sequences (below diagonal) and seven microsatellite loci genotypes (above diagonal) (    Hokkaido Island in Japan [3]. The elevated levels of withinpopulation genetic diversity observed for C. hangiongensis might reflect their larger N e [2,60], which is supported by our LD-based estimates of contemporary N e for populations of C. hangiongensis and C. koreanus (Table 2). Our N e estimates indicate that the median N e in C. hangiongensis (mean N e = 485.23 when the SCH and MC populations [N e = infinity] were excluded) is approximately 8.4 times greater than that in C. koreanus (mean N e = 57.80 when the CA, IJ and GP populations [N e = infinity] were excluded), providing direct evidence supporting our hypothesis. Amphidromous C. hangiongensis typically produces many small eggs (752-1376 per clutch; mean = 1005) in the lower reaches of rivers [24]. Newly hatched larvae float to the surface of the river as a phototactic response and drift downstream to the river mouth [25,33]. The planktonic larvae spend approximately one and a half months in coastal environments, in which they move with oceanic currents and return to the estuary and settle on the riverbed. This lifestyle may facilitate genetic exchange between geographically disconnected populations via larval dispersal (i.e. gene flow), at least partly contributing to the observed higher levels of withinpopulation genetic diversity in C. hangiongensis [10,61]. The home range of this amphidromous species at an adult life stage was observed to be no larger than 92 m [29]. By comparison, fluvial C. koreanus produces fewer but larger eggs (538-880 per clutch; mean = 744) in the upstream rivers, and hatched larvae undergo direct development and maintain a benthic life [24,62]. In particular, the home range of this fluvial species is found to have a lifetime migration distance within 10 m, according to a transponder telemetry-based survey [28]. These findings suggest extremely low dispersal capacity for C. koreanus that leads to restricted gene flow among populations over geographic scales of tens to hundreds of metres. Smaller N e along with a poorer dispersal ability makes C. koreanus populations more susceptible to genetic drift, resulting in the lower extent of withinpopulation genetic diversity [63]. A previous study [5] estimated and compared the level of microsatellite diversity between amphidromous and fluvial landlocked populations within Rhinogobius sp. in Okinawa Island, Japan, and found much greater genetic diversity in the amphidromous populations. As such, the higher extent of genetic diversity in C. hangiongensis is likely attributable to high levels of gene flow and larger N e resulting from its amphidromous life history with a pelagic larval phase. In contrast, the low genetic diversity in C. koreanus is presumed to be due to genetic isolation among individual populations resulting from its fluvial landlocked life history with direct development. The greater fecundity (approximately 1.4 times higher fecundity) of C. hangiongensis may also contribute to its higher intraspecific genetic variation [2].
Our results also imply that while genetic variation from samples at a locality in C. hangiongensis may represent a complete picture of the genetic diversity due to high gene flow from surrounding sites, genetic variation in C. koreanus may be partitioned among different localities due to very low gene flow [64].
Alternatively, amphidromous and fluvial species have undergone very different demographic population histories during the last glacial maximum, as sea level fluctuations could have affected these species differently [65]. If glaciations had an effect on genetic diversity, we would expect to see large population bottlenecks for both species. However, this is not true for either species, given that only one (CA) of ten populations in C. koreanus and none of the populations in C. hangiongensis showed evidence for population bottleneck (Table S2). Therefore, the observed differences in within-population genetic diversity between these species can be most likely attributed to their differing early life histories. Most populations of both species appear to be in HWE, except for the SCH population of C. hangiongensis and the PC 1, JS, IJ, GP and GG populations of C. koreanus, which suggests non-random mating at the loci tested. Somewhat ironically, the SCH population of amphidromous C. hangiongensis shows a genetic signal of inbreeding, albeit weak, whereas the IJ and GP populations of fluvial C. koreanus show a detectable signal of outbreeding (Table 1). In theory, small N e is likely to result in high rates of inbreeding, causing an increase in homozygosity [66]. Nevertheless, the limited dispersal and reported polygynous mating system of the fluvial sculpin may account for outbreeding within the IJ and GP populations. After spawning, the male chases the female partner away from its nest, while guarding the eggs and embryos [27]. Our finding that a critically endangered species, C. hangiongensis, has high genetic diversity is fairly good news for conservation and its future persistence. None of the five populations of this species analysed shows a genetic signal of a population bottleneck. On the other hand, C. koreanus might be more threatened than previously thought, as it shows lower levels of genetic diversity and limited connectivity when compared to C. hangiongensis.
Differences in patterns of population connectivity and spatial population structure between amphidromous and fluvial sculpins Our results from both mtDNA and microsatellite markers reveal noticeable differences in the magnitude and patterns of population connectivity between the two sculpin species that vary in the presence of a pelagic larval phase. Amphidromous C. hangiongensis, with a planktonic larval period, shows apparently less spatial population structure relative to the fluvial C. koreanus, with a benthic lifestyle. The lower spatial genetic variation of the amphidromous species is most likely due to the greater larval dispersal than the fluvial species, given their restricted ranges as adults (tens of metres). The magnitude of spatial genetic differentiation, as suggested by pairwise F-statistics, was markedly lower within C. hangiongensis than within C. koreanus (Table 3). The F ST values for C. koreanus were all highly significant, which is consistent with a previous population genetics study of this species with AFLP (Amplified Fragment Length Polymorphism) marker [39]. These results suggest that gene flow is restricted among the 10 populations of C. koreanus examined, and that this is most likely caused by limited dispersal ability. A recent study of the Manchurian trout Brachymystax lenok tsinlingensis from South Korea inhabiting watersheds that overlap like those of C. koreanus also showed significant genetic divergence among populations [67]. Examining another species pair in the genus Cottus that differs in the presence/absence of the post-hatching larval stage on Hokkaido Island in Japan, a previous study [3] also demonstrated that the amphidromous C. amblystomopsis populations are less structured than the fluvial C. nozawae. These previous findings also support the hypothesis that the divergent patterns of genetic structure between C. hangiongensis and C. koreanus may have been caused by differences in life history traits.
IBD analysis reveals no correlation between geographic and genetic distances among the populations in both species, indicating that the geographic proximity of populations is not responsible for the observed spatial population structure of either species. For amphidromous C. hangiongensis, geographic distance among localities does not seem to play a role in shaping the population structure over the geographical scales analysed (9-140 km), perhaps due to the effects of connectivity during larval phases promoted by the oceanic currents. However, a study of short ninespine stickleback Pungitius kaibarae without a larval phase at the east-flowing river basins in South Korea revealed strong inter-population genetic differentiation [68], which supports the notion that the absence of genetic structuring in C. hangiongensis is highly likely due to the presence of the larval phase. The lack of population structure (i.e. high population connectivity) observed for C. hangiongensis suggests that the juveniles of this species do not necessarily return to their home rivers where they hatched. For fluvial C. koreanus, its very low dispersal ability may mask the role of geographic distance that contributes to the population structure over the geographic scales we tested (1-177 km). Further studies on smaller geographic scales will be required to determine the spatial scales at which C. koreanus populations are structured.
The evidence that amphidromous C. hangiongensis comprises a single random mating population, but fluvial C. koreanus consists of genetically distinct individual populations is further supported by the results of our STRUCTURE and FCA analyses. While genotyped individuals of C. hangiongensis from five localities represent a single population with a high degree of admixture, C. koreanus individuals from thirteen sites were divided into five (or eight) genetically unique clusters (Fig. 4). However, the genetic groups within C. koreanus do not correspond well to river basins, such as the South Han River, the North Han River and east-flowing river basins, as suggested in a previous study of the Korean lenok [67].
Given the similarity of traits such as diet composition [37,38], mating system [27,69], spawning season [27,69] and adult's home range [28,29] between the two species, it is highly likely that the observed differences in the magnitude and pattern of genetic diversity and geographic population structure between C. hangiongensis and C. koreanus is driven by their divergent life histories [3,5]. These closely related species are very similar to each other, to the point where they are not fully ecologically or genetically isolated. We found evidence for hybridization on our microsatellite dataset, and we will evaluate the extent and outcomes of these events in future research projects by applying ecological and genomic approaches.

Conclusions
This study suggests that an amphidromous species with a planktonic larval phase, C. hangiongensis, shows higher levels of within-population genetic diversity than a fluvial species with direct development, C. koreanus. On the contrary, the fluvial sculpin exhibits much greater levels of geographic population structure than the amphidromous sculpin. While C. hangiongensis represents a single genetically homogeneous population, C. koreanus reflects genetically distinct local populations. We conclude that contrasting early life historiesincluding the presence or absence of the larval phasemay have led to divergent patterns of genetic diversity, population connectivity and spatial population structure between these two sister-species of sculpin. This study highlights the potentially important role of life history traits in the evolution of population genetic structure and suggests the need for further investigations into how life history features contribute to speciation more generally.

Additional file
Additional file 1: Table S1. Information of sampling localities, population codes, coordinate (latitude/longitude), and river basins. Table S2. Statistical tests for a recent bottleneck in each of the five and 10 populations of C. hangiongensis and C. koreanus, respectively from South Korea. P-values are based on the Wilcoxon test. Allelic frequency distribution shape was normal or shifted for mode-shift distortion. Population abbreviations as in text and Additional file 1: Table S1. (ZIP 26 kb)