Local environment‐driven adaptive evolution in a marine invasive ascidian (Molgula manhattensis)

Abstract Elucidating molecular mechanisms of environment‐driven adaptive evolution in marine invaders is crucial for understanding invasion success and further predicting their future invasions. Although increasing evidence suggests that adaptive evolution could contribute to organisms’ adaptation to varied environments, there remain knowledge gaps regarding how environments influence genomic variation in invaded habitats and genetic bases underlying local adaptation for most marine invaders. Here, we performed restriction‐site‐associated DNA sequencing (RADseq) to assess population genetic diversity and further investigate genomic signatures of local adaptation in the marine invasive ascidian, Molgula manhattensis. We revealed that most invasive populations exhibited significant genetic differentiation, low recent gene flow, and no signal of significant population bottleneck. Based on three genome scan approaches, we identified 109 candidate loci potentially under environmental selection. Redundancy analysis and variance partitioning analysis suggest that local environmental factors, particularly the salinity‐related variables, represent crucial evolutionary forces in driving adaptive divergence. Using the newly developed transcriptome as a reference, 14 functional genes were finally obtained with potential roles in salinity adaptation, including SLC5A1 and SLC9C1 genes from the solute carrier gene (SLC) superfamily. Our findings confirm that differed local environments could rapidly drive adaptive divergence among invasive populations and leave detectable genomic signatures in marine invaders.

genetic divergence among populations (Moran & Alexander, 2014;White et al., 2013). Additionally, both the rate and range of environmental changes during biological invasions are often orders of magnitude faster and greater than what species experience in natural processes (e.g., seasonal fluctuations) in their native habitats . In order to cope with various environmental challenges, some invasive species have considerable potentials for rapid adaptation (Batz et al., 2020;Lee, 2002). Such a common feature indicates that invasive species can provide promising materials to investigate how organisms adapt to changing environments Zhan et al., 2015). Indeed, empirical studies have recently evidenced that environment-driven microevolution, such as local adaptation, could play an important role in organisms' adaptation to different environments Ni et al., 2018). However, there remain knowledge gaps regarding local adaptation in marine invaders in their invasive range and genetic/genomic bases underlying local adaptation have not well been elucidated (Tepolt, 2015). In particular, it remains largely unexplored how local environmental conditions influence population genetic variation in a set of invasive populations. A better understanding of genetic/genomic bases of environmental adaptation will provide useful insights into evolutionary dynamics of invasive species in response to environmental changes, as well as predict their future invasion dynamics into new habitats (Jeffery et al., 2018;Prentis et al., 2008).
Marine organisms, especially sessile tunicates, are more vulnerable to be influenced by local environmental challenges, mainly owing to their sessile lifestyle (see the review by Zhan et al., 2015).
Different local environmental factors, such as water salinity, temperature, and dissolved oxygen stressors, can greatly affect their growth, development, reproduction, and geographical distributions.
Environmental heterogeneity could potentially lead to genetic divergence among populations, and thus investigating genomic signatures imposed by local environments would provide insights into mechanisms and dynamics of environmental adaptation (Attard et al., 2018;Hecht et al., 2015). However, it remains challenging to disentangle the effects of environmental selection from population demographic history in wild populations (Ferchaud & Hansen, 2016), mainly because population demography-associated processes, such as genetic drift and gene flow, are often accompanied by natural selection and further may obscure, or even eliminate, signals of natural selection in introduced populations (Li et al., 2012;. Genetic drift may contribute to random changes of allele frequency, further complicating the identification of loci under selection (Lotterhos & Whitlock, 2014). For many sessile tunicates, although they have limited natural dispersal ability owing to short pelagic larval durations and sessile adults, human-mediated gene flow is likely to drive genetic homogeneity Zhan et al., 2010). Thus, it poses challenges for the accurate estimation of population genetic divergence when merely using "lowresolution" genetic markers (e.g., a limited number of microsatellites).
Additionally, spatial distributions may confound the environmental effect in the analysis of local adaptation if environmental factors are spatially autocorrelated (Benestan et al., 2016;Rellstab et al., 2015).
It is therefore important to choose a promising model system and state-of-the-art techniques to investigate the mechanisms of local adaptation in marine invasive species.
The notoriously invasive ascidian, Molgula manhattensis, which is native to the east coast of North America, has invaded large geographical scales globally and caused significantly negative ecological, environmental, and economic impacts in invaded marine and coastal ecosystems (Chen, Li, et al., 2017;Haydar et al., 2011). In China, M. manhattensis was firstly recorded at Tanggu Harbor in Tianjin in 1976 (Zheng, 1995). Since then, it has widely spread along the Chinese coast, especially in Bohai and Yellow seas, where this species has been reported as the primary fouling organism and caused huge damages (Chen, Li, et al., 2017). Owing to its biological features including a short pelagic larval phrase and sessile adults as well as a high self-recruitment rate, this species is characterized by a low natural dispersal ability, and their range expansion in China is likely to be transported with aquaculture transfers and/or via shipping as fouling organisms. Across its wide distribution range along the Chinese coast (spanning > 15 degrees of latitude), this species has encountered heterogeneous environmental conditions, especially water temperature and salinity, which may impose strong selective pressures during range expansions. In addition, some biological characteristics, such as high fecundity, rapid growth rate, and short life cycle, can facilitate this ascidian to rapidly colonize in novel environments (Chen, Li, et al., 2017). Altogether, the rapid and widespread colonization by M. manhattensis along the Chinese coast provides us a promising model system to study mechanisms of local adaptation driven by divergent local environments, particularly by water temperature and salinity.
As shown in one previous study with one mitochondrial COI marker and 12 nuclear microsatellites, they revealed weak genetic structure at large geographical scales but significant genetic differentiation among some neighboring M. manhattensis populations along the Chinese coast (Chen, Li, et al., 2017). We propose the effect of local environmental selection on significant population genetic differentiation, but direct evidence remained unavailable, owing to the limited number of genetic markers used. In addition, the observed inconsistent results of population genetic diversity by mtDNA and nuclear markers were likely derived from their different inheritance natures responding to natural selection (Chen, Li, et al., 2017). As suggested by a growing number of related studies, the genome-wide survey would provide more power in uncovering subtle population genetic structure (Benestan et al., 2015;Blanco-Bercial & Bucklin, 2016) and more importantly detecting signatures of environment-driven local adaptation Waterhouse et al., 2018).
With technical advances in next-generation sequencing technologies, restriction-site-associated DNA sequencing (RADseq) provides a relatively economical and efficient approach for genome-wide genotyping of DNA polymorphisms (Andrews et al., 2016). Population genomics offer an opportunity to investigate genetic structure and gene flow more accurately and comprehensively, further contributing to the differentiation of locus-specific effects (e.g., selection, mutation) from genome-wide effects (e.g., gene flow, genetic drift, bottleneck) (Luikart et al., 2003). In addition, the utilization of this robust genotyping method, as well as efficient theoretical and analytical approaches such as appropriate genome scan methods, can contribute to the investigation of genomic basis of local adaptation, particularly the relative contribution of environmental factors to population genetic divergence (de Villemereuil et al., 2014;François et al., 2016).
Using the RADseq on M. manhattensis populations collected from a wide geographical range along the Chinese coast, we aim to test the hypothesis that differed local environmental factors could drive adaptive genetic divergence among populations and leave detectable signatures in the genome of marine invaders. Based on a neutral single nucleotide polymorphisms (SNPs) dataset, we firstly assessed population genetic variation and gene flow patterns. Subsequently, we combined population differentiation (PD) and environmental association (EA) tests to search for loci potentially under selection.
We further conducted multiple analyses to investigate the effects of environmental variables in driving putatively adaptive genetic variation. Finally, with the aid of the transcriptome sequencing, we annotated a set of functional genes potentially related to environmental adaptation.

| Sampling and genotyping
A total of 118 M. manhattensis individuals were collected from eight sites along the Chinese coast in 2013 (Table 1; Figure 1).
Specifically, QJZ was sampled in the Yueqing Bay, which is a typical semi-enclosed bay with low salinity (approximately 20‰) due to inflows of rivers. The others were collected along the inshore of Chinese coast (Figure 1). Six environmental factors associated with sea surface salinity and temperature (among 1955-2012) were obtained from the World Ocean Atlas 2013 of the National Oceanic and Atmospheric Administration (www.nodc.noaa.gov/OC5/SELEC T/woase lect/woase lect.html), including three sea surface salinity matrices: annual average salinity (AveS), the highest monthly average salinity (MaxS), and the lowest monthly average salinity (MinS), and three sea surface temperature matrices: annual average temperature (AveT), the highest monthly average temperature (MaxT), the lowest monthly average temperature (MinT) ( Table 1).
Total genomic DNA was extracted using the DNeasy Blood and Tissue Kit (69,506, QIAGEN). The quality and quantity of isolated genomic DNA were examined by a NanoDrop 2000 spectrophotometer (Thermo Fisher Scientific, USA). The 2b-RAD library construction was conducted as described by Wang et al., (2012) with minor modifications. Briefly, we added NNRW adapters (N refers to any nucleotide, R refers to A or G, W refers to A or T) in order to remove PCR duplicates, which can reduce downstream genotyping errors (Dixon et al., 2015). In order to pool different libraries into a single lane for sequencing, we adopted the multiplexing 2b-RAD library method by combining 12 barcoded adapters with 24 indexed PCR TA B L E 1 Sampling sites information and population genetic variation based on 3,616 putatively neutral single nucleotide polymorphisms (SNPs)  primers (Guo, Yuan, et al., 2014). For all eight populations, we sequenced three lanes of single-end sequencing (SE50) on the Illumina Hiseq 4,000 platform (Illumina, USA).
After the high-throughput sequencing, these libraries were de-multiplexed based on barcodes and further removed duplicates using a modified perl script (Gao et al., 2018). In addition, the terminal 2 bp nucleotides were removed from each read in order to avoid possible artifacts derived from the ligation sites.
After trimming and quality filtering with the "process_radtags" module in STACKS v.2.1 (Catchen et al., 2013), we obtained highquality reads with the uniform length of 32 bp. Filtered reads were then implemented with STACKS for de novo genotyping due to the lack of a reference genome of M. manhattensis. After parameter testing, the optimal parameters set in "denovo_map.pl" pipeline program in STACKS were as follows: minimum depth of three to create a stack (m = 3) in ustacks, maximum number of two mismatches allowed between stacks (M = 2) in ustacks, two mismatches allowed between sample loci within a catalog (n = 2) in cstacks and the default for sstacks to match stacks of all samples in the population against the catalog (Andrews et al., 2018).
The "populations" programs in STACKS were used to perform SNP filtering and produce a series of population genetics statistics. To minimize the occurrence of false SNPs, we performed the following strict filtering options as suggested by Rochette and Catchen (2017). Firstly, we discarded RAD tags with stacks depth less than four (m = 4), retained bi-allelic SNP per locus with minor allele frequency (MAF) > 0.05 to reduce false-positive SNP calls (Rochette & Catchen, 2017). Secondly, we screened for SNPs genotyped in at least 75% of all populations (p = 6) and 75% of individuals in a population (r = 0.75). Thirdly, we performed the "write_single_snp" option with only the first SNP per RAD tag retained in the final set to avoid physical-linked loci (Andrews et al., 2018). We finally removed SNPs and individuals with more than 20% missing data.

| Identifying loci potentially under selection
The main objectives of our study are firstly to provide a dataset of putatively adaptive loci, and then test their relevance to local environmental factors through genotype-environment matching analysis and RDA analysis, and finally conduct blast alignments to annotate candidate genes under selection. Therefore, our final panel of putatively adaptive loci was composed of loci identified by at least two of three genome scan methods, which can maximize the power of the analysis while reducing the FDR, especially for invasive species with complex demographic histories in different landscapes (Dalongeville et al., 2018;. Here, we adopted two population differentiation (PD) approaches (BAYESCAN, OUTFLANK) and one environmental association (EA) approach (LFMM) to test for natural selection. Firstly, we identified loci potentially under selection (i.e., outlier SNPs) using a Bayesian approach in BAYESCAN v.2.1 (Foll & Gaggiotii, 2008). The parameter settings of BAYESCAN analysis were as follows to reduce the number of false positives: 20 pilot runs with a length of 5,000 iterations after a burn-in of 50 000 steps. Candidate outlier loci were defined as those SNPs with false discovery rate (FDR) lower than 5% (q-value < 0.05). Another PD approach, OUTFLANK, which estimated the distribution of F ST for neutral loci based on the trimmed distribution of F ST values, was adopted to detect F ST outliers using the "OutFLANK" package in R . This method has been reported with a low false discovery rate and a high power in populations with complex demographic histories. We ran F I G U R E 1 Sampling localities of Molgula manhattensis along the Chinese coast. MinS, the lowest monthly average salinity (see more details in Table 1) OUTFLANK with the following options: RightTrimFraction = 0.05, LeftTrimFraction = 0.05, Hmin = 0.1, q-threshold = 0.05. As PD analyses just used the "blind" genome scanning without incorporating environmental variation (Rellstab et al., 2015;Villemereuil & Gaggiotti, 2015), we further employed one environmental association (EA) approach, a latent factor mixed-effect model (LFMM) in the R package "LEA" (Frichot & François, 2015). LFMM analysis tests correlate between loci and environmental factors while correcting for unobserved latent factors, which has been shown to be robust in complex models (Caye et al., 2019;Frichot et al., 2013). Firstly, we estimated the number of latent factors using a cross-entropy criterion with "snmf" function in LEA package with 10 repetitions for each K value between 1 and 8. Secondly, we performed an environmental association test using the LEA function "lfmm" with the best number of latent factor according to "snmf" results. In order to reduce the collinearity among environmental variables (Pearson'r test, Table S4), we summarized the environmental variables by running a principal component analysis (PCA) for three temperature variables and three salinity variables, respectively, and retaining only higher order PCs (i.e., PC1-T, PC1-S). The "lfmm" function was conducted with 100 000 iterations, 10 000 burn-in cycles and 5 replicates. We then combined z-scores over five runs and used adjusted p-values using the genomic control method implemented in LEA function "lfmm.pvalues", and this recalibration procedure could decrease the false discovery rate in the LFMM test (Frichot & François, 2015).
We finally obtained the list of candidate SNPs using the Benjamini-Hochberg procedure with a false discovery rate at 0.05.

| Neutral genetic variation and gene flow analyses
Natural selection can potentially influence estimates of neutral genetic diversity and gene flow (Luikart et al., 2003;Waterhouse et al., 2018), and we therefore chose to build a reliable "neutral dataset" by removing all potential outliers obtained from three approaches (BAYESCAN, OUTFLANK and LFMM; see below). We then used PLINK v.1.90 (Purcell et al., 2007) to test for departure from Hardy-Weinberg equilibrium (HWE, p <.01) at putatively neutral loci. Genome-wide population genetic diversity, including observed heterozygosity (Ho), expected heterozygosity (H E ), and nucleotide diversity (P i ), was calculated with the Populations program in STACKS. We further tested whether these populations had experienced recent genetic bottleneck using BOTTLENECK v.1.2.02 (Piry et al., 1999) under the two-phased mutation model (TPM) with 90% single-step mutations. Significant heterozygosity excess was evaluated by Wilcoxon signed-rank test with 1,000 iterations (Yang et al., 2014). These 100 random loci were sampled from the neutral dataset with 1,000 bootstrap iterations.
Population genetic differentiation was examined by calculating the pairwise F ST values with 10 000 permutations in ARLEQUIN.
The significance level of pairwise F ST values was adjusted after sequential Bonferroni correction. Population genetic structure was firstly examined by the discriminant analysis of principal components (DAPC) in the R package "ADEGENET" (Jombart, 2008). This multivariate method uses K-means of principal components and Bayesian information criterion (BIC) to infer the best-supported number of genetic clusters. We used the "optim.a.score" function to determine the number of retained principal components (PCs) and discriminant functions (DAs). To complement the DAPC analysis, we also performed one admixture analysis with the Bayesian clustering in STRUCTURE v.2.3 (Pritchard et al., 2000). For the STRUCTURE analysis, we performed 10 independent runs for each K (K was from one to eight), with 5,000,000 Markov chain Monte Carlo iterations preceded by 500,000 burn-in steps. The optimal K value was identified using the ΔK method of Evanno et al., (2005) in STRUCTURE HARVESTER (Earl, 2012), and we used DISTRUCT v.1.1 (Rosenberg, 2004) to visualize population assignments.
We estimated the recent migration rate between populations based on all neutral loci using the modified BayesAss analysis "BA3-SNPs" (Mussmann et al., 2019). The analysis was run for 10 000 000 iterations after a burn-in of 1 000 000 steps and sampled every 100th iteration. Mixing parameters for migration rate (m), allele frequency (a), and inbreeding coefficient (f) were adjusted to ensure the final MCMC acceptance rate between 0.35 and 0.45 for each variable (Mussmann et al., 2019). A rough 95% confidence interval (CI) was constructed as (mean ± 1.96 * sdev).

| Redundancy analysis and variance partitioning analysis
Redundancy analysis (RDA) and variance partitioning analysis (VPA) were used to investigate the relative contribution of spatial distribution and environmental factors to both putatively neutral and adaptive genetic variation among populations. Because both genetic variation and environmental factors may be spatially autocorrelated (Rellstab et al., 2015), we considered the effect of spatial distribution on population structure using the principal coordinates of neighbor matrices (PCNM) analysis (Benestan et al., 2016;Borcard and Legendre 2002). Spatial factors were assessed by the Euclidian distance matrix with the transformed latitude and longitude of each sampling site with the PCNM method. The PCNM analysis followed the procedure as described in Benestan et al., (2016), using the PCNM package in R. Before the RDA analysis, six environmental factors were firstly standardized using the "decostand" function with the standardize method implemented in the R package "vegan". We ran a PCA with "prcomp" procedure on six environmental factors to reduce their collinearity. The response variables were characterized with the minor allele frequency (MAF) of putatively adaptive and neutral loci calculated in PLINK v.1.90, respectively. The MAF data were detrended with the Hellinger method implemented in "vegan" package in R. We identified significant explanatory variables using "forward.sel" function in R. The partial RDA (pRDA) and variance partitioning analysis (VPA) were then used to distinguish the effect of environmental variables from spatial factors. The pRDA analysis was implemented with the R package "vegan". Variance partitioning analysis (VPA) was performed to assess the pure effects of environments and geography as well as their joint effects using the "varpart" function with 1,000 permutations.

| Transcriptome sequencing and gene annotation
Due to the limited genomic resources available for M. manhattensis, we used the transcriptome of this species as a reference to annotate outlier loci, which was a cost-effective method for nonmodel species without reference genomes (Benestan et al., 2016;Lamichhaney et al., 2012;Thomas et al., 2017).

| SNP genotyping
After raw data processing, a total of 769 983 RAD tags were obtained from three lanes of Illumina high-throughput sequencing (Table S1). The depth of RAD tags for each sample varied from 10 to 22 x. After the strict SNP filtering steps in STACKS, we finally obtained a panel of 6,635 high-quality genome-wide SNPs for subsequent analyses (Table S1).

| Neutral genetic variation and gene flow
After removing all potential non-neutral SNPs (i.e., a total of 2,335 loci based on three genome scan methods), 3,616 candidate neutral SNPs with Hardy-Weinberg equilibrium were finally retained, con-  Table S2). For the whole SNP dataset, eight populations were assigned into five genetic clusters (optimal K = 5; Figure 2a (Table S3).

| Candidate adaptive loci detection
For the BAYESCAN and OUTFLANK analyses, we identified 258 and 52 loci, respectively (Figure 3a). For the LFMM, the best number of latent factors was set to five with the minimal cross-entropy according to "snmf" results ( Figure S2), which was also consistent with the DAPC and STRUCTURE results (Figure 2a). The first axis of PCA for temperature-related variables (PC1-T) and salinity-related variables (PC1-S) both explained over 96% of the total variation (PC1-T: 96.96%; PC1-S: 96.22%). After a recalibration procedure with the genomic control method, we finally obtained a total of 2,140 loci, including 1,005 associated with PC1-T and 1,780 associated with PC1-S ( Figure 3b). Interestingly, 645 SNPs (30.14%) were shared among these two groups of environmental variables, whereas the others were group-specific, including 360 (16.82%) temperaturespecific loci and 1,135 (53.04%) salinity-specific loci. Collectively, 109 SNPs identified by at least two of three genome scan methods were considered as putatively adaptive loci for subsequent analyses. We found that the distribution of allele frequency at putatively adaptive loci largely varied among populations. Based on the allele frequency of 109 candidate loci, eight populations were divided into five groups ( Figure S3), which was consistent with the genetic patterns based on putatively adaptive SNPs as revealed by the DAPC analysis ( Figure S1b). Several loci, such as 8700_26, 32625_5, 17568_3, 18471_28, exhibited higher allele frequencies in the two south populations (ND and XIM) when compared to the others ( Figure S3). For another south population QJZ, allele frequency at some loci, such as 31311_0, 30284_22, 4585_28, was the highest (allele frequency > 0.8), whereas a large number of loci (79 out of 109) had an allele frequency of zero ( Figure S3).  Figure 4a). When partitioning the relative effect of environmental variables and spatial distribution on putatively adaptive genetic variation, the variation purely explained by environments was higher (15.05%) than that explained by spatial distribution (11.75%) (Table 2; Figure 4b). We  Table 2).  (Table 3). For instance, two nonsynony-  (Table 3). Another nonsynonymous SNP 19292_8 (G/A, Glu/Lys) is located in the titin (TTN) gene, which plays a role in the assembly and functioning of muscles as well as chromosome condensation and segregation during mitosis (Table 3). Among these candidate SNPs with synonymous substitutions, the SNP 8499_18

| Gene annotation of candidate adaptive loci
is located in the solute carrier family 9 member C1 gene (SLC9C1, also called sodium-hydrogen exchanger 10 gene), which is a member of the Na + /H + exchanger family involved in the electroneutral exchange of Na + and H + . Another SNP 37385_3 was also located in one gene of the solute carrier (SLC) subfamily, the solute carrier family 5 member 1 (SLC5A1). This gene can get involved in sodium ion transport, glucose transmembrane transport, and intestinal hexose absorption (Table 3).

| D ISCUSS I ON
In this study, we employed genome-wide SNPs to assess popula-

| Environment-driven adaptive genetic divergence
Dissecting the effect of environmental conditions on population genetic variation across introduced populations would provide us a better understanding on how invasive species respond to different environments during range expansions (Papacostas et al., 2017;White et al., 2013). For newly introduced populations, demographic processes, such as genetic drift and gene flow, might confound the effect of environmental factors when investigating signals of natural selection (Ahrens et al., 2018;Lotterhos & Whitlock, 2014 (Chen, Li, et al., 2017). Gene flow could play a role in influencing population divergence and may weaken, or even completely eliminate, signatures of local adaptation (Lenormand, 2002). We detected low recent migration rates, which is consistent with the observed significant genetic divergence among most populations. The observed population genetic heterogeneity at large geographical scales and low gene flow among distant populations were in contrast to the findings in Chen, Li, et al., (2017), which uncovered low population differentiation and high genetic similarity among several geographical-distant sites (i.e., DAD versus ND). This discordance of population genetic pattern may result from the adopted molecular makers among these two studies. Compared with population genetic studies with a limited number of genetic makers, genome-wide SNPs here could contribute to the detection of subtle genetic structure and accurate estimation of gene flow (Benestan et al., 2015;Blanco-Bercial & Bucklin, 2016).
M. manhattensis has limited natural dispersal ability, mainly owing to a short pelagic larval phrase, sessile adults, and a high selfrecruitment rate (Chen, Li, et al., 2017). Considering that M. manhattensis is commonly found in the aquaculture ponds and facilities with a high density (Zheng, 1995), human-mediated dispersal via frequent aquaculture transfers might exist as suggested by the previous study (Chen, Li, et al., 2017) and studies with other invasive ascidians (i.e., Zhan et al., 2012;Botrylloides violaceus in Bock et al., 2011;Botryllus schlosseri in Bock et al., 2012). Therefore, the significant genotype-environment matching pattern in this study suggests that the underlying gene flow among populations is low enough to allow for local adaptation.

Ciona intestinalis in
In our study, we revealed a large number of candidate environment-related loci (i.e., 2,140 in LFMM) and the key role of environmental factors in affecting putatively adaptive genetic variation.
Similar genotype-environment matching patterns have been recently uncovered in a variety of marine organisms, including the American lobster (Benestan et al., 2016), Chinese sea bass (Zhao et al., 2018), and greenlip abalone . For example, Sandoval-Castillo et al., (2018) detected population adaptive divergence in Haliotis laevigata along the coast in southern Australia and the observed adaptive variation was significantly correlated with the minimum sea surface temperature and oxygen concentration.
Natural selection derived from differed environments can contribute to the genotype-environment matching pattern with an indication of local adaptation (Rellstab et al., 2015;Villemereuil & Gaggiotti, 2015).
However, for invasive species, both adaptation occurring before introduction in the native range (pre-introduction adaptation) and after introduction in the invasive range (postintroduction adaptation) may create the observed environment-associated genetic variation in recently introduced populations Colautti & Lau, 2015). Previous studies suggested that local adaptation and invasion success would be more likely to occur with pre-adapted genotypes if there was a close match between source and recipient environments, such as anthropogenically induced adaptation to invade (AIAI) (Hufbauer et al., 2012). However, such parallel introduction alone is less likely to establish the obvious genotype-environment

TA B L E 3 (Continued)
matching pattern and local adaptation in diverse introduced regions in our study (Colautti & Lau, 2015). In addition, postintroduction evolution, that is, contemporary evolution, can urge non-native species to rapidly respond to different environments and increase invasion success Prentis et al., 2008). For example  detected rapid evolution after a recent introduction of C. robusta to harsh environments in the Red Sea. Therefore, further efforts by sampling the source regions of M. manhattensis are supposed to be conducted to distinguish the relative contribution of these two categories of adaptation in driving the observed genotypeenvironment matching pattern (Colautti & Lau, 2015).
Our sampling range spans more than 15 degrees of latitude along the Chinese coast and the sampling sites are characterized with contrasting environmental conditions, such as salinity and temperature.
Multiple lines of evidence in our study suggest that environmental conditions, especially salinity-related variables, significantly contribute to the adaptive genetic variation among populations. Our RDA results revealed that environmental variables accounted for more in explaining adaptive genetic variation when compared with spatial distribution (Table 2; Figure 4b). Furthermore, salinity-related variables were significantly selected and provided a higher proportion than temperature-related ones (Table 2). Compared with temperature-related variables, a larger number of candidate loci was also significantly associated with salinity (PC1-S), most of which belonged to the salinity-specific ( Figure 2b). As shown in the variance partitioning analysis, there was a large proportion unexplained by the environmental variables (PC1) and spatial variables (V2), suggesting that many other biotic and abiotic variables, such as dissolved oxygen, pH, various pollutants, and biological interactions, may also affect the observed adaptive genetic variation. Seawater salinity stressors, that is, high and low salinity, can both directly and indirectly influence various aspects of marine organisms, such as reproduction, osmoregulatory, metabolic processes, and immune system response (Renborg et al., 2014). Empirical evidence revealed the link between physiological and metabolic traits with salinity challenges in M. manhattensis. For example, Cao et al., (2012) documented that oxygen consumption was significantly correlated with salinity conditions in M. manhattensis and a low value was detected at high salinity. Salinity stressors can also largely reduce the ingestion rate of this species (Cao et al., 2012). In addition, different salinity conditions have an effect on the developmental process of tunicates, such as C. intestinalis, especially for early life stages including embryos and larvae (Renborg et al., 2014). Likewise, changes of DNA methylation profiles associated with responses to different salinity have been evidenced in other invasive ascidians (e.g., Huang et al., 2017;Ni et al., 2018).

| Candidate key genes potentially involved in environmental adaptations
Uncovering candidate genes potentially for adaptation would provide direct evidence for genetic basis of adaptive evolution in response to changing environments . Given the short length of tags generated by the RADseq approach and limited genomic resources of nonmodel species, the use of transcriptome as a reference has been demonstrated to be a credible and cost-effective way to annotate outlier loci (Ahrens et al., 2018;Benestan et al., 2016;Christmas et al., 2016;Lamichhaney et al., 2012). However, those loci located in noncoding regions could fail to get annotated but may be important to selection. For example, candidate variants in gene promoter regions would regulate gene expression that influence phenotypic variation (Grubert et al., 2015). In addition, the RADseq approach usually only collects a fraction of the genome and we may miss those loci or genes that are not closely linked to these identified SNPs in the genome (Lowry et al., 2017;Hoban et al., 2016). These technical issues are major obstacles when using the RADseq and transcriptome approaches to investigate the genomic basis of local adaptation. In order to avoid false positives in this study, we used extremely strict criteria to identify 14 protein-coding genes with nonsynonymous/ synonymous polymorphisms, and these genes are candidates for further detailed investigation of salinity adaptation mechanisms in marine invasive species.
Two candidate synonymous SNPs 37385_3 and 8499_18 are located in the SLC5A1 and SLC9C1 genes, respectively, both belonging to the solute carrier (SLC) gene superfamily. Adaptive response to salinity stressors of marine species requires the effective ion transport and osmoregulatory processes (Fiol & Kültz, 2007). The SLC5A1 gene is a member of the Na + /glucose cotransporter gene family, which can mediate the transport of glucose and structurally related substances across cellular membranes by Na + cotransport. This gene was reported to get involved in salinity adaptation and osmoregulation in the liver of spotted sea bass (Zhang et al., 2017). Another candidate gene, the sodium/hydrogen exchanger 1 (SLC9C1), encodes Na + /H + exchanger (NHE) protein to regulate intracellular pH and acid-base homeostasis. Multiple studies also indicated the crucial role of sodium/ hydrogen exchanger genes in osmoregulation in a variety of teleost species (Verri et al., 2012). In addition, we also obtained several functional genes associated with immune responses to cope with environmental stressors, such as the integrin alpha-L (ITGAL), multidrug resistance protein 1 (MDR1), and DNA repair protein RAD50 (RAD50). Environmental pressures such as temperature and salinity challenges can make organisms more vulnerable to pathogen infections and diseases, and these tunicates reply on the innate immunity to cope with these challenges (Shida et al., 2003).
For instance, the integrin alpha-L (ITGAL) has been reported as a receptor of the complement component 3 (C3)-like molecule in the complement system during innate immune response of Ciona intestinalis (Shida et al., 2003) and Halocynthia roretzi (Miyazawa et al., 2001). Further research should be conducted to confirm their adaptive significance and functional mechanisms in salinity adaptation by using gene expression analyses or directly performing genetic manipulation (i.e., RNA interference, CRISPR/Cas9 gene) after challenge experiments.

| CON CLUS ION
Utilizing a wide geographical coverage for M. manhattensis populations along the Chinese coast, we found that local environmental factors, that is, salinity-related variables, might drive putatively adaptive genetic variation among varied invasive populations. Furthermore, we identified a set of putatively adaptive loci and functional genes, which are good candidates in studying the evolutionary dynamics of environmental adaptations in marine species. Collectively, our findings confirm the hypothesis that differed local environments could drive adaptive genetic divergence among populations and leave detectable signatures in the genome of marine invaders. Our results suggest that local environment-driven adaptation might contribute to range expansions and invasion success into varied environments. We further highlight the necessity to incorporate local adaptation dynamics into the prevention and control strategies of biological invasions.

CO N FLI C T O F I NTE R E S T
The authors declare no conflict of interest. Writing-review & editing (lead).