From the LC-MS analysis, over 3000 distinct features were detected that represent a mixture of primary metabolites, specialized metabolites, and fragments of abundant proteins. As this initial dataset contains isotopomers and other redundancies, the ramclustR package was applied to reduce the dataset to unique signatures via an examination of the retention time and correlational similarity between all pairs of features in the dataset (Broeckling et al., 2014). After ramclustR was applied, 981 features remained and were used for further analysis. All features were log-transformed to obtain a normal distribution before calculating Best Linear Unbiased Prediction (BLUP) values. BLUPs were calculated from both biological and technical replicates using the Lmer package in R to include all available data in a single estimate for each feature. Biological rep (clone) was used as a random effect while injection and rep were used as fixed effects. Since potato is an autotetraploid, methods that take its ploidy into account are needed. The GWASpoly package was used specifically to account for SNP marker allele dosage (nulliplex, simplex, duplex, triplex, quadruplex). Previous studies involving the SolCAP population have shown that population structure also needs to be accounted for (Rosyara et al., 2016). To control for population structure, several methods were employed. For GWAS, a relationship matrix explicitly defines population structure within the GWASpoly package, so that population structure is addressed during those analyses. For WGCNA, we estimated population structure effects on individual features using the residuals from general linear models. To perform that calculation, principal components (PCA) for each clone were calculated from the genetic data and used to model the BLUP scores using ANOVA, such that BLUP residuals (resBLUPs) were used as inputs for network construction and genetic mapping analyses. Scripts used in this study are reported in Supplemental Materials .

WGCNA was implemented in R through a package developed by Langfelder and Horvath (Langfelder and Horvath, 2008). Data were automatically scaled, and outliers were removed before the start of network construction. Correlation networks were constructed with a threshold of six, based on analysis of this dataset, and minimum module membership was set at five to encourage the detection of small groups of highly correlated features. Potato clone MSL512-6 was removed by the package on account of being an outlier. Module behavior was summarized by WGCNA with an eigenvalue estimated for each group of co-expressed features; eigenvalues were then used for phenotypic correlation analysis and genetic mapping via genome-wide analysis. The WGCNA package labels each module with the name of an arbitrary color, where the grey module refers to those features that were so poorly correlated with others that they are not assigned to another group. The grey module can be thought of as a residual module. Node and edge relationships were visualized using Cytoscape version 3.4.0 (Shannon et al., 2003).

The WGCNA package derives several variables from a data set including module membership and gene trait significance (Langfelder and Horvath, 2008). Module membership refers to the correlation between individual features and the eigenvalue of each module and ranges from -1 to 1 with most values being close to 0. A feature is assigned to the module where its module membership is highest. Gene trait significance is the correlation between an individual feature and the value for each trait in the phenotype dataset, also ranging from -1 to 1 with most values being close to 0. Features with higher gene trait significance scores may be causally related to a trait of interest and thus identifying such features is of value.

We have previously reported SNPs associated with individual features in potato tubers (Levina et al., 2021). In this study, 3521 SNPs were evaluated for association with 44 module eigenvalues using genome-wide association mapping as implemented by GWASpoly (Rosyara et al., 2016). The SolCAP population has previously been genotyped with an Infinium array of 8303 SNP markers; SNPs were not used here if there were more than 20% missing data, or if the SNP was not able to distinguish between the different heterozygous classes (Hirsch et al., 2013). Of the 185 clones summarized in Table 1 , only 184 clones were used for GWAS, since one clone was removed as an outlier during WGCNA. Population structure was taken into account by calculating the K matrix and including it in the model calculations. We performed GWAS with both BLUPs and residual corrected BLUPs. We used the procedure suggested by Lipka et al., 2012 to address population structure. We calculated the first three principal components and then asked what fraction of the signal was explained by principal components as opposed to markers. If the population structure was highly associated with a specific phenotype then we removed the effect of that population structure. A genome-wide significant threshold was determined using Bonferroni correction (0.05 divided by the number of markers used) which addressed the problem of testing multiple hypotheses and allowed us to focus on modules and features with very strong correlations.

Previously, we used UPLC-MS/MS to conduct metabolic profiling on a subset of the SolCAP potato diversity panel to study tuber composition (Levina et al., 2021). We measured 981 small molecules and peptide fragments (or “features”), where we obtained structural information on 6.8% of them and used GWAS to map genetic determinants for 48%. From that analysis, we identified a large number of mappable features, detected four regions of the genome associated with unexpectedly high numbers of features, and highlighted new candidate loci for glycoalkaloids (Levina et al., 2021).

?>In a preliminary analysis, using data uncorrected for population structure, we found that WGCNA assembled a network composed of 42 modules, each containing between five and 220 members, with only five features assigned to the grey (residual) module (data not shown). Because of the potentially confounding aspect of population structure on co-expression networks, we also calculated the network after correcting individual feature BLUP scores for population structure using genetic marker information as analyzed by principal components analysis. In this improved network, we observed 44 modules with five to 171 members each, similar to the initial network ( Supplemental Table 1 ). However, the grey module expanded from 5 to 73 members, which suggests that an appreciable fraction of spurious associations exist in the naïve analysis and that addressing population structure explicitly is appropriate. The following correlation and genetic analyses utilize the population structure-adjusted co-expression network.

Once the relationships between all members of an expression data set were calculated, WGCNA was used to conduct a correlation analysis between traits and module eigenvalues ( Figure 1 ). For this study, we used phenotype data collected by SolCAP ( Supplemental Table 2 ) (Douches, 2008; Rosyara et al., 2016). All traits were significantly correlated with at least one module, with the magnitudes varying from r = |0.18| for tuber length and free sucrose to r = |0.46| for chip color ( Figure 1 ).

Modules that exhibited a significant positive correlation (r>0.2, p value< 0.05) value with the fry quality traits chip color, log(glucose), and log(fructose) were tan, brown, and salmon, while modules that negatively correlated with fry quality were thistle1, midnightblue, royalblue, brown4, and lightcyan ( Figure 1 ). We also observed correlations between modules and other traits collected by SolCAP: yield was correlated with yellow and sienna3, tuber size was correlated with mediumpurple3, midnight blue and blue, while vine maturity was correlated with darkgrey ( Figure 1 ).

To identify genetic regions that influence tuber composition module (rather than an individual feature) behavior, GWAS was performed on module eigenvalues ( Supplemental Table 3 ) using GWASpoly, an R package that allows tetraploid allele dosage to be considered. Twenty-two out of 44 modules were significantly associated with at least one SNP marker ( Table 2 ).

Supplemental Table 4 lists all modules significantly associated with one or more SNPs, the genetic model used to make each association, the identity of significant SNPs, and their chromosome positions. SNPs associated with modules were found on chromosomes 3, 5, 6, 7, 8, and 12, with the preponderance of associations found on chromosomes 3, 7, and 8 – a similar distribution to that seen for SNPs associated with individual features in (Levina et al., 2021). All of the QQ plots and Manhattan plots are collated in Supplemental Figure 1 . Of the modules with genetic associations, three were associated (r>0.2) with chip color: thistle1, midnightblue, and tan.

Selection for chip color in the Cornell breeding program is currently based solely on phenotypic evaluation as we have no useful molecular markers to identify superior or inferior clones. Accordingly, we sought to better understand connections between chip score, potato composition modules, and specific features within those modules to create useful molecular markers (genetic, biochemical, or both) that can facilitate our selection efforts. To find such markers, we evaluated module membership (MM) and gene-trait significance (GS) for modules that exhibited the strongest correlation with chipping scores. MM and GS are among the derived variables created by the WGCNA package, where both range from -1 to 1 and most values are near 0. MM is the correlation between individual features and the module eigenvalue of each module, while GS is the correlation between each feature and each phenotype. The plot of module membership vs. gene trait significance of chip color yields several modules with substantial correlations (r² >0.2) ( Figure 2 ).

Specifically, brown (r^{2 =} 0.81) and tan (r^{2 =} 0.26) modules show a very strong positive correlation while midnightblue (r^{2 =} 0.36) and thistle1 (r^{2 =} 0.38) show a strong negative correlation ( Figure 2 ). Note that lower chip score values indicate superior performance, thus features with the largest negative GS values indicate predictors of superior chipping ( Supplemental Table 5 ). Many features with extreme positive and negative GS scores have a sufficiently simple inheritance to be mapped using GWAS. Of the features with the 25 most extreme module membership values for each of brown, tan, midnight blue, and thistle1, the locations of those that could be genetically mapped are summarized in Table 3 . The highest value of MM within each module was used to select those features because this allowed the selection of features that also had the highest absolute value of GS of chip color for each module. The map locations of those features were not evenly distributed across the genome but were instead clustered on chromosomes 3, 6, and 8 ( Table 3 ). The list of all features and their correlation with chip color is presented in Supplemental Table 6 .

As reported by (Levina et al., 2021), loci associated with individual features of tuber composition are not evenly distributed throughout the potato genome; some regions are associated with a large number of features. With network analysis, it became clear that these genetic hotspots represent the detection of many members of a small number of modules ( Figure 3 ).

For example, SNP c2_20259 was associated with the largest number of individual features (108, using the additive genetic model) according to GWAS. Network analysis revealed that these 108 features are members of just seven modules (of the 44 detected; Figure 3 ). While the turquoise module is the largest in the network with 170 features, only one feature from this module (C327) is associated with SNP c2_20259. The black module is the second largest, with 71 members, but only five of its features are correlated with this SNP. Neither the turquoise nor black module eigenvalues were associated with SNP c2_20259 by GWAS. The remaining 106 features associated with SNP c2_20259 are all members of modules that are also associated with c2_20259: blue (where 45 of 59 members are associated with this SNP), dark orange (32/34), midnight blue (16/28), medium purple (4/7) and violet (3/10); Table 3 ; Figure 3 ). The representation of features associated with c2_20259 in these modules is far higher than expected by chance, with a nearly 3x enrichment for the violet module and 5x or more for the other four (chi-square test, p ≤ 0.0001). To simplify the presentation in Figure 3 , only those correlations (edges) between features (nodes) with r² > 0.45 are shown. At this threshold, not all nodes are connected to other members of their own module. For the dark orange module, five of the nodes are polypeptides ( Figure 3 ), one is a fatty acid, and one is a polyphenol. One member of the midnight blue module has been identified as a polypeptide.

The amount of structural information about individual features varies among modules, with only half of the modules having any information available. In terms of absolute number, the brown module (4^th largest) has the most information available, for 15 of 43 members. In terms of density, the bisque4 module (15^th largest) is best with 12 of 18 members annotated. While WGCNA estimates connectivity between nodes (features) via the topological overlap measure (TOM), here we use r² to report a correlation between nodes, using the population structure-corrected residuals of BLUPs. All members of the bisque4 module are highly correlated with one another, with a median r² of 0.43. We focused on feature C433 in the bisque4 module that has previously been identified as β-chaconine (Levina et al., 2021), a glycoalkaloid commonly found in potato tubers. Glycoalkaloids impart a bitter taste and breeders want to minimize their levels. Unexpectedly, as shown in Figure 4 , β-chaconine and its immediate biosynthetic precursor α-chaconine were not in the same module. β-chaconine and two other glycoalkaloids (feature C38, a tetrahydropentoxyline-like compound, and feature C176, α-chaconine) are highlighted with yellow ovals in Figure 4 ). α-chaconine is highly (r² = 0.97) correlated with the unknown feature C170 while essentially uncorrelated with β-chaconine (r² <0.01. β-chaconine is highly correlated with C38, C223 (a terpene), C578 (a polypeptide likely formed by hydrolysis of an abundant protein), and C915 (no structural information). β-chaconine is moderately (0.5<r²<0.75) correlated with eight other members of the bisque4 module, including C83 (a fatty acid), C150 (a polypeptide), C609 (a polyphenolic compound), C780 (choline) and C797 (a polypeptide). β-chaconine is modestly (0.2<r²<0.5) correlated with other members of the bisque4 module, including C627 (a napthofuran) and C795 (a polypeptide).