{'fraction_non_alphanumeric': 0.030047604489671136, 'fraction_numerical': 0.026903087740752063, 'mean_word_length': 4.904589994842702, 'pattern_counts': {'":': 0, '<': 8, '<?xml version=': 0, '>': 7, 'https://': 2, 'lorem ipsum': 0, 'www.': 0, 'xml': 0}, 'pii_count': 0, 'substrings_counts': 0, 'word_list_counts': {'cursed_substrings.json': 0, 'profanity_word_list.json': 0, 'sexual_word_list.json': 0, 'zh_pornsignals.json': 0}}
Natural variation at FLOWERING LOCUS T2 mediates local adaptation in a key life history trait in European aspen

The timing of growth and dormancy represent critical life-history trade-offs in perennial plants and often show strong local adaptation. Despite their importance, the genetic architecture of phenological traits remains poorly understood. Here we identify a ~700 Kbp region mediating local adaptation in the timing of bud set in Populus tremula across a latitudinal gradient (~56-66° N) in Sweden, where extensive gene flow has removed almost all traces of population structure. The strongest genomic associations with bud set are centered on a P. tremula homolog of FLOWERING LOCUS T2 (PtFT2) and explain 32-52% of the variation in bud set across sites and years. This region also shows multiple signs of a recent selective sweep restricted to the northernmost populations. Under field and greenhouse conditions, variation in bud set among accessions originating from different latitudes is strongly associated with variation in PtFT2 gene expression. Finally, transgenic down-regulation of PtFT2 yields a phenotype that, under field conditions, closely mimics variation observed between phenotypic extremes in natural populations. Our results thus provide evidence of one genomic region in P. tremula, centered on PtFT2, with major effects on adaptation to shorter growing seasons and colder climates following post-glacial colonization.

the obligately outcrossing, woody perennial plant European aspen (Populus tremula). We 48 performed whole-genome re-sequencing on 94 unrelated P. tremula trees that display a strong 49 latitudinal cline in bud set (Fig. 1A,B; Supplementary Table 1) [2][3][4] . After stringent variant 50 calling and filtering (see "Online methods"), we identified a total of 4,425,109 high-quality 51 single nucleotide polymorphisms (SNPs) with a minor allele frequency (MAF) greater than 52 5%. We found very weak population structure across the entire range using principal 53 component analysis (PCA) 5 , with a single significant axis separating individuals according to 54 latitude (r=0.889, P-value <0.001) but explaining only 1.3% of the total genetic variance (Fig. 55 1C; Supplementary Table 2). Consistent with this, a Mantel test also showed a weak pattern 56 of isolation by distance (r=0.210; P-value=0.047, Supplementary Fig. 1). Although recent 57 admixture of divergent postglacial lineages following the Last Glacial Maximum (LGM) is 58 capable of generating a genome-wide pattern of clinal variation 6  We used three complementary approaches to identify candidate SNPs involved in local 73 adaptation (See "Online Methods"). We identified SNPs that were most strongly associated 74 with observed population structure using PCAdapt 7 , SNPs showing strong associations with 75 environmental variables based on a latent factor mixed-effect model (LFMM) 8  performed genome-wide association mapping (GWAS) on the timing of bud set, our target 77 adaptive trait, using GEMMA ( Fig. 2A, 9 ). SNPs identified as significant (false discovery 78 rate<0.05) by the three methods showed a large degree of overlap ( Supplementary Fig. 3) and 79 for subsequent analyses we consider SNPs that were identified as significant by at least two of 80 the three methods to be involved in local adaptation. 99.2% of the 910 SNPs identified by all 81 three methods and 89.1% of the additional 705 SNPs identified by two methods were located 82 in a single region spanning c. 700 Kbp on chromosome 10 ( Fig. 2A PtFT2 is known to be involved in controlling seasonal phenology in perennial plants 10 and 91 has previously been implicated in regulating short-day induced growth cessation, bud set and 92 dormancy induction in Populus [11][12][13] . Although both copies of PtFT2 appear to be expressed 93 In order to further understand the evolution of functional differences between northern 118 and southern PtFT2 alleles, we analysed sequence diversity in the region surrounding PtFT2 119 and identified multiple signs of a strong and recent selective sweep that has been restricted to 120 the northern-most populations. We observed a strong reduction in genetic diversity in the 121 northern populations and high genetic differentiation between southern and northern 122 populations around the PtFT2 locus (   Table 5). Two haplotype homozygosity-based tests, iHS 15 and nSL 16 , also reveal a high 130 concentration of significant SNPs with signals of positive selection within this region 131 ( Supplementary Fig. 10). Using the most strongly associated SNP (Potra001246:25256) as a 132 proxy we observed that haplotypes carrying the derived allele have elevated extended 133 haplotype homozygosity (EHH, 17 ) relative to the ancestral allele (Fig. 2H,I) and the high 134 EHH allele is largely restricted to high-latitude populations (Fig. 3D). We assumed that the 7 allele present in northern populations has been driven to fixation by a hard selective sweep 136 and used an Approximate Bayesian Computation (ABC) method 18  When considered together with the very weak population structure in P. tremula across 184 Sweden, our findings provide empirical evidence supporting recent theoretical predictions 185 that local adaptation in the face of high levels of gene flow tends to favor a genomic 186 architecture of one or a few tightly clustered loci of large-effect as opposed to the more 187 traditional architecture of many loci of small effect expected for most complex traits 19 . This 188 is likely because a single large-effect or a group of tightly linked loci of smaller effects are 189 better able to resist the homogenizing effect of gene flow than independent small-effect loci 190 do 20 . While small-effect loci likely also contribute to local adaptation, variation at such loci is 191 expected to be transient and/or difficult to distinguish from neutral background variation, 192 especially when recent demographic history coincides with local adaptation, as is the case in 193 the present study ( Supplementary Fig. S12) 21,22 . 194 In summary, we identified a recent and strong selective sweep at PtFT2, a P. tremula 195 homolog of FT, which is regionally restricted to the northernmost populations and likely 196 occurred in response to the substantially shorter growing seasons that are characteristic of 197 northern Scandinavia. While FT is mainly involved in controlling flowering in plants, 198 functional diversification of FT has occurred repeatedly in many plants 23   Online Materials and Methods

Sample collection and sequencing
We collected material from all available trees in the Swedish Aspen (SwAsp), which consists of 116 individuals collected from 12 different locations spanning the distribution range in Sweden 1 (Fig. 1A). Leaf material was sampled from one clonal replicate of each individual growing at a common garden experiment located in Sävar, (iv) sites that were in genomic scaffolds with a length shorter than 2 Kbp.

SNP and genotype calling
SNP calling in each sample was performed using the GATK HaplotypeCaller, and GenotypeGVCFs were then used to perform the multi-sample joint aggregation, regenotyping and re-annotation of the newly merged records among all samples. We performed several filtering steps to minimize SNP calling bias and to retain only highquality SNPs: (i) Remove SNPs at sites not passing all previous filtering criteria; (ii) Retain only bi-allelic SNPs with a distance of more than 5 bp away from any indels; (iii) Remove SNPs for which the available information derived from <70% of the sampled individuals after treating genotypes with quality score (GQ) lower than 10 as missing; (iv) Remove SNPs with an excess of heterozygotes and deviates from

Relatedness, population structure and isolation-by-distance
To identify closely related individuals and to infer population structure among the sampled individuals, we discarded SNPs with missing rate >10%, MAF < 5% and that failed the Hardy-Weinberg equilibrium test (P <1×10 -6 ) after all filtering steps as shown above. We also generated LD-trimmed SNP sets by removing one SNP from each pair of SNPs when the correlation coefficients (r 2 ) between SNPs exceed 0.2 in blocks of 50 SNPs using PLINK v1.9 9 . This yielded 217,489 independent SNPs that were retained for downstream analyses of population structure. First, we used PLINK

Screening for SNPs associated with local adaptation
We used three conceptually different approaches to test for genome-wide signatures of local adaptation.

Population differentiation
First, we detected candidate SNPs involved in local adaptation using the principal component analysis as implemented in PCAdapt 14 . PCAdapt examines the correlations (measured as the squared loadings ρ 2 jk , which is the squared correlation between the jth SNP and the kth principal component) between genetic variants and specific principal components (PCs) without any prior definition of populations. As only the first principal component was significant from the PC analysis (Supplementary Table 2), we only estimated the squared loadings ρ 2 j1 with PC1 to identify SNPs involved in local adaptation. Our results showed that most outlier SNPs that were highly correlated with the first population structure PC also had high F ST values between populations (Supplementary Fig. 15). Assuming a chi-square distribution for the squared loadings ρ 2 j1 , as suggested by 14 , we used PCAdapt to compute P-values for all SNPs, and then calculated the false discovery rate (FDR) using the method of 15 to generate a list of candidate SNPs showing significant associations to population structure. Only SNPs with FDR < 5% were retained as those significantly involved in local adaptation.

Environmental associations
We tested for the presence of candidate SNPs that exhibited high correlations with environmental gradients. To do this, a total of 39 environmental variables were analysed. Precipitation and temperature values were retrieved from WorldClim version 1 16 Supplementary Fig.16A), we performed a PCA on these variables using the 'prcomp' function in R to identify PCs that best summarized the range of environmental variation. The first environmental PC, which explained > 60% of the total variance ( Supplementary   Fig.16B,C) and had the strongest loadings for the length of growing season ( Supplementary Fig.16D), was kept to represent our target environmental variable for further analyses. We then used a latent factor mixed-effect model (LFMM) implemented in the package LEA in R 17 to investigate associations between SNPs and the first environmental PC while simultaneously accounting for population structure by introducing unobserved latent factors into the model 18 . Due to the weak population structure found in the SwAsp collection (see Results), we ran the LEA function lfmm with the number of latent factors (K) ranging from one to three, using 5000 iterations as burn-in followed by 10,000 iterations to compute LFMM parameters for all SNPs. This was performed five times for each value of K, and we observed identical results both across different values of K and across independent runs within each value of K (data not shown). We only showed the results using K=2 to account for the background population structure. LFMM outliers were detected as those SNPs with FDR < 0.05 after using the method of 15 to account for multiple testing.

Genome-wide association mapping
We obtained previously published measurements of the timing of bud set, which is a highly heritable trait that shows strong adaptive differentiation along the latitudinal gradient 19,20 . To measure phenotypic traits, all SwAsp individuals have previously been clonally replicated (four ramets per individual) and planted at two common garden sites in 2004 (Sävar, 63ºN, and Ekebo, 56ºN) (Fig. 1A). The common garden setup is described in detail in 1  been shown to outperform methods that try to correct for population structure by including it as a fixed effect in the GWAS analyses 23 . Given the extremely weak population structure we observe in our GWAS population (see Results) we did not pursue any further corrections for population structure in the association analyses as this likely would severely reduce our power to detect significant associations. As described previously, we used a FDR < 5% 15 to control for the multiple testing across the 4,425,109 SNPs. To test for the accuracy of imputation, and its relationship with the MAF cutoff and the missing rate of genotypes in our dataset, we selected 346,821 SNPs with a rate of missing genotypes lower than 10% from the pseudo-chromosome 2 (~32.6

Testing for positive selection on the region surrounding the
Mb) for the simulation analysis. We randomly masked out varying proportions (5-50%) of SNPs, which were treated as missing. BEAGLE v 4.1 was then used to impute genotypes at the masked positions. We found high imputation accuracy (>0.97) across a wide range of MAF when rates of missing genotypes were less than 30% ( Supplementary Fig.17), suggesting imputation and phasing by BEAGLE should not bias the accuracy of our results. We therefore phased and imputed genotypes of the SNPs anchored on pseudo-chromosomes using BEAGLE v 4.1.

Estimation of ancestral states for all SNPs
Since the ancestral states of SNPs are usually used for selection detection, for each SNP, we classified alleles as either ancestral or derived on the basis of comparisons with two outgroup species: P. tremuloides and P. trichocarpa. We obtained publicly available short read Illumina data for one P. tremuloides (SRA ID: SRR2749867) and one P. trichocarpa (SRA ID: SRR1571343) individual from the NCBI Sequence Read Archive (SRA) 26 . We individually aligned the reads from these two samples to the de novo P. tremula assembly (Potra v1.1, available at PopGenIE.org) and used UnifiedGenotyper in GATK to call SNPs at all sites (--output_mode EMIT_ALL_SITES). For each SNP, two procedures were performed to define their ancestral states: (1) because P. trichocarpa is more distantly related to P.
tremula compared to P. tremuloides 27 and from our previous study there were less than 1% polymorphic sites shared between P. tremula and P. trichocarpa 26 , we inferred the ancestral state as the P. trichocarpa allele at sites where the P.
trichocarpa individual was homozygous and matched one of the P. tremula alleles; Otherwise, (2) we inferred the ancestral state as the P. tremuloides allele at sites where the P. tremuloides individual was homozygous and matched one of the P. tremula alleles. If the above two requirements were not met, the ancestral state was defined as missing. In total, we obtained information of ancestral states for 96.3% of all SNPs.

Refine the anchoring and orientation of SNPs associated with local adaptation to a single region on chromosome 10
As we found that a large majority of significant SNPs (>90%) detected by at least two of the three methods (PCAdapt, LFMM, and GEMMA) were clustered in a single genomic region on pseudo-chromosome 10, we performed several further steps to refine the anchoring and orientation of these SNPs within this region. First, we used ColaAlignSatsuma from the Satsuma packages 25 to align the genomes of P. tremula and P. trichocarpa using default settings. The output was then converted and filtered into GBrowse synteny compatible format that was available at http://popgenie.org 3 . Based on the alignment of the two genomes, 15 scaffolds from the P. tremula assembly that contain SNPs inferred to be associated with local adaptation were completely or partially mapped to a single region on chromosome 10 of P. trichocarpa genome (Supplementary Table 3). We then retained only seven scaffolds that were completely mapped to the region and with length longer than 10 Kbp. The seven scaffolds contained more than 95% (1465 out of 1528) of the total number of significant SNPs in the single region of chromosome 10. Lastly, according to the alignment results between the genome of P. tremula and P. trichocarpa, we reordered and re-oriented the seven scaffolds to a ~700 Kbp region for all downstream selection tests (Supplementary Fig. 4).

Linkage disequilibrium
A possible explanation for the occurrence of the clustering of SNPs involved in local adaptation in the ~700 Kbp region is the presence of genomic rearrangements, such as chromosomal inversions, that suppress recombination and create extensive linkage disequilibrium (LD) spanning the region within the rearrangement 28 . We calculated correlations (r 2 ) between all pairwise common SNPs (MAF>5%, 9149 SNPs) in the ~700 Kbp region using PLINK 1.9 9 . To explore and compare patterns of LD between the ~700 Kbp region on chromosome 10 and genome-wide levels, we first used PLINK 1.9 to randomly thin the number of common SNPs across the genome to 200,000, and then calculated the squared correlation coefficients (r 2 ) between all pairs of SNPs that were within a distance of 100 Kbp. The decay of LD against physical distance was estimated using nonlinear regression of pairwise r 2 vs.
the physical distance between sites in base pairs 29 . In contrast to the expectation of the presence of an inversion, we did not observe blocks of elevated LD across the region on chromosome 10 ( Supplementary Fig. 18A,B). Although LD is slightly higher in this region compared to genome-wide averages, LD still decays rapidly and falls to background levels within a few thousand bases ( Supplementary Fig. 18A).
This suggests that frequent recombination has occurred in this region and that the clustering of SNPs involved in local adaptation more likely arose from selective sweeps 30 . Nonetheless, owing to the limited ability to detect inversions using shortinsert paired reads, future characterization of structural variation across the genome is clearly required to determine whether, and how, genomic rearrangements are involved in mediating adaptation in Populus.

Positive selection detection
We measured two haplotype-based tests, integrated haplotype score (iHS) 31  To test for whether there is significant concentration of selection signals on the region surrounding the PtFT2, we divided the 19 pseudo-chromosomes (without the seven scaffolds around the PtFT2 locus) into non-overlapping windows of 700 Kbp and calculated the proportion of SNPs with |iHS| > 2 or with [nS L |>2 in each window.
Statistical significance was assessed using the ranking of genome-wide windows, with windows having fewer than 100 SNPs being excluded.

Population-specific selective sweeps
Several standard methods were further applied to search for signs of selective SweepFinder2 is most efficient when information on the ancestral and derived states is available for SNPs and we therefore polarized SNPs as described in Section 5.2. The small fraction of SNPs (~3.7%) that could not be polarized were excluded from further analysis using SweepFinder2. CLRs were calculated using non-overlapping windows with a spacing of 2 Kbp, and the empirical site frequency spectrum across the whole P. tremula genome was estimated using the -f option in SweepFinder2 after including all polymorphic sites in the genome (a total of 8,007,303 SNPs). As recommended by 39 we only used sites that were polymorphic or that represented fixed substitutions in each group of populations to scan for sweeps.
To determine whether there are significant differences of the above statistics between the 700 Kbp region around PtFT2 gene on chromosome 10 and genome-wide estimates, we use the same strategy to divide the genome into the windows with the same size for each test and calculated the above statistics across the genome (Results are shown in Supplementary Fig. S7 and Table S5). Significance for the above statistical measurements was evaluated using Mann-Whitney tests.
To assess the scale of a genomic region that is affected by a selective sweep, we ran coalescent simulations modeling a selective sweep in the Northern populations.
Simulations were run assuming that the selected site was located at the centre of the simulated region. Parameters for the simulations were taken from ABC calculations dating the selective sweep inferred in the North populations (See Section 5.7). Briefly, we used a scaled population mutation rate (4N e µ) of 0.0081/bp, which corresponds to the average observed diversity in the North populations. Similarly we set the scaled population recombination rate (4N e r) to 0.0019 to match the genome-wide ratio of r/µ=0.229 in P. tremula 26 . Analyses of the simulated data using SweepFinder2 showed that a single selective sweep often yields multiple significant peaks across a region spanning up to, and even exceeding, 100 Kbp (95% quartile: 148221 bp; Supplementary Fig. 19).

Dating the selective sweep in the North populations
To date the inferred selective sweep in the North populations we used the Approximate Bayesian Computation (ABC) method described in 40  replicates to serve as biological repeats, flash-frozen in in liquid nitrogen and stored at -80°C until sample preparation. Samples were collected at 2h intervals for a total period of 24 h. RNA extraction for all samples was performed using a CTAB-LiCl method 43 . CDNA synthesis was performed using the iScript cDNA Synthesis Kit (BIO-RAD) according to the manufacturer's instructions. Quantitative real-time PCR analyses were performed using a Roche LightCycler 480 II instrument, and the measurements were obtained using the relative quantification method 44 . We used primers qFT2F (5'-AGCCCAAGGCCTACAGCAGGAA-3') and qFT2R (5'-GGGAATCTTTCTCTCATGAT-3') for amplifying the transcript of FT2 and qUBQF (5'-GTTGATTTTTGCTGGGAAGC-3') and qUBQR (5'-GATCTTGGCCTTCACGTTGT-3') for UBQ as the internal control.

Field experiment with transgenic PtFT2 lines
Construction of the PtFT2 RNAi lines are described in detail in 45 Fig. 20) and classified as active growth (score 3), growth cessation (score 2), bud formation (score 1) and bud set (score 0).
Scoring was continued until all plants had completely senesced in late October. Bud scoring data was converted to Julian date of bud set and analysed using the following linear model: where µ is an overall mean, α i is the effect of treatment i (where i is either PtFT2 RNAi or WT) and β j is the effect of block j and ε ij are individual residual errors.

Sequence and SNP data
The raw sequencing reads have been deposited in NCBI's short-read archive,

Other data
The following data is available from zendo.org (https://doi. org/10.5281/zenodo.844372) • Bud set genetic values (BLUPs) for all clones used in the GWAS

BASH, Perl, Python and R scripts
All scripts used for the analysis described are available in an online repository at https://github.com/parkingvarsson/PhotoperiodLocalAdaptation