Data analyses 1. Genetic diversity

  number of mismatches allowed between loci to build the catalogue is defined by the ‐n parameter. 

 

The optimal values for the different parameters in the denovo_map.pl must be chosen with an appropriate methodology, because they frequently have a significant effect on the building and quality of the resulting loci. Values for the -M and -n parameters, were identified following the r80 optimization approach (Paris et al., 2017) with a subset of one individual per population (10). The de novo pipeline was run six times with different combinations of these parameters

(M1n1, M2n2, M3n3, M4n4, M5n5, M6n6). Paris et al. (2017) found that the approach when setting the -M and -n parameters to the same value obtained the best results. Setting -m parameter as 3 (Table 8) seems to be an optimal value for other studies based on different biological systems (Paris et al., 2017; Rochette

& Catchen, 2017). The r80 optimization method includes a

filtering step using the population program in Stacks, consisting of retaining only loci shared among a minimum of 80% of samples across all populations. The

number of polymorphic loci (Figure 27) shared among at least 80% of samples across all populations were compared between runs. After optimal parameters were identified (M=n=3), the de novo pipeline was run again with optimal settings. The final output was filtered using population setting the -R parameter as 0.80, the minimum number of populations a locus must be present in so as to process this locus (-p parameter) as 10 and with a minimum allele frequency (MAF) of 0.05. Only the first single SNP per RAD tag was called using –write-single-snp to reduce linkage disequilibrium effects. The SNP dataset was exported to Bayescan (Foll & Gaggiotti, 2008) format using PGDSpider version v 2.1.1.5

Program Parameters Description Parameter

trials

ustacks

Assembled identical raw sequence reads into Stacks

-m The minimum number of raw reads

required to form an initial stack 3 -N The maximum nucleotide distance

allowed between stacks and secondary reads

Default (M+2) -M The maximum nucleotide distance

allowed between stacks 1-6 cstacks

Creates a catalogue of consensus loci by compiling stacks identified across

individuals.

-n The number of mismatches allowed

between loci to build the catalogue 1-6 Table 8. Parameters chosen for the de novo pipeline in Stacks (Catchen et al., 2013). 

-3000 (R80) obtained for different parameter trials in stacks.

best results. Setting -m parameter as 3 (Table 8) seems to be an optimal value for other studies based on different biological systems (Paris et al., 2017; Rochette & Catchen, 2017).

The r80 optimization method includes a filtering step using the population program in Stacks, consisting of retaining only loci shared among a minimum of 80% of samples across all populations. The number of polymorphic loci (Figure 27) shared among at least 80% of samples across all populations were compared between runs. After optimal parameters were identified (M=n=3), the de novo pipeline was

run again with optimal settings. Figure 27. Plot of the number of new polymorphic loci (R80) obtained for different parameter trials in stacks.

The final output was filtered using population setting the -R parameter as 0.80, the minimum number of populations a locus must be present in so as to process this locus (-p parameter) as 10 and with a minimum allele frequency (MAF) of 0.05. Only the first single SNP per RAD tag was called using –write-single-snp to reduce linkage disequilibrium effects. The SNP dataset was exported to Bayescan (Foll &

Gaggiotti, 2008) format using PGDSpider version v 2.1.1.5 (Lischer & Excoffier, 2012) and to Admixture v1.3.0, (Alexander et al., 2009) format using PLINK v1.90 (Chang et al., 2015).

3.2.3.4. Data analyses 3.2.3.4.1. Genetic diversity

Population specific genetic parameters were calculated for the dataset with first-single SNPs. Using STACKS v4.2 (Catchen et al., 2013), number of private alleles (PA), nucleotide diversity (π), observed heterozygosity (Ho), expected heterozygosity (He), and inbreeding coefficient (F_IS) were obtained for all positions and only for variant position. Moreover, the F_ST parameter that estimates genome-wide divergence between populations was calculated with a p-value correction (p<0.001) to exclude non-statically significant F_ST measures. Hierfstat R package (Goudet & Jombart, 2015) was used as another approach to calculate some genetic diversity parameters and corroborated the patterns of variability among these populations. Ho and He were obtained using the basic.stats function and population specific F_IS and F_ST were determined using the boot.ppfis and boot.ppfst functions respectively with a bootstrap confidence interval of 1,000 replicates. Pairwise F_ST following Weir & Cockerham (1984) was calculated using the genet.dist function and represented on heatmaps created with ggplot2 in R (Wickham et al., 2016) for both datasets, single SNPs and only outliers SNPs.

57

Materials and Methods

3.2.3.4.2. Population structure

The genetic relationship and population structure across all individuals were analysed following three approaches with different assumptions. These analyses were performed based on all the variability represented by the first-single SNPs dataset and the variability given by SNPs under selection (outliers).

Discriminant Analysis of Principal Components (DAPC) is a method to obtain a global representation of divergence between populations and consequently know their structure. This method finds genetic clusters in populations and is executed using the R package adegenet (Jombart & Ahmed, 2011). Optimal number of clusters (k) was determined by the lowest Bayesian Information Criterion (BIC) and a successive K‐means algorithm to compare different numbers of groups (k=1 to k=10). In the first step, a Principal Components Analysis (PCA) was applied to examine the differentiation among individuals irrespective of population of origin. All biallelic genome sites coded as 0, 1, or 2, corresponding to homozygosity of the reference allele, heterozygosity, or homozygosity of the alternative allele, respectively. In the second step, individuals were grouped by population for the Discriminant Function Analysis (DFA), which maximises among-group relative to within-group variation. The optimal number of principal components (PCs) retained for the DAPC analysis was assessed determining the highest Mean Successful Assignment Rate and the lowest Root Mean Square Error (RMSE) using the cross-validation (CV) procedure in adegenet with a 10% hold‐out set and 100 replicates (Table 9 and Figure 28).

Admixture v1.3.0 (Alexander et al., 2009) was used to estimate population structure assuming different hypothetical numbers of clusters (k). The programme includes a CV procedure to identify the value of K associated with the best predictive accuracy model. CV was set to 10‐fold to compare different number of K (ranging from 2 to 10) in which lower CV error value indicates the most likely number of clusters.

Neighbor-Joining (NJ) trees were inferred using MEGA7 (Kumar et al., 2016) based on pairwise F_ST distances obtained by both approaches.

   

  validation (CV) procedure to identify the value of k associated with the best predictive accuracy model. CV was set to 10‐fold to compare different number of k (ranging from 2 to 10) in which lower CV error value indicates the most likely number of clusters.

Neighbor‐Joining (NJ) trees were inferred using Mega 7 (Kumar et al., 2016) based on pairwise FST distances.

1.2.3.4.3. Adaptive divergence: Test of selection

The whole SNPs were analysed to find signals of divergent selection among the ten P. lilfordi populations, by means of two different approaches. (i) FST method to detect candidate loci under selection using BayeScan (Foll & Gaggiotti, 2008); and (ii) association analysis to test the correlation between genetic variation and environmental variables through a Redundancy Analysis (RDA) using vegan R package (Oksanen et al., 2018).

PCs Mean Successful Assignment Root Mean Squared Error

All SNPs Outlier SNPs All SNPs Outlier SNPs

5 0.8683 0.7963 0.1460 0.2136

Table 9 . Optimal number of PCs used in DAPC analysis resulting from the cross‐validation method. 

600

Figure 28. Graphics with optimal K number (above) and proportion of variance explained by each eigenvalue (below) for all SNPs (left) and only outliers (right).  

Table 9 . Number of PCs used in DAPC analysis resulting from the cross-validation method. Optimal values are highlighted in bold.

58

Materials and Methods

3.2.3.4.3. Adaptive divergence: Test of selection

The single SNPs dataset was analysed to find signals of divergent selection among the ten P. lilfordi populations, by means of two different approaches: (i) F_ST method to detect candidate loci under selection using BayeScan (Foll & Gaggiotti, 2008); and (ii) association analysis to test the correlation between genetic variation and environmental variables through a Redundancy Analysis (RDA) using vegan R package (Oksanen et al., 2018).

   

  Admixture v1.3.0 (Alexander et al., 2009) was used to estimate population structure assuming different hypothetical numbers of clusters (k). The programme includes a cross-validation (CV) procedure to identify the value of k associated with the best predictive accuracy model. CV was set to 10‐fold to compare different number of k (ranging from 2 to 10) in which lower CV error value indicates the most likely number of clusters.

Neighbor‐Joining (NJ) trees were inferred using Mega 7 (Kumar et al., 2016) based on pairwise FST distances.

1.2.3.4.3. Adaptive divergence: Test of selection

The whole SNPs were analysed to find signals of divergent selection among the ten P. lilfordi populations, by means of two different approaches. (i) FST method to detect candidate loci under selection using BayeScan (Foll & Gaggiotti, 2008); and (ii) association analysis to test

PCs Mean Successful Assignment Root Mean Squared Error

All SNPs Outlier SNPs All SNPs Outlier SNPs

5 0.8683 0.7963 0.1460 0.2136

Table 9 . Optimal number of PCs used in DAPC analysis resulting from the cross‐validation method. 

600

Figure 28. Graphics with optimal k number (above) and proportion of variance explained by each eigenvalue (below) for all SNPs (left) and only outliers (right).

It has been suggested that BayeScan (Foll & Gaggiotti, 2008) is more efficient in recognition of outlier loci with low false positive rate (Pérez-Figueroa et al., 2010). The method implements a rjMCMC approach that can move between a selection model, containing a population-specific component and a locus-specific component, and a model with just a population-locus-specific component (no selection). The posterior probability for selection at a locus is determined by the proportion of MCMC samples that include the model with the locus-specific component. Twenty pilot runs of 5,000 steps were used in the MCMC approach. A burn-in of 50,000 iterations was used, followed by 5,000 iterations using a sampling interval of 10. The prior odds specified for the ratio of neutral: selected sites were set at 100, reflecting that the model with selection might be around 100 times less likely than the model without selection (Lotterhos

& Whitlock, 2014). This prior can have considerable influence on the number of sites detected, so runs were also carried out with 1:10 and 1:1000 proportions. Outliers were identified from the results using the R code supplied with BayeScan and significance was determined under a false discovery rate of 5% to avoid overestimating the percentage of outliers.

59

Materials and Methods

3.2.3.4.4. Correlation with environmental variables (RDA)

Redundancy analysis (RDA) (Forester et al., 2018) using vegan R package (Oksanen et al., 2018) was performed based on the entire genome variability, represented by the single SNPs dataset and using the subset of outlier identified with BayeScan. This analysis was carried out to estimate the proportion of genetic variance that is explained by the history of populations (mtDNA genetic distances), genetic drift (effective population size), and/or divergence selection (ecology variables) (Table 10). Genetic distances based on mtDNA (Terrasa et al., 2009a) alignment (2,382 bp) for all the populations studied were calculated with MEGA7 (Kumar et al., 2016) using Tamura-Nei model with a bootstrap of 1,000.

Resulting distances were transformed into vectors using the function pcoa in the ape package in R with Lingoes correction for negative eigenvalues (Legendre & Anderson, 1999).

In RDA analysis, missing genotypes were imputed on the SNP loci dataset by replacing them with the most common genotype across all individuals. The global significance of the model and the percentage of genetic variance explained by each variable separately were tested with an ANOVA using 1,000 permutations.

Multicollinearity occurs when two or more predictors in a model are correlated providing redundant information. Multicollinearity was measured by variance inflation factor (VIF), which measures the proportion by which the variance of regression coefficient is inflated in the presence of other explanatory variables. VIFs above 20 indicate strong collinearity and above 10 should be examined and avoided if possible. Therefore, variables with VIF > 10 were ruled out (Liu et al., 2017; Borcard et al., 2018).

Table 10. Population name, subspecies designation, presence or absence of melanism, number of samples used in this study (N); mtDNA clade to which they belong according to Terrasa et al. (2009a), geographical localization (latitude and longitude), population density (individual per hectarea (id/ha)) and population size about all popu-lations used in RDA analysis.

Two types of P. lilfordi samples were used in this analysis, tail tissue and eggs at different incubation times. The fieldwork on these protected species was authorised by the Regional Ministry of Agriculture, Environment and Territory (Conselleria d’Agricultura, Medi Ambient i Territori), Government of the Balearic Islands, who specifically authorised capture, release of lizards, and removal of a tail tip for DNA/RNA analysis.

Tail samples were caught by careful noosing following the same method described in section 3.1.1. After sample collection, tissues were immediately stored in RNAlater^® Tissue Protect

Tubes (Qiagen) to preserve RNA integrity and subsequently stored at -20°C when they arrived at the laboratory. Thirty-three tail samples of P. lilfordi populations with different skin colouration (Table 11) were obtained. Ten of them, five individuals from

Foradada and five from Esclatasang islets, both located in the Cabrera archipelago, presented a melanic phenotype characterised by a dark blue abdomen and black dorsal region. The remaining 23 individuals belonged to populations with non-melanic pigmentation and included 11 individuals from Dragonera and 12 from Cabrera harbour (Figure 29). In the Dragonera population, we found individuals with an orange abdomen and brown dorsal while individuals from Cabrera harbour presented a green or blue colouration in the dorsal region and a bright blue abdomen (Figure 30).

Population Subspecies Melanism N mtDNA

clade Latitude Longitude Density (id/ha)

Estimated Table 11. Number of analysed samples (N), phenotype and location.

3.3. MC1R gene expression analyses

In document Genetic analyses of adaptive processes in the podarcis genus (sider 80-83)