The plant materials used in this project include 385 accessions that were collected from more than 35 flax-growing countries. Approximately 60% of the germplasm originated from the Old-World, i.e., where flax has been cultivated for millennia, including the postulated centers of origin and secondary centers of diversity.

Genome-wide SNP datasets were extracted from a 1.7M SNP dataset originally generated by resequencing the flax core collection (n = 407) using the Illumina HiSeq 2000 platform [31, 33]. SNPs were filtered with a 90% call rate and a minimum minor allele frequency (MAF) of 5%. To reduce the number of missing SNPs, imputation was performed using LinkImpute [34] implemented in TASSEL v5 [35] with default parameters, but with the maximum distance between sites set to 100 kb. Individuals with $>$10% missing SNPs after imputation were omitted. The resulting SNP dataset is referred to as the 277K dataset. A second dataset was prepared by extracting from the 277K dataset only the SNPs without missing data and a MAF $>$10%. The resulting dataset is later defined as the 23K dataset. A third dataset containing only SNPs located in exons was filtered from the original SNP dataset using an 80% call rate and a MAF $>$5%. This dataset was also imputed. The resulting exon dataset was further filtered to retain only SNPs with a 90% call rate after imputation. This set is referred to as the 34K exon dataset.

In order to define the appropriate number of ancestral populations, estimates were obtained using the pcadapt tool [36]. Here eigenvalues were computed for 50 PCs and visualized into scree plots. By applying the cattle rule [37], the PC to the left of the last PC with eigenvalues that deviates from the smooth line is considered to be the most appropriate number of populations. To corroborate the estimated number of populations, a cross validation [38] was performed using ADMIXTURE [39] which uses a Bayesian clustering approach. PCA was performed based on the estimated number of populations (K) and a neighbor-joining (NJ) phylogenetic analysis was carried out using TASSEL [35]. The PCA based on the first three PCs and the NJ tree were visualized using ggplot2 [40] in R and interactive Tree Of Life (iTOL) [41], respectively.

Following the PCA clustering using the 23K dataset, individuals were assigned to one of 12 populations. To estimate the genetic variation between populations, pairwise genetic differentiation between populations was estimated with the fixation index (F${}_{S\mathit{}T}$) based on population size estimated after 10000 permutations at $\alpha$ = 0.05 and 0.01 using Arlequin 3.5 [42]. To quantify the contribution of each SNP to the variation, the F${}_{S\mathit{}T}$ of individual SNPs was obtained for each dataset using the R package LEA which performs Landscape and Ecological Association (LEA) analyses [43]. Haplotypes were assessed using the web-based tool SNiPlay v3 [44] based on SNPs with significant (p 0.05/n, where n = number of SNPs in the data set) F${}_{S\mathit{}T}$ values at $\alpha$ = 0.05 from the 23K dataset for each chromosome. Haplotypes were defined as those present in at least three individuals within the overall germplasm collection or one of its populations. To determine the distribution of private haplotypes for each population, haplotypes observed in at least three individuals were considered for each chromosome. To understand the effect of bottlenecks in haplotype, the per population gene diversity was calculated for each of the 12 populations using GENEPOP v4.7.5 [45] in R package genepop [46].

The geographic region/country that best represents the population was assigned for each population based on the passport data, which indicates the origin of dominant members of a population. The representative geographic areas were inferred based on the history of cultivation of flax in each of the regions (Supplementary Table 1). Therefore, geographic coordinates of a representative district or a province with a long history of flax cultivation were used as environmental data tags. To get insight into some environmental factors, mean annual temperature data downloaded from www.wroldclim.org was overlaid onto the world map using DIVA-GIS v7.5 (LizardTech Inc, Portland, OR, USA) [47]. The representative coordinates of each population were positioned onto the world map along with the annual temperature data (Fig. 1).

Fig. 1.

World map overlaid with mean annual temperatures and showing the geographic locations assigned to the populations (Pop1-12).

The monthly mean temperature and day-length data for each population were downloaded from 22 years of records available in the NASA database (https://power.larc.nasa.gov/). The average records of these periods were considered. The cropping period for the selected regions was determined based on the major crop calendars from FAO (http://www.amis-outlook.org/amis-about/calendars/soybeancal/en/) for all populations that follow similar northerly cropping patterns. Because flax is a Rabi crop in some parts of the world such as India and Pakistan, different cropping calendars (https://nfsm.gov.in/nfmis/rpt/calenderreport.aspx and http://namc.pmd.gov.pk/crop-calender.php) were used for genotypes of Indian and Pakistani origins, respectively.

To detect outlier loci associated with local adaptation signatures, two genome-scan methods were used: individual SNP F${}_{S\mathit{}T}$ using LEA [43] and the principal component-based to detect local adaptation (PCADAPT) as implemented in the R package pcadapt [36]. Each dataset was analyzed separately using both methods and results were compared to identify robust genetic signature loci. To capture adaptive loci associated with specific environmental conditions, a GEA analysis was performed following the latent factor mixed model 2 (LFMM2) using lfmm in R [48] based on a PC of average cropping season temperature, day-length and latitude of the representative locations of each population for the 23K dataset. The use of a representative PC was chosen because the environmental factors were highly correlated (Supplementary Fig. 1). For each population and all three datasets, allelic frequencies were calculated and summarized for all loci detected in all the datasets by at least two methods.

Genomic regions that consistently showed a strong association with an environmental factor across the datasets for both of the first two methods (PCADAPT, and F${}_{S\mathit{}T}$) were retained as potential signatures of local adaptation. For the GEA, genome association with the first PC, which explained more than 97% of the variation of the three environmental factors of latitude, cropping season average temperature and day-length among the populations, was applied (Supplementary Fig. 2).

Genomic regions spanning 20 kb up and downstream of the most significant SNP loci were examined for the presence of candidate genes for local adaptation using the flax reference genome annotation [33]. Linkage disequilibrium (LD) between candidate genes and their associated marker was calculated using gpart package in R [49]. Candidate gene putative functions were further assessed through the identification of their Arabidopsis orthologs (www.arabidopsis.org) and through a literature search evidencing their role(s) in adaptation.

Filtering of the datasets yielded a total of 277,399, 23,592 and 34,451 SNPs for the two genome-wide (277 and 23K) and the exon-based (34K) datasets, respectively. The former two contained 385 genotypes, while the exon-based dataset included 393. The estimated number of populations [36] was 12 (Supplementary Fig. 3) which is also consistent with the result from cross validation technique where the lowest error was obtained at K = 13 (Supplementary Fig. 4). The populations tended to follow geographic gradient where Pop1-5 were dominated by Eurasian accessions (Supplementary Table 1. Pop6, Pop10 and Pop12 contained mostly Canadian, Abyssinian and Mediterranean accessions, respectively, whereas the majority of the South Asian accessions grouped into the remaining three populations (Fig. 2A,B). The NJ phylogenetic analysis clustered the accessions slightly differently, i.e., reflecting both geographic and historical-use patterns of variation (Fig. 2C). Accessions from Old-World flax-growing regions tended to have longer branches compared to those from the New-World regions. In addition, the majority of the fiber types clustered in the single clade Pop1_FIB (Fig. 2C).

Fig. 2.

Population genetic structure based on the 23K SNP dataset. (A) and (B) principal component analysis (PCA) clustering of the accessions based on the first three principal components (PC1, PC2 and PC3). (C) Neighbor-joining (NJ) phylogenetic tree where accessions from the Old-World flax growing regions tended to have longer branches. Accession names in NJ indicate the type (O, oil; F, fiber; U, unknown), the country of origin, the breeding status (C, cultivar; B, breeding material; L, landrace) followed by the accession name.

Pairwise population differentiation was significant (p 0.01) for all comparisons except between one of the temperate populations (Pop4) and a population dominated by Canadian cultivars (Pop6). Most populations dominated by Eurasian and Canadian accessions displayed low differentiations (Table 1). Populations with the strongest (F${}_{S\mathit{}T}$ = 0.76) differentiation were the fiber-dominated Pop1 and the South-Asian Pop11 (Table 1).

Pop1 and Pop2 harbored the most common haplotype from each of the 15 chromosomes (Table 2). All accessions considered, chromosome 1 displayed the lowest percentage of common haplotype with ~19% while chromosome 4 had the highest with 78%. Some populations, such as Pop8, 9 and 11, contained very few common haplotypes while others, such as Pop1, 2, 4 and 6, comprised a large proportion of individuals with common haplotypes across all chromosomes (Table 2).

Private haplotypes are those found in a single population. Based on this definition, 30 private haplotypes were observed in seven of the populations (Supplementary Table 2). The South Asian (Pop11) and Abyssinian (Pop10) populations contained 48% and 20% of the private haplotypes, respectively (Fig. 3). The most frequently observed private haplotype was Chr14: Hap1 which was present in 97% of the accessions of the South Asian population (Pop11) representing 8.6% of all accessions (Supplementary Table 2). The highest gene diversity was in Pop12 followed by Pop5 which putatively originated from Mediterranean Portugal and Turkey regions respectively (Supplementary Table 1). In contrast, Pop7 and Pop10 displayed no diversity (Table 2).

Fig. 3.

Frequency distribution of private haplotypes by population. Only seven of the 12 populations had private haplotypes.

A total of 19 outliers were detected, including six on chromosomes 1, 8, 9 and 12 that were detected in all three datasets with both PCADAPT and F${}_{S\mathit{}T}$ methods (Fig. 4 and Supplementary Table 3). The allele frequencies at outlier loci detected in all the three datasets differed among populations and where the frequencies were higher in populations from old flax growing regions for either of the two alleles at a locus (Fig. 5). Candidate gene searches revealed several genes with putative function in ecological adaptation. For example, the Chr9:9310330 locus harbors gene Lus10006147, an ortholog to ArabidopsisAT4G15530, which encodes a pyruvate orthophosphate dikinase. The locus marked by Chr1:5466653 harbors the predicted flax gene Lus10011967, whose ortholog AT4G18130 encodes a phytochrome E (PHYE). The other locus on chromosome1 (Chr1:10413007) includes two genes: Lus10022627 and Lus10022628, which are orthologs of AT2G36800 and AT2G36780, encoding URIDINE 5’-DIPHOSPH (UDP)-GLYCOSYLTRANSFERASE 73C5 (UGT73C5) and UDP-glycosyltransferase 73C3 (UGT73C3), respectively. The Chr8:2932993 locus contained the tandemly repeated genes Lus10012356, 7 and 8, orthologous to the Arabidopsis gene AT3G07870 that encodes the F-BOX PROTEIN92 (FBX92). The Chr7:10407650 locus, one of the most significant SNPs of the 277K (p~3.48E-52 for PCADAPT; p 8.1 $\times$ 10${}^{-224}$ for F${}_{S\mathit{}T}$) and 23 K (p~6. 7 $\times$ 10${}^{-28}$ for PCADAPT; p~1.94 $\times$ 10${}^{-66}$ for F${}_{S\mathit{}T}$) datasets contains the predicted flax gene Lus10000371, which is orthologous to the early flowering gene AT1G17455 encoding ELF4-L4.

Fig. 4.

Manhattan plots (left panels) and quantile-quantile (Q-Q) plots (right panels) of genome-wide associations using PCADAPT (A,C,E) and F${}_{\mathrm{S}\mathit{}\mathrm{T}}$ (B,D,F) of the 277K (A,B), 23K (C,D) and exon-34K (E,F) datasets. Outliers with the same number indicate the same SNP. The numbers before and after the hyphen indicate the chromosome number and the outlier SNP as listed in Supplementary Table 3, respectively.

Fig. 5.

Distribution of alleles at loci detected in the three SNP datasets among the populations. Populations from old flax growing region (Pop3, 5, 7–12) displayed a higher frequency of one of the alleles at each locus than the remaining populations from relatively newer flax growing regions.

PC-based genome scan of the 34K dataset detected several significant SNPs on chromosome 1 between positions 3612570 and 3671655 (Fig. 4) with p-values between 2.1 $\times$ 10${}^{-16}$ and 6.9 $\times$ 10${}^{-12}$. This locus contains six genes of which five were orthologous to AT2G30140 encoding UDP-GLUCOSYL TRANSFERASE 87A2 (UGT87A2), and one was orthologous to AT2G30150 encoding UDP-GLUCOSYL TRANSFERASE 87A1 (UGT87A1). All UGT genes at this locus were in strong LD with their associated SNPs (Supplementary Fig. 5).

GEA based on PC1 associated with latitude, average day-length and temperature during the cropping season performed using the assigned coordinates of the population resulted in significant SNPs on chromosomes 3, 7, 9 and 13 (Fig. 6). The significant locus at Chr3:16799360 harbored multiple candidate genes for flowering time regulation, including the AT2G29950 ortholog Lus10040667 that is predicted to be an EARLY FLOWERING LOCUS-LIKE1 (ELF4-L1) gene (Table 3). The locus Chr7:16799360 contains AT1G67170 ortholog Lus10015495 that is predicted to encode a FLOWERING LOCUS C EXPRESSOR-LIKE 2 (FLL2) gene (Table 3).

Fig. 6.

Manhattan and quantile-quantile (Q-Q) plots of the first principal component (PC) for latitude, day-length and temperature of the cropping seasons using the genome-wide 23K SNP dataset showing outliers on four of the 15 flax chromosomes.

The genetic structure observed in this experiment is consistent with what was establish in previous works, which demonstrated the clustering of accessions attributed to their eco-geographic origin [31, 54]. The relatively high differentiation and high concentration of private haplotypes in populations dominated by South Asian accessions may suggest the adaptation of the crop to warm regions [55, 56]. This is unlike the cool and/or temperate habitats where flax is widely grown [57]. Flax in this warm South Asian region might also be adapted to short-day photoperiod given that flax accessions from this region differ from the common flax adapted to grow under the long-day seasons of most cooler regions [56, 58]. In a similar way, the higher rate of private haplotypes in the Abyssinian population may be attributable to adaptation of the crop to the equatorial region, with equinox effects, elevation-induced cool temperate climate and windstorms [31, 59]. Despite the high private haplotype concentration, the low gene diversity in the Abyssinian population may suggest effect of genetic drift [60]. The lack of private haplotypes in some populations such as Pop1, 2, 4, and 6, dominated by Eurasian and Canadian accessions, can be due to introduction of materials of different origins, selection of dual flax types (both fiber and oil) or through hybridization with wider germplasm in breeding programs [31, 61].

The outlier loci harbored genes of known roles for local adaptation. Most loci detected using the PCADAPT and F${}_{S\mathit{}T}$ methods are linked to genes that mediate responses to stresses. The flax genes Lus10011967, Lus10022627 and Lus10022628 that are predicted to encode UGT73Cs might be involved in Fusarium oxysporum wilt tolerance in flax [62, 63, 64] and other plant species [65, 66, 67]. Dmitriev et al. [68] demonstrated the up regulation of UGT73C3 in flax in response to Fusarium infection. By and large, UGTs play an important role in Fusarium wilt resistance [69] including in flax [62, 64]. Fusarium fungi are the most common diseases in many crop plants and can cause devastating losses. Fusarium wilt in flax is one of the most severe biotic stresses and may lead to a complete loss of flax production [70]. Hence, Fusarium diseases can be one of the natural selection forces that result in genetic divergence in many crops [71], including flax [72]. The consistent outlier SNP Chr1:10413007, that marks the locus harboring both UGT73C3 and UGT73C5, is likely an important genetic signature for local adaptation [73, 74]. The other locus marked by Chr1:5466653 included a gene predicted to encode PHYE. This gene was hypothesized to be involved in regulating responses to light quality and temperature [75]. As such, this locus may be an essential genetic signature of divergent adaptation to different eco-geographic regions.

The ELF4-L4 gene at Chr7:10407650 locus plays a significant role in circadian clocking [76] and flowering time [77], which are both affected by day-length, and consequently, latitude. The tandemly repeated FBX92 genes at Chr8:2932993 locus can also be important signatures of adaptation of flax to varying environmental factors along its latitudinal gradient. FBX92 mediates responses to different abiotic stresses including light [78]. The FBX92 protein affects leaf sizes in Arabidopsis [79], which may contribute to adaptation to latitude-induced abiotic stresses such as temperature [80]. The multiple copies of the UGT87A gene at locus Chr1:3612570 and Chr1:3671655 might also play an important role in flowering time regulation and abiotic stress response. In Arabidopsis, UGT87A2 mutants overexpressing the flowering repressor FLOWERING LOCUS C (FLC) had delayed flowering times regardless of the duration of the day-length [81]. UGT87A2 is also involved in plant adaptation to osmotic stresses, including drought and salinity, via regulation of multiple genes mediating responses to these stresses [82]. The Chr1:3612570-3671655 locus harbors multiple copies of the UGT87A2 gene which may suggest its involvement in flax adaptation through regulation of stress responses [81, 82]. Other outlier loci also harbored genes of importance for local adaptation. For example, Lus10042995 and Lus10042996 were predicted to be orthologous to ethylene response factors ERF106 (AT5G07580) and ERF105 (AT5G51190), respectively. The former is DECREASE WAX BIOSYNTHESIS2 (DEWAX2) that negatively regulates cuticular biosynthesis [83], and, as such, adversely impacts cuticular wax-related tolerance to abiotic stresses such as drought [84]. In contrast, the latter (ERF105) is a transcription factor involved in freezing tolerance [85]. Moisture availability and temperature are among the major factors that shaped the genetic structure of flax, and they constitute major determinants of the success of the crop in its current eco-geographic regions [31], which span the warm semi-arid Indian subcontinent [55, 86] to the relatively cool and humid temperate Eurasian regions [87].

Most of the outlier SNP loci that are associated with the PC1 of major environmental factors such as cropping season day-length, temperature and latitude, also harbored important genes with known roles in adaptive divergence of plants along eco-geographic gradients. The ELF4-L1 orthologous gene at Chr3:6169770 locus is predicted to regulate circadian rhythm and flowering time [88, 89], suggesting the adaptive signature of this locus in flax. Chr7:16799360, which marks the locus harboring the predicted FLOWERING LOCUS C EXPRESSOR-LIKE 2 (FLL2) gene, is another crucial latitudinal gradient adaptation signature. FLLs have been reported for their role in the regulation of FLOWERING LOCUS C (FLC) and vernalisation [90, 91]. The AT2G02540 orthologous genes Lus10010693 and Lus10010694 at the Chr13:10377859 locus also mediate flowering time via positive regulation of FLC [92, 93]. Flax is a “long-day” plant whose flowering time can be determined by both photoperiod and vernalisation [94]. The strong association between loci harboring genes that regulate flowering time and vernalisation suggests the importance of these loci for genetic divergence, which allowed flax to expand to vast geographic regions. Since GEA was performed based on representative putative locations, validation using materials of known and precise geographic coordinates and associated site factors such as climatic, edaphic and biodiversity records, remains warranted.