*2.1. Phenotypic Evaluation*

Two hundred flax accessions from the Canadian flax core collection were planted in two environments. Evaluation of MC and HC showed a normal distribution across the two environments according to the Shapiro–Wilk normality test and normality plots (Table S1, Figure S1). Variance component analysis indicated significant effects of genotype, environment, and genotype × environment interaction, according to the Wald statistic (*p* < 0.001). The phenotypic variation for MC in Vilcún location 2014 (CAR2014) ranged from 23.52 to 103.57 mg g−<sup>1</sup> with an average of 58.67 mg g−<sup>1</sup> . A lower variation was observed for MC in Huichahue location 2015 (HU2015), which ranged from

18.88 to 91.90 mg g−<sup>1</sup> with an average of 55.04. mg g−<sup>1</sup> HC variation ranged from 35.56 to 48.59% in CAR2014 and from 35.73 to 48.59% in HU2015 (Figure 1, Table S1). MC and HC were significantly positively correlated in CAR2014 and HU2015 with coefficients of 0.28 and 0.29, respectively. Narrow sense heritability (*h* 2 ) for MC attained 70.7 and 73.8% in CAR2014 and HU2015, respectively. Lower *h* <sup>2</sup> of 51.4 and 46.2% for HC at CAR2014 and HU2015 were observed. MC did not differ statistically between flax morphotypes nor seed color classes, according to the Kruskal–Wallis non-parametric test. The average MC was 55.33 and 56.63 mg g−<sup>1</sup> for the fiber and oil morphotypes, respectively (*p* = 0.651) (Figure S2a). The average MC registered values of 56.63 and 59.22 (*p* = 0.517) for the brown and yellow seeded classes, correspondingly. The average HC did not differ statistically between flax morphotypes (fiber = 43.41%, oil = 42.79%; *p* = 0.373). On the other hand, yellow-seeded genotypes averaged 2.66% less HC than brown seeded accessions (*<sup>p</sup>* = 3.2 <sup>×</sup> <sup>10</sup>−<sup>5</sup> ) (Figure S2b). ranged from 35.56 to 48.59% in CAR2014 and from 35.73 to 48.59% in HU2015 (Figure 1, Table S1). MC and HC were significantly positively correlated in CAR2014 and HU2015 with coefficients of 0.28 and 0.29, respectively. Narrow sense heritability (*h*2) for MC attained 70.7 and 73.8% in CAR2014 and HU2015, respectively. Lower *h*2 of 51.4 and 46.2% for HC at CAR2014 and HU2015 were observed. MC did not differ statistically between flax morphotypes nor seed color classes, according to the Kruskal–Wallis non-parametric test. The average MC was 55.33 and 56.63 mg g−1 for the fiber and oil morphotypes, respectively (*p* = 0.651) (Figure S2a). The average MC registered values of 56.63 and 59.22 (*p* = 0.517) for the brown and yellow seeded classes, correspondingly. The average HC did not differ statistically between flax morphotypes (fiber = 43.41%, oil = 42.79%; *p* = 0.373). On the other hand, yellow-seeded genotypes averaged 2.66% less HC than brown seeded accessions (*p* = 3.2 × 10−5) (Figure S2b).

*Int. J. Mol. Sci.* **2018**, *19*, x FOR PEER REVIEW 3 of 16

**Figure 1.** Mucilage (MC) and hull (HC) contents distribution in the association panel in two environments: CAR2014 = Vilcún location 2014, HU2015 = Huichahue location 2015. Values represent the mean of three biological replicates for each trait. **Figure 1.** Mucilage (MC) and hull (HC) contents distribution in the association panel in two environments: CAR2014 = Vilcún location 2014, HU2015 = Huichahue location 2015. Values represent the mean of three biological replicates for each trait.

## *2.2. Population Structure and Linkage Disequilibrium 2.2. Population Structure and Linkage Disequilibrium*

The dendrogram based on 771,914 SNPs and the STRUCTURE plot grouped the 200 individuals into two major clusters arbitrarily assigned as "red" and "blue" (Figure 2a,b). In the *K* against *ΔK* plot, a break in the slope was clearly observed at *K* = 2 (Figure 2b). The red cluster comprised almost exclusively genotypes belonging to the oil morphotype, while the blue cluster included both flax morphotypes. The coefficient of population differentiation (*FST* = 0.08) indicated a weak population The dendrogram based on 771,914 SNPs and the STRUCTURE plot grouped the 200 individuals into two major clusters arbitrarily assigned as "red" and "blue" (Figure 2a,b). In the *K* against ∆*K* plot, a break in the slope was clearly observed at *K* = 2 (Figure 2b). The red cluster comprised almost exclusively genotypes belonging to the oil morphotype, while the blue cluster included both flax morphotypes. The coefficient of population differentiation (*FST* = 0.08) indicated a weak population structure between the two clusters.

structure between the two clusters. The genome-wide linkage disequilibrium (LD) decayed to *r*2 = 0.1 at a physical distance of ~100 kb (Figure 2c). Intrachromosomal LD decayed to *r*2 = 0.1 at a distance between marker pairs ranging from ~40 kb on chromosome 8 to ~210 kb on chromosome 13 (Figure S3). A highly significant positive correlation (*r* = 0.75, *p* = 0.0012) between marker density and the intrachromosomal LD blocks was observed (Figure 3). For example, chromosomes 4 and 8 with the smallest number of markers, and chromosomes 6 and 13 with the largest, displayed the fastest and slowest LD decays, respectively (Figures 3 and S3). The fast LD decays observed in this association panel are indicative of its advantageous potential for reducing QTL intervals and fine mapping of candidate genes for MC and HC. The genome-wide linkage disequilibrium (LD) decayed to *r* <sup>2</sup> = 0.1 at a physical distance of ~100 kb (Figure 2c). Intrachromosomal LD decayed to *r* <sup>2</sup> = 0.1 at a distance between marker pairs ranging from ~40 kb on chromosome 8 to ~210 kb on chromosome 13 (Figure S3). A highly significant positive correlation (*r* = 0.75, *p* = 0.0012) between marker density and the intrachromosomal LD blocks was observed (Figure 3). For example, chromosomes 4 and 8 with the smallest number of markers, and chromosomes 6 and 13 with the largest, displayed the fastest and slowest LD decays, respectively (Figure 3 and Figure S3). The fast LD decays observed in this association panel are indicative of its advantageous potential for reducing QTL intervals and fine mapping of candidate genes for MC and HC.

*Int. J. Mol. Sci.* **2018**, *19*, x FOR PEER REVIEW 4 of 16

*Int. J. Mol. Sci.* **2018**, *19*, x FOR PEER REVIEW 4 of 16

**Figure 2.** Population structure and genome-wide linkage disequilibrium decay. (**a**) Neighbor-joining (NJ) tree for 200 flax accessions based on 779,914 single nucleotide polymorphisms (SNPs); (**b**) Model-based population structure of 200 flax accessions belonging to two clusters predefined by the STRUCTURE software. Each accession is represented by a vertical bar. The color subsections within each vertical bar indicate membership coefficient (Q) to different clusters; (**c**) Genome-wide linkage disequilibrium decay of *r*2 values (red line), against physical distance (kb) using the Hill and Weir (1988) function in *L. usitatissimum*. Yellow line indicates the cutoff value (*r*2 = 0.1) used to determine the genome-wide linkage disequilibrium (LD) block size. **Figure 2.** Population structure and genome-wide linkage disequilibrium decay. (**a**) Neighbor-joining (NJ) tree for 200 flax accessions based on 779,914 single nucleotide polymorphisms (SNPs); (**b**) Model-based population structure of 200 flax accessions belonging to two clusters predefined by the STRUCTURE software. Each accession is represented by a vertical bar. The color subsections within each vertical bar indicate membership coefficient (Q) to different clusters; (**c**) Genome-wide linkage disequilibrium decay of *r* <sup>2</sup> values (red line), against physical distance (kb) using the Hill and Weir (1988) function in *L. usitatissimum*. Yellow line indicates the cutoff value (*r* <sup>2</sup> = 0.1) used to determine the genome-wide linkage disequilibrium (LD) block size. **Figure 2.** Population structure and genome-wide linkage disequilibrium decay. (**a**) Neighbor-joining (NJ) tree for 200 flax accessions based on 779,914 single nucleotide polymorphisms (SNPs); (**b**) Model-based population structure of 200 flax accessions belonging to two clusters predefined by the STRUCTURE software. Each accession is represented by a vertical bar. The color subsections within each vertical bar indicate membership coefficient (Q) to different clusters; (**c**) Genome-wide linkage disequilibrium decay of *r*2 values (red line), against physical distance (kb) using the Hill and Weir (1988) function in *L. usitatissimum*. Yellow line indicates the cutoff value (*r*2 =

0.1) used to determine the genome-wide linkage disequilibrium (LD) block size.

**Figure 3.** Single nucleotide polymorphism (SNP) density plot across the *L. usitatissimum* genome. Numbers of SNPs and LD blocks are also indicated for each of the 15 chromosomes. **Figure 3.** Single nucleotide polymorphism (SNP) density plot across the *L. usitatissimum* genome. Numbers of SNPs and LD blocks are also indicated for each of the 15 chromosomes. **Figure 3.** Single nucleotide polymorphism (SNP) density plot across the *L. usitatissimum* genome. Numbers of SNPs and LD blocks are also indicated for each of the 15 chromosomes.

Three GWA models were tested, including GLM-Q, GLM-PCA, and MLM-K. According to the quantile-qualtile (Q–Q) plot results, the GLM-Q model showed a strong skew toward significance

Three GWA models were tested, including GLM-Q, GLM-PCA, and MLM-K. According to the quantile-qualtile (Q–Q) plot results, the GLM-Q model showed a strong skew toward significance

*2.3. Genome-Wide Association Analysis* 

*2.3. Genome-Wide Association Analysis* 
