Skip to content

Advertisement

  • Research article
  • Open Access

Linkage analysis and QTL mapping in a tetraploid russet mapping population of potato

BMC Genetics201819:87

https://doi.org/10.1186/s12863-018-0672-1

  • Received: 9 March 2018
  • Accepted: 12 September 2018
  • Published:

Abstract

Background

Genome-wide single nucleotide polymorphism (SNP) markers coupled with allele dosage information has emerged as a powerful tool for studying complex traits in cultivated autotetraploid potato (Solanum tuberosum L., 2n = 4× = 48). To date, this approach has been effectively applied to the identification of quantitative trait loci (QTLs) underlying highly heritable traits such as disease resistance, but largely unexplored for traits with complex patterns of inheritance.

Results

In this study, an F1 tetraploid russet mapping population (162 individuals) was evaluated for multiple quantitative traits over two years and two locations to identify QTLs associated with tuber sugar concentration, processing quality, vine maturity, and other high-value agronomic traits. We report the linkage maps for the 12 potato chromosomes and the QTL location with corresponding genetic models and candidate SNPs explaining the highest phenotypic variation for tuber quality and maturity related traits. Significant QTLs for tuber glucose concentration and tuber fry color were detected on chromosomes 4, 5, 6, 10, and 11. Collectively, these QTLs explained between 24 and 46% of the total phenotypic variation for tuber glucose and fry color, respectively. The QTL on chromosome 10 was associated with apoplastic invertases, with ‘Premier Russet’ contributing the favorable allele for fry processing quality. On chromosome 5, minor-effect QTLs for tuber glucose concentration and fry color co-localized with various major-effect QTLs, including vine maturity, growth habit, tuber shape, early blight (Altenaria tenuis), and Verticillium wilt (Verticillium spp.).

Conclusions

Linkage analysis and QTL mapping in a russet mapping population (A05141) using SNP dosage information successfully identified favorable alleles and candidate SNPs for resistance to the accumulation of tuber reducing sugars. These novel markers have a high potential for the improvement of tuber processing quality. Moreover, the discovery of different genetic models for traits with overlapping QTLs at the maturity locus clearly suggests an independent genetic control.

Keywords

  • Potato
  • Solanum tuberosum
  • Single nucleotide polymorphism
  • Quantitative trait locus
  • Vine maturity
  • Tuber processing quality

Background

With a complex polyploid genetics the cultivated autotetraploid potato (Solanum tuberosum L., 2n = 4× = 48) presents unique challenges for linkage mapping and quantitative trait loci (QTL) analysis. The recent development of potato single nucleotide polymorphism (SNP) arrays provides a set of genome-wide markers that can be used in mapping and quantitative trait locus (QTL) analysis [14]. Combined with the availability of statistical tools that enable the utilization of allele dosage information (i.e., the number of copies of each allele at a given polymorphic locus), it is now possible to significantly increase the power of SNP information to detect recombination events and to incorporate this information into QTL mapping [57]. To date, this approach has been successfully used to detect novel QTLs and candidate SNPs associated with late blight resistance [4, 8], internal heat necrosis [9], Potato virus Y [10], and agronomic traits [11], but largely unexplored for studying processing quality and other high-value traits with complex inheritance.

There is a need for molecular markers linked to tuber processing quality traits to assist the selection of germplasm for cultivar development. Processing quality is largely influenced by tuber sugar content. The accumulation of fructose and glucose (reducing sugars) in potato tubers during cold storage negatively impacts the end-use quality characteristics. Reducing sugar content is a complex trait controlled by multiple genetic and environmental factors [12, 13]. The complex genetic architecture of this trait is consistent with the activity of numerous enzymes involved in the metabolic pathways linking starch breakdown and synthesis to sugar formation and utilization [14, 15]. Genetic mapping studies have shown that fry processing quality and reducing sugars are associated with a large number of small effect QTLs [1620]. Soltis-Kalina et al. [21] identified two major QTLs for chip color on chromosomes 1 and 6 in a diploid mapping population. More recently, Braun et al. [22] reported two QTLs for resistance to cold-induced sweetening on chromosomes 4 and 6 of diploid potato. Most of these studies have been performed on diploid potato germplasm. However, further studies are needed to better understand the genetic architecture of these traits at the tetraploid level in populations derived from russet-skinned potato varieties.

The objective of this research was to identify QTLs in a tetraploid F1 mapping population derived from a cross between ‘Rio Grande Russet’ × ‘Premier Russet’; two varieties representative of the U.S. market class characterized by long tubers with russet skin suitable for use in processing or fresh consumption. The population was evaluated for multiple agronomic and processing quality traits in two contrasting environments over 2 years. Using SNP markers and the recently developed TetraploidSNPMap software [6] we were able to generate high-resolution linkage maps and provide the QTL location, with corresponding genetic models and candidate SNPs for potato tuber processing quality and vine maturity related traits in adapted potato germplasm utilized by the commercial sector.

Results and discussion

SNP genotyping

A total of 7157 SNPs from the Infinium 8303 Potato Array map to unique positions in the Potato Genome Sequence Consortium (PGSC) version 4.03 Pseudomolecules of the reference potato genome [4]. After initial filtering to remove SNPs with missing theta values and monomorphic markers, 3739 segregating SNPs were available for further analyses. Additional filtering steps removed SNPs with double reduction, segregation distortion, and missing values. The resulting 2479 SNPs were used to run TetraploidSNPMap. Duplicate and near-duplicate SNPs (i.e., those that differed for at most 2/162 offspring from another SNP) were further removed during the twopoint ordering routine. A final set of 1545 SNPs was used to construct the linkage maps (Table 1).
Table 1

Number of SNPs and length of linkage groups in the two parental genetic maps of ‘Rio Grande Russet’ (RGR) and ‘Premier Russet’ (PR)

Chr

No. Mapped SNPs

Map Length (cM)

Map Length (Mb)

PGSC v4.03 PM (Mb)

Coverage

Average Inter-loci Distance (cM)

 

Total

RGR

PR

RGR

PR

RGR

PR

DM

RGR

PR

RGR

PR

1

174

146

131

89.0

89.0

88.3

88.3

88.7

99.5

99.5

0.61

0.68

2

158

126

126

73.5

75.9

36.3

36.7

48.6

74.7

75.5

0.58

0.60

3

125

97

106

76.5

76.5

61.4

59.7

62.3

98.6

95.8

0.79

0.72

4

143

122

125

83.6

80.1

71.0

71.0

72.2

98.3

98.3

0.69

0.64

5

110

89

140

56.8

58.2

49.7

49.7

52.1

95.4

95.4

0.64

0.42

6

132

106

115

65.5

63.8

59.4

58.9

59.5

99.8

99.0

0.62

0.55

7

91

68

126

65.8

64.4

50.2

49.8

56.8

88.4

87.7

0.97

0.51

8

129

110

105

74.4

77.1

56.8

55.2

56.9

99.8

96.9

0.68

0.73

9

132

107

114

83.4

85.8

58.3

61.1

61.5

94.8

99.3

0.78

0.75

10

109

81

101

70.2

70.3

59.3

59.2

59.8

99.2

99.0

0.87

0.70

11

151

123

108

72.0

70.7

44.7

44.7

45.5

98.2

98.2

0.59

0.65

12

91

76

128

77.4

77.4

58.8

58.8

61.2

96.1

96.1

1.02

0.60

Total

1545

1251

1425

888.1

889.1

694.2

693.1

725.1

95.2

95.1

0.71

0.62

Chr chromosome, cM centimorgan. Map length (Mb) and map coverage values are based on the PGSC version 4.03 Pseudomolecules of the reference potato Solanum tuberosum group Phureja DM1–3516 R44 (DM)

The most informative parental genotype classes, simplex (AAAB × AAAA, ABBB × BBBB), duplex (AABB × AAAA, AABB × BBBB) and double-simplex (AAAB × AAAB, ABBB × ABBB), accounted for more than 50% of the 1545 SNP loci. Among higher dosages, simplex-duplex configurations were the most represented (Additional file 1: Table S1). Overall, the parental genotypes shared alleles at a large number of SNP loci. This is consistent with previous results obtained with the Infinium 8303 Potato Array suggesting that french-fry processing and table russets are members of a single russet-type germplasm [23]. From the perspective of genetic mapping, the incorporation of shared alleles enables the alignment of the parental maps so that allele effects from both parents can be studied simultaneously [5].

Linkage map construction

Of the 1545 SNPs used in the construction of linkage maps, 1251 were heterozygous in ‘Rio Grande Russet’ and 1294 in ‘Premier Russet’, one thousand of which were heterozygous in both parents. Constructed linkage maps spanned a genetic distance of 888 and 889 cM, for ‘Premier Russet’ and ‘Rio Grande’, respectively. Map length of individual chromosomes ranged from 64.42 to 89.03 cM, with an average of 106 SNP markers per chromosome and a marker density of 1.4 SNP per cM. On average, both maps covered 95% of the PGSC v4.03 Pseudomolecules (Table 1, Additional file 2: Figure S1). The overall maps (Chr 1–12) and corresponding homologs (H1-H8) are listed in Additional file 3: Table S2.

Phenotypic trait data

Considering locations and years individually, a total of 23 traits were mapped (Fig. 1). For any given trait, the abbreviated name is followed by location, Idaho (ID) and North Carolina (NC) and year (2010, 2011). Traits included, bud-end fry color (FBE), stem-end fry color (FSE), citric acid (CA), early blight (EB), growth habit (GH), tuber glucose (GLU), tuber shape (TSHA), vine maturity at 95 days after planting (MAT95), vine maturity after 120 days after planting (MAT120), and Verticillium wilt (VW). When the mean from 2 years of observations was used, only the abbreviated name for the trait and location is indicated (e.g., MAT95ID). All traits exhibited a broad range of phenotypic variation (Additional file 4: Table S3; Additional file 5: Table S4). The broad-sense heritability (H2) for both years at each location was statistically significant for all traits, with H2 estimates ranging from 0.45 for GLUNC to 0.93 for MAT120ID (Table 2). For traits with higher heritability, the means from 2 years were used to perform the QTL analysis (Table 3).
Fig. 1
Fig. 1

Correlations between traits and locations, based on least square means. EB, early blight; CA, citric acid; FBE, bud-end fry color; FSE, stem-end fry color; GH, growth habit; GLU, glucose; MAT95, vine maturity at 95DAP; MAT120, vine maturity at 120DAP; TSHA, tuber shape; VW, Verticillium wilt; ID, Idaho; NC, North Carolina. Cells with no significant correlations (P-value < 0.01) are left blank. Pearson’s correlation coefficients are indicated in the cells

Table 2

Broad-sense heritability and correlation for seven phenotypic traits for potatoes grown in Idaho (ID) and North Carolina (NC)

Trait

Location

Heritability

Correlation

Vine maturity 95

ID

0.90

0.82

Vine maturity 95

NC

0.83

0.75

Vine maturity 120

ID

0.93

0.86

Vine maturity 120

NC

0.86

0.77

Growth habit

ID

0.85

0.73

Growth habit

NC

0.57

0.41

Fry bud end

ID

0.83

0.71

Fry stem end

ID

0.74

0.61

Glucose

ID

0.74

0.73

Glucose

NC

0.45

0.30

Tuber shape

ID

0.77

0.62

Table 3

Quantitative trait loci (QTLs) for the phenotypic traits evaluated on the A05141 population during 2010 and 2011

Trait

Location (year)

Chr

Position (cM)

LOD

R 2

95% PT

Model; Homologous Chromosome; and Candidate SNPs1

Maturity 95 DAP

Idaho (2010–11)

5

15.0

25.2

44.3

3.4–3.9

Double-simplex; H2, H8; c2_11605

Maturity 120 DAP

Idaho (2010–11)

5

15.0

34.1

54.5

3.4–3.9

Double-simplex; H2, H8; c2_11605

Maturity 95 DAP

North Carolina (2010–11)

5

15.0

24.0

42.7

3.4–3.9

Double-simplex; H2, H8; c2_11605

Maturity 120 DAP

North Carolina (2010–11)

5

15.0

21.2

38.5

3.2–3.7

Double-simplex; H2, H8; c2_11605

Growth habit

Idaho (2010–11)

5

16.0

17.9

32.7

3.4–3.9

Double-simplex; H2, H8; c2_11605

Growth habit

North Carolina (2010–11)

5

15.0

8.0

14.8

3.2–3.7

Double-simplex; H2, H8; c2_11605

Tuber shape

Idaho (2010–11)

5

9.0

4.3

7.5

3.2–3.6

Double-simplex; H2, H8; c2_11605

Early Blight

Idaho (2011)

5

15.0

20.6

41.3

3.3–3.6

Double-simplex; H2, H8; c2_11605

Verticillium wilt

Idaho (2011)

5

15.0

24.3

46.5

3.3–3.8

Double-simplex; H2, H8; c2_11605

Citric acid

North Carolina (2010)

5

9.0

5.3

10.7

3.2–3.7

Double-simplex; H2, H8; c2_11605

Glucose

Idaho (2011)

5

21.0

3.9

6.2

3.1–3.5

Double-simplex; H2, H5

Fry bud end

Idaho (2010–11)

5

14.0

4.7

8.3

3.1–3.5

Duplex, H5, H7; c2_50301

Glucose

Idaho (2010–11)

4

66.0

6.0

11.3

3.7–4.1

Double-simplex; H2, H5; c2_12976

Fry bud end

Idaho (2010–11)

4

66.0

5.2

9.4

3.6–4.1

Double-simplex; H2, H5; c2_12976

Fry stem end

Idaho (2010–11)

4

66.0

4.5

7.2

3.6–4.0

Double-simplex; H2, H5; c2_12976

Fry bud end

Idaho (2010–11)

6

43.0

5.4

10.8

3.6–4.0

Double-simplex; H2, H5; c2_56145

Fry stem end

Idaho (2010–11)

6

43.0

7.6

16.1

3.5–4.0

Double-simplex; H2, H5; c2_56145

Glucose

Idaho (2010–11)

6

53.0

5.7

11.3

3.3–3.8

Double-simplex; H2, H5; c2_56145

Glucose

Idaho (2010–11)

10

55.0

4.5

7.5

3.5–4.1

Simplex; H5; c2_48118

Fry bud end

Idaho (2010–11)

10

55.0

4.3

6.8

3.4–3.9

Simplex; H5; c2_48118

Fry bud end

Idaho (2010–11)

11

39.0

5.3

10.5

3.1–3.4

Double-simplex; H4, H5

1Candidate SNPs correspond to the Infinium 8303 Potato Array designations (solcap_snp)

QTL position is indicated by chromosome (Chr) and map position in centimorgan (cM). LOD logarithm of the odds; R2, explained variance; PT, 95% permutation threshold

QTL analysis

Vine maturity

‘Premier Russet’ has a late vine maturity [24] and ‘Rio Grande Russet’ has a medium late maturity. In the present study, both parents scored as expected, with ‘Rio Grande Russet’ showing slightly earlier maturity scores (5.0–6.5) than ‘Premier Russet’ (6.0–7.5). For the progeny, mean trait values displayed a wide range of phenotypic variation (Additional file 4: Table S3), with a high correlation between years (> 0.70, P < 0.001) and across environments (> 0.55, P < 0.001). The broad-sense heritability for both years at each location was statistically significant, with H2 estimates ranging from 0.83 to 0.93 (Table 2). These results were consistent for all vine maturity traits. When mean values for maturity were mapped to the linkage maps, a significant association was detected on chromosome 5 (Fig. 2, Additional file 6: Figure S2). This association was consistent between years and environments. Logarithm of the odds (LOD) scores ranged from of 21.2 and 34.1, well above the permutation thresholds. Collectively, they explained between 38 and 54% of the trait variance (Table 3). The peak of the LOD profile was consistently detected at 15.0 cM (Table 3, Fig. 2). The best simple model for maturity was a double-simplex configuration (AAAB × AAAB) on homologous chromosomes H2 and H8 of Rio Grande russet (RGR) and Premier Russet (PR), respectively (Additional file 6: Figure S2). This model implies a QTL with allele effects from both parents, and with the B allele associated with earlier maturity (P < 0.001). Thus, progeny with two copies of the B allele (AABB) had on average earlier maturity (3.1–4.0) than progeny with a single copy (AAAB, 4.5–6.2) or lacking the B allele (AAAA, 6.4–7.5).
Fig. 2
Fig. 2

LOD score distribution of QTLs for traits at the maturity locus region on chromosome 5. Each panel in the figure represents a different trait. LOD threshold 95% confidence intervals are listed in Table 3. CA, citric acid 2010; EB, early blight; FBE, bud-end fry color Idaho; GH, growth habit; GLU, glucose; M95, vine maturity at 95DAP; M120, vine maturity at 120DAP; TSHA, tuber shape; VW, Verticillium wilt; ID, Idaho; NC, North Carolina; DAP, days after planting; numbers 10 and 11 in the trait names indicate years 2010 and 2011, respectively

The SNP closest to the QTL for maturity that had a double-simplex configuration was solcap_snp_c2_11605, located at 9.1 cM. Regression of the maturity mean values on the genotype of this SNP explained between 38 and 50% of the phenotypic variance. The peak of the LOD curve (15 cM) was positioned within a gap between markers that spanned 8 to 10 cM, with the solcap_snp_c2_11605 located on one side of the gap. This region harbors the Solanum tuberosum cycling DOF factor gene 1 (StCDF1, PGSC0003DMG400018408), which was reported to be a major regulator of potato plant maturity [25]. Graphical representation of the genetic and physical (Mb) distance of homologous chromosomes H2 and H8 displayed the expected curve shape and fit well with chromosome structure indicating a high degree of concordance between the genetic position (cM) and physical location (Mb) (Additional file 7: Figure S3).

Multiple traits co-localize at the maturity locus region on chromosome 5

Eight traits including citric acid (CANC10), growth habit (GHID, GHNC), bud-end fry color (FBEID), early blight score (EBID11), tuber glucose (GLUID11), tuber shape (TSHAID), and Verticillium wilt score (VWID11) co-localized with vine maturity traits (MAT95ID, MAT95NC, MAT120ID, MAT120NC) on chromosome 5 (Table 3, Fig. 2). The best genetic model for all traits, except GLUID11 and FBEID was a double-simplex effect on H2 and H8 with the SNP, solcap_snp_c2_11605, having this configuration (Table 3). When the residuals from the regression of all six traits on vine maturity were analyzed, no significant effect of a QTL was found for EBID11, TSHAID, and GHNC. For growth habit (GHID) and VWID11, a significant QTL was observed only when corrected for maturity at 95DAP but was not significant at 120DAP. Analysis of glucose concentrations (GLUID11) showed a LOD of 3.9 at 21 cM, with a QTL effect on H2 and H5. After adjusting for maturity, a peak in the LOD curve for GLUID11 was observed just below the 95% permutation threshold.

For FBEID, a duplex model (BBBB × AABB) with additive effects of the A allele on H5 and H7 fit the data significantly better than a model with a double-simplex with additive effects of the A allele (Additional File 6: Figure S2). In our analysis, this model (H5, H7) corresponds to the low temperature sweetening resistant parent ‘Premier Russet’, with the A allele associated with an increase in light reflectance. As reflectance readings are inversely proportional to fry color and sugar content, this implies that progeny with two A alleles (AABB) had on average a higher reflectance measurement, and therefore a lighter bud-end fry color than those with a single or no A allele (Fig. 3a). A candidate SNP for the effect on H5 and H7 is the solcap_snp_c2_50301, which maps to the distal end of the short arm of chromosome 5 at position 5.1 Mb on the methyltransferase gene (PGSC0003DMG400031262) of the reference potato genome (Solanum tuberosum group Phureja DM1–3516 R44) PGSC v4.03 Pseudomolecules. This region harbors two genes involved in starch/sugar inter-conversion, alpha amylase (AMY-5, PGSC0003DMG401017626) and sucrose transporter (Sut2, PGSC0003DMG400025610). Both are located within 1 Mb of the SNP solcap_snp_c2_50301. Analysis of the residual FBE after regression on maturity (95DAP, 120DAP) indicated a peak of the LOD curve that was just below the 95% permutation threshold. When FBEID was analyzed per year, the residuals for FBE 2011 increased the LOD above the upper permutation threshold. Collectively, these results indicate that QTLs for EBID11, VWID11, TSHAID, and GHID could be explained by the effect of maturity. However, results for FBEID and GLUID11 suggest an independent genetic control.
Fig. 3
Fig. 3

Box plots showing the effects of the marker genotypic classes (AAAA, AAAB, AABB, ABBB, BBBB) of 4 candidate SNPs on tuber glucose concentration (GLU); bud-end fry color (FBE); stem-tuber fry color (FSE); ID, Idaho

The co-localization of multiple QTLs at the maturity locus region in this study was consistent with the observed patterns of phenotypic correlations. Vine maturity was significantly correlated (P-value < 0.01) with several traits overlapping the same genetic position (Fig. 1). For instance, EBID11 and VWID11 scores both showed a strong negative correlation with MAT95ID and MAT120ID. This trend may explain the major effect QTL with later maturity closely associated with delayed onset of VW and EB, or alternatively, early maturity leading to early onset of VW and EB. A relatively high correlation was also observed between vine maturity and growth habit, with early maturity associated with more spreading growth habit. Conversely, a low correlation (r ≤ 0.2) between vine maturity and both tuber glucose concentration and bud-end fry color was consistent with different models for the QTL effects.

Tuber glucose, bud-end and stem-end fry color

In addition to chromosome 5, significant QTLs for glucose, bud-end and stem-end fry color were detected in overlapping regions of chromosomes 4, 6, 10, and 11. On chromosome 4, all three QTLs were at 66.0 cM with the LOD scores ranging from 4.5 to 6.0 (Table 3, Fig. 4). These QTLs explained between 7.2 and 11.3% of the trait variance. Analysis of different simple models indicated a single best-fitting model for all three traits, with a double-simplex allele (AAAB × AAAB) on homologous chromosomes H2 and H5 of RGR and PR, respectively (Table 3, Additional file 8: Figure S4). For glucose, the B allele was associated with higher glucose concentration, while for bud-end and stem-end fry color, the presence of the B allele indicated lower reflectance values (Fig. 3b).
Fig. 4
Fig. 4

LOD score distribution of QTLs for three traits across the 12 potato chromosomes (Chr). Abbreviations are as in Fig. 1

A candidate SNP for a double-simplex configuration on H2 and H5 is solcap_c2_12976 at position 67.4 cM. On the reference PGSC v4.03 Pseudomolecules, this SNP maps to the long arm of chromosome 4 at position 69.1 Mb on the F-box and WD40 domain protein gene (PGSC0003DMG400003339), which is close (< 2 Mb) to several genes associated with starch/sugar inter-conversion as annotated in the potato reference genome, and/or associated with tuber quality traits [26]. They include alpha amylase (AMY-4/1 PGSC0003DMG400007974, AMY-4/2 PGSC0003DMG400009891), hexose transporter (HT-4/3, PGSC0003DMG400009994), invertase (NV-4, PGSC0003DMG400009936), and starch branching enzyme I (SBE I, PGSC0003DMG400009981).

On chromosome 6, significant QTLs mapped at 53 cM for glucose and at 43 cM for FBE and SBE, respectively. The LOD profiles ranged from 5.4 to 7.6 and explained between 10.8 and 16.1% of the trait variance, respectively (Table 3, Additional file 9: Figure S5). The best model for all three traits was a double-simplex configuration with allele effects on H2 and H5. As in chromosome 4, the presence of the B alleles was associated with an increase in glucose concentration and a decrease in reflectance readings (Fig. 3c). There are four candidate SNPs with a double simplex configuration nearby (solcap_snp_c2_33777, solcap_snp_c1_10130, solcap_snp_c2_5769, solcap_snp_c2_56145, c1_3003). All were significantly (P-value < 0.0001) associated with the mean trait values, however, solcap_snp_c2_56145 had on average the highest percent of variance explained (9% to 16%). The SNP solcap_snp_c2_56145 maps at position 51.4 Mb of the potato reference genome on the ATP binding protein gene PGSC0003DMG400033089. This SNP is located in a region (~ 7 Mb) that contains several genes associated with starch/sugar inter-conversion including hexokinase (HXK-6, PGSC0003DMG400016521), fructokinase (FRK-6/1, PGSC0003DMG400027017, FRK-6/2, PGSC0003DMG400026916), and invertase (INV-na-6/1, PGSC0003DMG400026107, INV-6/2, PGSC0003DMG400033142) genes.

On chromosome 10, a QTL for GLU and FBE at 55.0 cM was detected (Table 3, Fig. 4). The best genetic model was for the low temperature sweetening resistant parent ‘Premier Russet’, having a simplex allele AAAB, with the B allele associated with an increase in reflectance readings and a decrease in glucose concentration (Fig. 3d). The closest SNP with this configuration was solcap_snp_c2_48118 (55.1 cM), which map to the beta-fructofuranosidase (invertase) gene (InvCD111, PGSC0003DMG401028252,) at position 55.8 Mb. InvCD111 and the adjacent InvCD141 gene (PGSC0003DMG402028252) are tandem-duplicated genes of the Invap-a locus, which encode cell wall bound (apoplastic) invertases [27, 28]. Apoplastic as well as vacuolar invertases have been associated with tuber quality traits [29, 30].

On chromosome 11, a QTL for FBEID at position 39.0 cM accounted for 10.5% of the trait variation. A dominant double-simplex model was the best fit for the QTL located on H4 and H5 at position 39 cM of RGR and PR, respectively. Results from the simpler model test also indicate that lines with the B allele have on average an increased light reflectance, which translate into lighter bud-end fry color. While there is no double-simplex SNP (H4, H5) near 39 cM, adjacent SNPs map on chromosome 11 of the reference genome close to hexokinase (HXK-11, PGSC0003DMG400013187) and invertase/pectin methylesterase inhibitor family protein (InvInh-3/1, PGSC0003DMG400045232) genes.

A recent study on cold-induced sweetening in diploid potato, Braun et al. [22] reported two QTLs associated with resistance to the accumulation of reducing sugars, one on chromosome 4 and another on chromosome 6. The flanking SNP markers of these QTLs (Table 3 in Braun et al. [22]) encompass the SNPs solcap_c2_12976 and solcap_snp_c2_56145 reported here for chromosomes 4 and 6, respectively. Considering the high degree of concordance between linkage maps based on the Infinium 8303 Potato Array and the reference potato genome sequence [2, 31], it is reasonable to assume that these are homologous QTL regions. Yet, the SNP alleles underlying positive effects in Braun et al. [22] need to be determined.

Four individuals within the mapping population (i.e., A05141–010, A05141–101, A05141107, A05141–208) carried the optimal theoretical allele configuration for the four candidate SNPs associated with glucose concentration and bud and stem end fry color. These are, AAAA for solcap_snp_c2_12976 on chromosome 4, AABB for solcap_snp_c2_50301 on chromosome 5, AAAA for solcap_snp_c2_56145 on chromosome 6, and AAAB for solcap_snp_c2_48118 on chromosome 10. On average, reflectance values for these individuals were in the upper quartile of the trait scores, while glucose concentrations were in the lower quartile (Additional file 4: Table S3). Although information about the actual phase of these alleles in the progeny cannot be determined with the current methodology, the range of phenotypic values for these genotypes support our findings.

This study identified six important regions on five chromosomes (4, 5, 6, 10, 11) that are associated with tuber sugar concentration and tuber processing quality. On chromosome 5 there are two non-overlapping, maturity independent QTLs for glucose and bud end fry color, respectively (Table 3). On chromosome 4, there are overlapping QTLs for tuber glucose content, and bud and stem end fry color. They have the same peak and the same SNP explains the majority of phenotypic variation. There was a similar pattern with a cluster of overlapping QTLs on chromosome 6. Although the different QTL profile peaks (separated by 10 cM) might suggest that there is some independence between glucose and fry color. The 2-LOD support intervals completely overlapped and a single SNP explained the highest proportion of the phenotypic variance for all traits. QTLs for glucose and bud end fry color co-localized on chromosome 10. In all these cases, the direction of allele effects for glucose and fry color was coherent.

Conclusions

Linkage analysis and QTL mapping in a russet mapping population (A05141) using SNP dosage information successfully identified favorable alleles and candidate SNP markers for resistance to the accumulation of tuber reducing sugars after cold storage. These novel markers have a high potential for the improvement of tuber processing quality at the tetraploid level, especially within the russet market class represented by the parents of the mapping population. The ability to process into French fries with the light color desired by consumers is a key attribute within the russet market class. Moreover, the discovery of different genetic models for traits with overlapping QTLs at the maturity locus on chromosome 5 clearly suggests an independent genetic control.

Methods

Plant material

The mapping population used in this study (A05141) consisted of 162 F1 progeny derived from a cross between ‘Rio Grande Russet’ (female parent) and ‘Premier Russet’ (male parent). ‘Premier Russet’ is a dual purpose (processing and fresh market) potato variety with significant resistance to the accumulation of reducing sugars after long-term storage at low temperatures [24]. Rio Grande Russet is a high yielding fresh market potato variety [32]. A05141 was developed at the USDA-ARS Small Grains and Potato Germplasm Research Unit (Aberdeen, ID) and was one of the mapping populations phenotyped and genotyped by the Solanaceae Coordinated Agricultural Project (SolCAP) [1].

Phenotypic data

Field experiments were conducted at the USDA-ARS, University of Idaho Research and Extension Center (Aberdeen, ID) and North Carolina State University Tidewater Research Station (Plymouth, NC) during the summers of 2010 and 2011. The progeny and parents were evaluated in a randomized complete block design with two replications of ten-hill plots at each location. Vine maturity was visually scored on a 1–9 scale at 95 (95DAP) and 120 (120DAP) days after planting, with 1 = very early, 100% vines in a plot dead, while 9 = very late, full bloom with buds evident. Plant growth habit was recorded at the beginning of flowering and expressed as 3 = erect (> 45° with ground), 5 = semi-erect (30–45° with ground), and 7 = spreading (Additional file 10: Table S5). Verticillium wilt and early blight (EB) (Altenaria tenuis) diseases were rated on a 0–9 scale. For EB, 0 indicates no disease and 9 > 90% of the leaf area affected in a plot. For VW, the scale is a combined rating for leaf area/stem flagging due to the disease (Additional file 10: Table S5). Plots were not inoculated, but Aberdeen fields have high inoculum of Verticillium in the soil, built up over years. Early blight is observed yearly at the ID and NC locations and for this trial the disease was allowed to develop naturally, with no fungicides applied for its control.

Storage conditions

Tuber samples were treated with Isopropyl Chlorpropham (N-(3-chlorophenyl) carbamate) at a rate of 22 ppm to prevent sprouting and stored at 42 °F and 95% relative humidity for 3 months.

Glucose determinations

Glucose concentrations for the treatments in this trial were determined from a ten-tuber sample using the method of Sowokinos et al. [33] with modifications. Tubers were cut using a Keen Kut Shoe Stringer French fry cutter. Two hundred grams of tuber tissue collected from the center of the ten tubers were macerated in an Acme Juicerator (Acme Equipment, Spring Hill, FL). During processing, tuber tissue was washed with 150 mL of sodium-phosphate buffer (0.05 M, pH 7.5) for a final homogenate volume of 275 mL. Glucose concentrations were determined using a YSI model 2700 Analyzer (Yellow Springs Instrument Co., Inc., Yellow Springs, OH) and expressed on a percent fresh weight basis.

Bud-end and Stem-end fry color

A strip of tissue (3.0 cm × 0.8 cm) removed from each of the ten tubers used in the sugar extraction procedure was used for fry color determination (10 strips per replicate). Strips were fried in canola oil at 191 °C for 3.5 min. Fry color was determined within 3 min using a model 577 Photovolt Reflection Meter (model 577, Photovolt Instruments Inc., Minneapolis, MN). A green filter was used and the instrument was calibrated using a black-cavity standard as 0.0% reflectance and a white plaque (Cat. No. 26-5Cat. No. 26–570-08) as 99.9% reflectance. Measurements were taken on the ends of each strip. A relationship between USDA fry color and photovolt reflectance is described in Additional file 10: Table S5. The lower the reflectance measurement, the darker the fry color and vice versa.

Statistical analyses

Phenotypic values for individual clones within each year and location were estimated using the lsmeans (least square means) statement in SAS GLM procedure (SAS Institute, NC). Least squared means were later used for QTL mapping. For traits with high heritability (> 0.70), the lsmeans were obtained for each location over the 2 years of observation. Trait correlations were calculated and plotted using the corrplot package [34] in R Core Team [35]. Variance components of the trait scores were estimated by the restricted maximum likelihood method with years as fixed effects. Results were used to calculate broad sense heritability according to equation.
$$ {H}^2=\frac{\sigma_g^2}{\left({\sigma}_g^2+\frac{\sigma_{g^{\ast }y}^2}{m}+\frac{\sigma_e^2}{m}\right)} $$

where \( \left({\sigma}_g^2\right) \), \( \left(\frac{\sigma_{g^{\ast }y}^2}{m}\right) \), and \( \left(\frac{\sigma_e^2}{rm}\right) \) are the genetic, genotype × year interaction, and residual variance components, m is the number of years, and r is the number of replications.

SNP genotyping

DNA was extracted from freeze-dried leaves using the QIAGEN DNeasy Plant Mini Kit (QIAGEN, Germantown, MD), quantified with the Quant-iT PicoGreen assay (Invitrogen, San Diego, CA), and adjusted to a concentration of 50 ng·μL− 1. Genotyping was performed with the Infinium 8303 Potato Array as previously described [2]. The Illumina GenomeStudio software (Illumina, Inc., San Diego, CA) was used for initial sample quality assessment and for obtaining SNP theta scores (which contain information about the allele dosages for the parents and offspring).

SNP data processing

The SNP positions of the Infinium 8303 Potato Array are based on the PGSC Version 4.03 Pseudomolecules [23]. Only SNPs uniquely mapped to the reference genome were used in this study. Three filtering steps were applied to identify SNPs suitable for the analysis [4, 7]. First, SNPs were retained if the trimmed range between the 2% and 98% quantiles was equal or greater than 0.1. This step excluded monomorphic markers and removed outliers. Second, the theta values were modeled as a smooth function using locally weighted regression to remove SNPs showing a significant trend and/or differences between plates (with P < 0.0001) [7, 36]. The last step removed SNPs with missing theta values.

Genotype calling and allele dosage estimation

Five genotype configurations of bi-allelic markers (AAAA, AAAB, AABB, ABBB, BBBB) are possible in a full-sib autotetraploid mapping population. To call genotypes based on allele dosage information, we fit mixture models to the theta scores distribution of the offspring, given the parental configurations as mixture proportions [7]. The five SNP genotypes were coded as 0, 1, 2, 3, and 4, respectively. To detect double reduction products, the theta scores with potential double reduction products were regressed on the dosage class and residuals were examined to see if possible double reduction products were present [7]. SNPs with maximum residual values greater than 0.25 were examined for potential double reduction products. An additional filtering step was performed to exclude SNPs with > 10% missing genotype calls.

Linkage map construction and QTL analysis

Analyses for linkage map construction and QTL interval mapping were performed using the software TetraploidSNPMap version 1.0.5 [6]. This program is designed to handle allele dosage from SNP data. Steps for linkage mapping including testing for segregation ratios, SNP clustering, estimation of recombination fractions and LOD scores, SNP ordering, and inference of the parental phase followed the methodology developed by Hackett and collaborators [5, 7, 37]. Briefly, a SNP was removed as distorted if the significance of the χ2 goodness-of-fit statistic was less than 0.001 for simplex SNPs and 0.01 for duplex and higher dosage SNPs. To identify the alleles aligned to each of the eight homologous chromosomes, we ran the phase analysis and manually completed the phase information. Homologous chromosomes were defined as H1-H4 for parent 1 (‘Rio Grande Russet’) and H5-H8 for parent 2 (‘Premier Russet’). A thousand permutations were run to assess the statistical significance for QTLs. All possible simple models were further explored to determine the best fit to the data. The model with the lowest value for the minimum Schwarz information criterion (SIC) was selected [5]. MapChart 2.3 software [38] was used to generate linkage maps and QTL positions.

A high degree of concordance has been documented between linkage maps based on the Infinium 8303 Potato Array and the reference potato genome sequence [2, 31]. Here, we tested the concordance between the reference genome (PGSC Version 4.03 Pseudomolecules) and the linkage maps by comparing the genetic position (cM) with the physical location (Mb) of each SNP marker using the Marey map approach as implemented in the R package MareyMap version 1.3 [39]. The slope of the curve was obtained using the ‘cubic splines’ interpolation method.

Abbreviations

Chr: 

Chromosome

CI: 

Confidence interval

cM: 

Centimorgan

LOD: 

Logarithm of Odds

Mb: 

Mega base

PGSC: 

Potato Sequence Genome Consortium assembly

QTL: 

Quantitative Trait Loci

R2

Explained variance

SNP: 

Single Nucleotide Polymorphism

SolCAP: 

Solanaceae Coordinated Agricultural Project

Declarations

Acknowledgments

The authors thank Dr. Christine Hackett at the James Hutton Institute, Invergowrie, Dundee, United Kingdom, for her helpful feedback on the utilization of TetraploidSNPMap.

Mention of proprietary product or company is included for the reader’s convenience and does not imply any endorsement or preferential treatment by the USDA-ARS.

Funding

This work was supported by funds from the U.S. Department of Agriculture, National Institute of Food and Agriculture, Agriculture and Food Research Initiative Plant Breeding, Genetics, and Genome grant 2009–85606-05673.

Availability of data and materials

All of the supporting data were included as additional files.

Authors’ contributions

DSD, RGN, and GCY planned and designed the project, ANM drafted the manuscript, ANM and NMC conducted linkage mapping and QTL analysis, ANM, NMC, JC, and KH conduced phenotypic data analysis, DSD, RGN, and GCY were involved field experiments. PCB, SG, and TLB analyzed sugar content and determined fry color reflectance. All authors contributed to the editing of the manuscript and approved the final draft.

Ethics approval and consent to participate

Not applicable.

Consent for publication

Not applicable.

Competing interests

The authors declare that they have no competing interests.

Publisher’s Note

Springer Nature remains neutral with regard to jurisdictional claims in published maps and institutional affiliations.

This is a U.S. Government work and not under copyright protection in the US; foreign copyright protection may apply.

Authors’ Affiliations

(1)
USDA, Agricultural Research Service, National Peanut Research Laboratory, 1011 Forrester DR., SE, Dawson, GA 39842, USA
(2)
Plant, Soil and Microbial Sciences Department, Michigan State University, 1066 Bogue St, East Lansing, MI 48824, USA
(3)
USDA, Agricultural Research Service, 10300 Baltimore Ave, Beltsville, MD 20705-1249, USA
(4)
USDA, Agricultural Research Service, Vegetable Crops Research Unit, 1575 Linden Drive, Madison, WI 53706, USA
(5)
J.R. Simplot Company, Boise, ID 83707, USA
(6)
Department of Soil, Water, and Climate, 1991 Upper Buford Circle, Saint Paul, MN 55108, USA
(7)
Department of Horticulture, North Carolina State University, 214 Kilgore Hall, Raleigh, NC 27695, USA
(8)
USDA Agricultural Research Service, Small Grains and Potato Germplasm Research Unit, 1693 S 2700 W, Aberdeen, ID 83210, USA

References

  1. Douches D, Hirsch CN, Manrique-Carpintero NC, Massa AN, Coombs J, Hardigan M, Bisognin D, De Jong W, Buell CR. The contribution of the solanaceae coordinated agricultural project to potato breeding. Potato Res. 2014;57(3–4):215–24.View ArticleGoogle Scholar
  2. Felcher KJ, Coombs JJ, Massa AN, Hansey CN, Hamilton JP, Veilleux RE, Buell CR, Douches DS. Integration of two diploid potato linkage maps with the potato genome sequence. PLoS One. 2012;7(4):e36347.View ArticleGoogle Scholar
  3. Hamilton J, Hansey C, Whitty B, Stoffel K, Massa A, Van Deynze A, De Jong W, Douches D, Buell CR. Single nucleotide polymorphism discovery in elite north American potato germplasm. BMC Genomics. 2011;12(1):302.View ArticleGoogle Scholar
  4. Massa AN, Manrique-Carpintero NC, Coombs JJ, Zarka DG, Boone AE, Kirk WW, Hackett CA, Bryan GJ, Douches DS. Genetic linkage mapping of economically important traits in cultivated tetraploid potato (Solanum tuberosum L.). G3 Genes Genomes Genet. 2015;5(11):2357–64.Google Scholar
  5. Hackett C, Bradshaw J, Bryan G. QTL mapping in autotetraploids using SNP dosage information. Theor Appl Genet. 2014;127(9):1885–904.View ArticleGoogle Scholar
  6. Hackett CA, Boskamp B, Vogogias A, Preedy KF, Milne I. TetraploidSNPMap: software for linkage analysis and QTL mapping in autotetraploid populations using SNP dosage data. J Hered. 2017;108(4):438–42.View ArticleGoogle Scholar
  7. Hackett CA, McLean K, Bryan GJ. Linkage analysis and QTL mapping using SNP dosage data in a tetraploid potato mapping population. PLoS One. 2013;8(5):e63939.View ArticleGoogle Scholar
  8. Lindqvist-Kreuze H, Gastelo M, Perez W, Forbes GA, de Koeyer D, Bonierbale M. Phenotypic stability and genome-wide association study of late blight resistance in potato genotypes adapted to the tropical highlands. Phytopathology. 2014;104(6):624–33.View ArticleGoogle Scholar
  9. Schumann MJ, Zeng Z-B, Clough ME, Yencho GC. Linkage map construction and QTL analysis for internal heat necrosis in autotetraploid potato. Theor Appl Genet. 2017;130(10):2045–56.View ArticleGoogle Scholar
  10. da Silva W, Ingram J, Hackett CA, Coombs JJ, Douches D, Bryan G, De Jong W, Gray S. Mapping loci that control tuber and foliar symptoms caused by PVY in autotetraploid potato (Solanum tuberosum L.). G3 Genes Genomes Genet. 2017;7(11):3587–95.Google Scholar
  11. Manrique-Carpintero NC, Coombs JJ, Pham GM, Laimbeer FPE, Braz GT, Jiang JM, Veilleux RE, Buell CR, Douches DS. Genome reduction in tetraploid potato reveals genetic load, haplotype variation, and loci associated with agronomic traits. Front Plant Sci. 2018;9:944.View ArticleGoogle Scholar
  12. Li L, Paulo MJ, Strahwald J, Lübeck J, Hofferbert HR, Tacke E, Junghans H, Wunder J, Draffehn A, Van Eeuwijk F. Natural DNA variation at candidate loci is associated with potato chip color, tuber starch content, yield and starch yield. Theor Appl Genet. 2008;116(8):1167–81.View ArticleGoogle Scholar
  13. Menendez CM, Ritter E, Schäfer-Pregl R, Walkemeier B, Kalde A, Salamini F, Gebhardt C. Cold sweetening in diploid potato: mapping quantitative trait loci and candidate genes. Genetics. 2002;162(3):1423–34.PubMedPubMed CentralGoogle Scholar
  14. Geigenberger P. Regulation of starch biosynthesis in response to a fluctuating environment. Plant Physiol. 2011;155(4):1566–77.View ArticleGoogle Scholar
  15. Sowokinos JR. Biochemical and molecular control of cold-induced sweetening in potatoes. Am J Potato Res. 2001;78(3):221–36.View ArticleGoogle Scholar
  16. Bradshaw J, Hackett C, Pande B, Waugh R, Bryan G. QTL mapping of yield, agronomic and quality traits in tetraploid potato (Solanum tuberosum subsp. tuberosum). Theor Appl Genet. 2008;116(2):193–211.View ArticleGoogle Scholar
  17. Douches DS, Freyre R. Identification of genetic factors influencing chip color in diploid potato (Solanum spp.). Am Potato J. 1994;71(9):581–90.View ArticleGoogle Scholar
  18. Gebhardt C, Li L, Pajerowska-Mukhtar K, Achenbach U, Sattarzadeh A, Bormann C. Candidate gene approach to identify genes underlying quantitative traits and develop diagnostic markers in potato. Crop Sci. 2007;47(Suppl 3):106–11.Google Scholar
  19. Menéndez CM, Ritter E, Schäfer-Pregl R, Walkemeier B, Kalde A, Salamini F, Gebhardt C. Cold sweetening in diploid potato: mapping quantitative trait loci and candidate genes. Genetics. 2002;162(3):1423–34.PubMedPubMed CentralGoogle Scholar
  20. Rak K, Bethke PC, Palta JP. QTL mapping of potato chip color and tuber traits within an autotetraploid family. Mol Breed. 2017;37(2):15.View ArticleGoogle Scholar
  21. Sołtys-Kalina D, Szajko K, Sierocka I, Śliwka J, Strzelczyk-Żyta D, Wasilewicz-Flis I, Jakuczun H, Szweykowska-Kulinska Z, Marczewski W. Novel candidate genes AuxRP and Hsp90 influence the chip color of potato tubers. Mol Breed. 2015;35:224.View ArticleGoogle Scholar
  22. Braun SR, Endelman JB, Haynes KG, Jansky SH. Quantitative trait loci for resistance to common scab and cold-induced sweetening in diploid potato. Plant Genome. 2017;10(3):1–9.View ArticleGoogle Scholar
  23. Hirsch CN, Hirsch CD, Felcher K, Coombs J, Zarka D, Van Deynze A, De Jong W, Veilleux RE, Jansky S, Bethke P, Douches DS, Buell CR. Retrospective view of north American potato (Solanum tuberosum L.) breeding in the 20th and 21st centuries. G3 Genes Genomes Genet. 2013;3(6):1003–13.Google Scholar
  24. Novy RG, Whitworth JL, Stark JC, Love SL, Corsini DL, Pavek JJ, Vales MI, James SR, Hane DC, Shock CC, Charlton BA, Brown CR, Knowles NR, Pavek MJ, Brandt TL, Olsen N. Premier russet: a dual-purpose, potato cultivar with significant resistance to low temperature sweetening during long-term storage. Am J Pot Res. 2008;85(3):198–209.View ArticleGoogle Scholar
  25. Kloosterman B, Abelenda JA, Gomez MMC, Oortwijn M, de Boer JM, Kowitwanich K, Horvath BM, van Eck HJ, Smaczniak C, Prat S, Visser RGF, Bachem CWB. Naturally occurring allele diversity allows potato cultivation in northern latitudes. Nature. 2013;495(7440):246–50.View ArticleGoogle Scholar
  26. Schreiber L, Nader-Nieto AC, Schönhals EM, Walkemeier B, Gebhardt C. SNPs in genes functional in starch-sugar interconversion associate with natural variation of tuber starch and sugar content of potato (Solanum tuberosum L.). G3 Genes Genomes Genet. 2014;4(10):1797–811.Google Scholar
  27. Draffehn AM, Meller S, Li L, Gebhardt C. Natural diversity of potato (Solanum tuberosum) invertases. BMC Plant Biol. 2010;10:271.View ArticleGoogle Scholar
  28. Fridman E, Zamir D. Functional divergence of a syntenic invertase gene family in tomato, potato, and arabidopsis. Plant Physiol. 2003;131(2):603–9.View ArticleGoogle Scholar
  29. Draffehn AM, Durek P, Nunes-Nesi A, Stich B, Fernie AR, Gebhardt C. Tapping natural variation at functional level reveals allele specific molecular characteristics of potato invertase Pain-1. Plant Cell Environ. 2012;35(12):2143–54.View ArticleGoogle Scholar
  30. Li L, Tacke E, Hofferbert HR, Lubeck J, Strahwald J, Draffehn AM, Walkemeier B, Gebhardt C. Validation of candidate gene markers for marker-assisted selection of potato cultivars with improved tuber quality. Theor Appl Genet. 2013;126(4):1039–52.View ArticleGoogle Scholar
  31. Sharma SK, Bolser D, De Boer J, Sønderkær M, Amoros W, Carboni MF. Construction of reference chromosome-scale pseudomolecules for potato: integrating the potato genome with genetic and physical maps. G3 Genes Genomes Genet. 2013;3:2031–47.Google Scholar
  32. Holm D, Davidson R, Essah S. Rio Grande Russet: A new high quality fresh market russet. Colorado Spud Items. 2004;1(2):5.Google Scholar
  33. Sowokinos JR, Shock CC, Stieber TD, Eldredge EP. Compositional and enzymatic changes associated with the sugar-end defect in russet Burbank potatoes. Am J Pot Res. 2000;77(1):47–56.View ArticleGoogle Scholar
  34. Wei T, Simko V. R package “corrplot”: Visualization of a Correlation Matrix (Version 0.84). 2017. https://github.com/taiyun/corrplot.
  35. Team RC. R: a language and environment for statistical computing. Vienna: R Foundation for Statistical Computing; 2018.Google Scholar
  36. Cleveland WS. Robust locally weighted regression and smoothing scatterplots. J Am Stat Assoc. 1979;74(368):829–36.View ArticleGoogle Scholar
  37. Preedy KF, Hackett CA. A rapid marker ordering approach for high-density genetic linkage maps in experimental autotetraploid populations using multidimensional scaling. Theor Appl Genet. 2016;129(11):2117–32.View ArticleGoogle Scholar
  38. Voorrips RE. MapChart: software for the graphical presentation of linkage maps and QTLs. J Hered. 2002;93(1):77–8.View ArticleGoogle Scholar
  39. Rezvoy C, Charif D, Gueguen L, Marais GA. MareyMap: an R-based tool with graphical interface for estimating recombination rates. Bioinformatics. 2007;23(16):2188–9.View ArticleGoogle Scholar

Copyright

© ᅟ 2018

Advertisement