Search
Submit Manuscript
2.8CiteScore
Advanced Search
MAP  >  Forestry Research
2024 Volume 4
Article Contents
Abstract
INTRODUCTION
RESULTS
Genetic parameters
Characteristics of DNA pools, and variant calling
FBD-associated genomic regions, and candidate gene screening
WCI-associated genomic regions, and candidate genes screening
Homologous regions associated with FBD and WCI
Differentially expressed genes in R6:MdMYB10 apple fruit
DISCUSSION
Apple Type 1 red flesh is a complex trait influenced by interplay between several TF families
Multiple transcription factors play key roles in FBD
The nexus between WCI and FBD
MATERIALS AND METHODS
Generation of red-fleshed seedling populations
Generation of R6:MdMYB10 apple fruit
Real-Time quantitative RT-qPCR analysis
Phenotypic data analysis
Identification of extreme phenotypes, pool-sequencing, and QC of sequence data
Association test using DNA pools
Candidate gene identification
ACKNOWLEDGMENTS
Conflict of interest
Rights and permissions
References
About this article
Next Previous
ARTICLE   Open Access    

Soil parameters affecting longleaf pine (Pinus palustris) site quality in east Texas

  • 1.

    Arthur Temple College of Forestry and Agriculture, Stephen F. Austin State University, Nacogdoches, Texas 75965, USA

  • 2.

    Texas Parks and Wildlife Department, Nacogdoches, Texas 75965, USA

More Information
  • Received: 25 August 2023
    Revised: 12 December 2023
    Accepted: 19 December 2023
    Published online: 12 January 2024
    Forestry Research  4,  Article number: e002 (2024)  |  Cite this article
  • The decline since European colonization in longleaf pine (Pinus palustris Mill.) within its range in the southeastern United States, attributed to factors including both site conversion and fire exclusion has spurred interest in the re-establishment of the species. Land that originally supported longleaf pine in the southeastern United States has often been converted for agricultural use, loblolly pine (Pinus taeda Mill.) plantations, and urban development. Longleaf pine was found on a wide range of soil properties due to frequent fires which kept many competing species suppressed; fire has often been excluded due to human health, safety, and liability concerns. Longleaf pine ecosystem restoration efforts might be best focused on soils that have characteristics that naturally restrain herbaceous and hardwood competition. Properties of three soil series in east Texas that historically or are currently supporting longleaf pine ecosystems were evaluated. Analysis of Variance, Principal Component Analysis, and regression techniques were used to compare soil properties; while all three soils historically supported longleaf pine, they vary in texture, depth to argillic horizons, nutrient availability, available water capacity, and other parameters which are likely related to site quality, as measured by site index. Longleaf pine site index is influenced by depth to E and the first argillic B horizons, B horizon texture and nutrients. B horizon physical and chemical variables appear to be the most influential for longleaf pine site index on these sites, and should be considered when evaluating potential sites for longleaf pine restoration efforts.

  • INTRODUCTION

    The competition for consumer preference for fresh apples (Malus domestica) from exotic and tropical fruits is intense. Red-fleshed (RF) apple may not only provide a novel point of differentiation and enhanced visual quality, but also a source of increased concentration of potentially health-benefiting compounds within both the fresh fruit and snack/juice markets[1]. Two different types of RF apples have been characterised: Type 1 RF apple has red colouration not only in the fruit core and cortex, but also in vegetative tissues, including stems and leaves; Type 2 RF apples display red pigment only in the fruit cortex[1, 2]. To facilitate trade and lengthen the supply-window, harvested fruit are usually cold stored, which can induce a series of disorders, including physiological breakdown manifesting as a flesh browning disorder (FBD) in RF apples[3, 4]. FBD in RF apples can be caused by senescence, and there is also some evidence to suggest that a large proportion of RF apples are chilling-sensitive (Jason Johnston, Plant & Food Research Hawke's Bay, personal communication).

    Earlier studies suggested that Type 1 RF colour was determined by a promoter mutation of MdMYB10 that has a tandem replication of a myeloblastosis (MYB) binding cis-element (R6) within the promoter, resulting in autoregulation of MdMYB10[5]. However, van Nocker et al.[6] observed a large variation in the degree and pattern of red pigmentation within the cortex among the accessions carrying MdMYB10, and concluded that the presence of this gene alone was not sufficient to ensure the RF colour. A genome-wide association study (GWAS)[3], reported that, in addition to the MdMYB10 gene, other genetic factors (e.g. MdLAR1, a key enzyme in the flavonoid biosynthetic pathway) were associated with RF colour, too. Wang et al.[7] reported that many of the up-regulated genes in RF apples were associated with flavonoid biosynthesis (e.g., chalcone synthase (CHS), chalcone isomerase (CHI), dihydroflavonol 4-reductase (DFR), anthocyanin synthase (ANS), UDP-glucosyltransferase (UGT) and MYB transcription factors). Recently, MdNAC42 was shown to share similar expression patterns in RF fruit with MdMYB10 and MdTTG1, and it interacts with MdMYB10 to participate in the regulation of anthocyanin synthesis in the RF apple Redlove®[8].

    Several transcription factor genes (e.g., MYB, WRKY, bHLH, NAC, ERF, bZIP and HSF) were reported to be differentially expressed during cold-induced morphological and physiological changes in 'Golden Delicious' apples[9]. A study by Zhang et al.[10] showed that ERF1B was capable of interacting with the promotors of anthocyanin and proanthocyanidin (PA) MYB transcription factors, and suggested that ethylene regulation and anthocyanin regulation might be linked in either direction. It was reported that ethylene signal transduction pathway genes or response genes, such as ERS (ethylene response sensor), EIN3 (ethylene-insensitive3) and ERFs (ethylene response factors), may play an important role in the regulatory network of PA biosynthesis[11].

    Espley et al.[12] observed no incidence of FBD in cold-stored fruit of 'Royal Gala', but over-expression of MdMYB10 in 'Royal Gala' resulted in a high rate of FBD in RF fruit, which was hypothesised to be caused by elevated fruit ethylene concentrations before harvest and more anthocyanin, chlorogenic acid (CGA) and pro-cyanidins in RF fruit. In addition, the MYB10 transcription factor was shown to elevate the expression levels of MdACS, MdACO, and MdERF106 ethylene-regulating genes[12]. To elucidate the mechanism regulating the FBD of RF apples, Zuo et al.[13] analysed the transcriptome of tree-ripe apples at 0, 0.5 and 4 h after cutting, and reported that the differentially expressed genes at different sampling points were mainly related to plant–pathogen interactions.

    GWAS is a powerful technique for mining novel functional variants. One of the limitations of GWAS, using SNP arrays, is that they require genotyping of large numbers of individuals, which may be expensive for large populations. DNA pooling-based designs (i.e., bulk segregant analysis) test differential allele distributions in pools of individuals that exhibit extreme phenotypes (XP) in bi-parental populations, large germplasm collections or random mating populations[14−16]. In addition to reducing the number of samples to be genotyped, the use of whole genome sequencing (WGS)-based XP-GWAS has the potential to identify small-effect loci and rare alleles via extreme phenotypic selection.

    In this WGS-based XP-GWAS, we investigated the genetic basis of RF and FBD by sequencing the pools of individuals that exhibited extreme phenotypes for these two traits, and analysed the differences in allele frequencies between phenotypic classes. This method combines the simplicity of genotyping pools with superior mapping resolution. We also examined the transcriptome from transgenic apple fruit harbouring the R6:MYB10 promoter as a model for red flesh in apple. Differences in gene expression of a highly pigmented line were compared with expression in control fruit and these genes were then used for comparison with the seedling population. Understanding the genetic basis of the link between RF and FBD will help in design of strategies for selection against FBD in high-quality Type 1 RF apple cultivars.

    RESULTS

    Genetic parameters

    A snapshot of visual variation in FBD and WCI is presented in Fig. 1. The average WCI and FBD across all ~900 seedlings ranged from 0 to 7, and from 0% to 58%, respectively. Based on the MLM analysis, the estimated narrow-sense heritability (h2) of WCI and FBD was 0.57 (standard error = 0.18) and 0.09 (standard error = 0.05), respectively. The estimated genetic correlation between WCI and FBD was 0.58, and several fruit quality traits displayed unfavourable correlations with WCI (Supplemental Fig. S1). Seedlings with higher WCI scores were generally characterised by poor firmness and crispness, plus higher astringency and sourness. Estimated phenotypic and genetic correlations between all pairs of traits are listed in Supplemental Table S1.

    Figure 1.  Transverse cross-sections of apple slices showing range in (a) flesh colouration, and (b) flesh browning disorder for Type 1 red-fleshed apple. The weighted cortical intensity (WCI) scores (0−9 scale) and the proportion of the cortex area showing symptoms of flesh browning disorder are also displayed.
    DownLoad: Full-Size Img PowerPoint

    Characteristics of DNA pools, and variant calling

    A few seedlings had no red pigment in the cortex, but the average WCI score across all seedlings was 2.25. About two-thirds of the seedlings did not display any FBD symptoms, but among the remaining seedlings, FBD ranged between 1% and 58% (Fig. 2). The average WCI score for the 'low' and 'high' WCI pool was 0.45 and 5.2, respectively, while the average FBD was 0% and 20.6% for the 'low' and 'high' FBD pool, respectively (Supplemental Table S2). The average WCI score of the seedlings in the FBD pools was similar (Low: 5.5; High: 4.6), while the average FBD of the high- and low-WCI pools was 4.5% and 0.1%, respectively.

    Figure 2.  The distribution of weighted cortex intensity (WCI) scores (a) and the internal flesh browning disorder IFBD%; (b) in the population of ~900 apple seedlings. The green and red circles highlight the individuals used to form the 'low' and 'high' pools of samples.
    DownLoad: Full-Size Img PowerPoint

    FBD-associated genomic regions, and candidate gene screening

    After filtering, about 204,000 SNPs were used and the average sequencing depth of SNP loci was similar for the two pools (42 vs 44). There was a near-perfect correlation between the Z-test statistics and G-statistics, so only the latter are discussed hereafter. A plot of the G' values, smoothed over 2 Mb windows, is shown for all 17 chromosomes (Chrs) in Supplemental Fig. S2. XP-GWAS identified genomic regions significantly associated with FBD on 12 out of the 17 Chrs (Fig. 3), and putative candidate genes within ±1.0 Mb distanceof the significant G' peaks were identified (Table 1). Additional genomic regions, which did not meet the significance threshold but displayed distinguished G' peaks, were also identified across all chromosomes (Supplemental Table S3).

    Figure 3.  G'-statistics across the linkage groups (LG) showing significant association with the flesh browning disorder (FBD) in apple The horizontal red lines indicate the significance threshold. The putative candidate genes (refer to Table 1) underpinning various G' peaks are also shown.
    DownLoad: Full-Size Img PowerPoint
    Table 1.  A list of the genomic regions associated with internal flesh browning disorder (FBD) in apples. Putative candidate genes residing within these regions are also listed using GDDH13v1.1 reference genome assembly.
    ChrGenomic region (Mb)Putative genes functions
    221.9–23.5Ethylene-responsive element binding factor 13 (MdERF13: MD02G1213600);
    33.1–4.9cinnamate 4-hydroxylase (C4H) enzyme (MD03G1051100, MD03G1050900 and MD03G1051000); MdWRKY2: MD03G1044400; MdWRKY33 (MD03G1057400)
    38.2–9.6ascorbate peroxidase 1 (MdAPX1: MD03G1108200, MD03G1108300)
    314.7–16.6senescence-related MdNAC90 (MD03G1148500)
    336.5–37.5Ethylene response sensor 1 (MdERS1: MD03G1292200); flavonoid biosynthesis protein MdMYB12 (MD03G1297100); Heat shock protein DnaJ (MD03G1296600, MD03G1297000); pectin methylesterase (MdPME) inhibitor protein (MD03G1290800, MD03G1290900, MD03G1291000).
    411.0–13.0phenylalanine and lignin biosynthesis protein MdMYB85 (MD04G1080600)
    422.1–24.4MYB domain protein 1 (MD04G1142200); HSP20-like protein (MD04G1140600); UDP-glucosyltransferase (UGT) proteins UGT85A7 (MD04G1140700, MD04G1140900); UGT85A3 (MD04G1140800); UGT (MD04G1141000, MD04G1141300); UGT85A2 (MD04G1141400); UGT85A4 (MD04G1141500); DNAJ heat shock protein (MD04G1153800, MD04G1153900, MD04G1154100)
    427.6–29.6HCT/HQT regulatory genes MD04G1188000 and MD04G1188400
    629.8–31.7Volz et al. (2013) QTL for IFBD; anthocyanin regulatory proteins MdMYB86 (MD06G1167200); triterpene biosynthesis transcription factor MdMYB66 (MD06G1174200); Cytochrome P450 (MD06G1162600; MD06G1162700, MD06G1162800, MD06G1163100, MD06G1163300, MD06G1163400, MD06G1163500, MD06G1163600, MD06G1163800, MD06G1164000, MD06G1164100, MD06G1164300, MD06G1164400, MD06G1164500, MD06G1164700)
    713.9–15.9heat shock protein 70B (MD07G1116300)
    717.4–19.0Drought-stress WRKY DNA-binding proteins (MdWRKY56: MD07G1131000, MD07G1131400)
    723.6–25.2MdPAL2 (MD07G1172700); drought-stress gene NGA1 (MD07G1162400); DNAJ heat shock family protein (MD07G1162300, MD07G1162200), stress-response protein (MdNAC69: MD07G1163700, MD07G1164000)
    97.9–11.8MD09G1110500 involved in ascorbate oxidase (AO); MdUGT proteins (MD09G1141200, MD09G1141300, MD09G1141500, MD09G1141600, MD09G1141700, MD09G1141800, MD09G1142000, MD09G1142500, MD09G1142600, MD09G1142800, MD09G1142900, MD09G1143000, MD09G1143200, MD09G1143400) involved in flavonoids biosynthesis; heat shock proteins 89.1 (MD09G1122200) and HSP70 (MD09G1137300); Ethylene-forming enzyme MD09G1114800; Anthocyanin regulatory protein MdNAC42 (MD09G1147500, MD09G1147600)
    914.9–16.5triterpene biosynthesis transcription factor protein MdMYB66 (MD09G1183800);
    111.7–3.0ethylene response factor proteins (MdEIN-like 3: MD11G1022400)
    1138.6–40.6Senescence-related gene 1 (MD11G1271400, MD11G1272300, MD11G1272000, MD11G1272100, MD11G1272300, MD11G1272400, and MD11G1272500); Chalcone-flavanone isomerase (CHI) protein (MD11G1273600) and MdbHLH3 (MD11G1286900, MDP0000225680); cytochrome P450 enzyme (MD11G1274000, MD11G1274100, MD11G1274200, MD11G1274300, MD11G1274500, and MD11G1274600); heat shock transcription factor A6B (MdHSFA6B: MD11G1278900) – involved in ABA-mediated heat response and flavonoid biosynthesis.
    122.2–3.5heat shock protein 70-1 ((MD12G1025600, MD12G1025700 and MD12G1026300) and heat shock protein 70 (MD12G1025800 and MD12G1025900 and MD12G1026000); ethylene (MD12G1032000) and auxin-responsive (MD12G1027600) proteins.
    127.2–8.4DNAJ heat shock domain-containing protein (MD12G1065200 and MD12G1067400)
    1327.5–30.5MD13G1257800 involved in polyphenol 4-coumarate:CoA ligase (4CL)
    1337.5–39.5pectin methyl esterase inhibitor superfamily protein MdPMEI (MD13G1278600)
    1425.4–27.5Drought-stress response gene MdWRKY45 (MD14G1154500); chalcone synthase (CHS) family proteins (MD14G1160800 and MD14G1160900); triterpene biosynthesis transcription factor MdMYB66 (MD14G1180700, MD14G1181000, MD14G1180900); anthocyanin biosynthesis protein (MdMYB86: MD14G1172900); cytochrome P450 proteins (MD14G1169000, MD14G1169200, MD14G1169600, MD14G1169700)
    150–1.5Ethylene synthesis proteins (MD15G1020100, MD15G1020300 and MD15G1020500); dihydroflavonol reductase (DFR) gene (MD15G1024100)
    154.6–6.8MdMYB73 (MD15G1076600, MD15G1088000) modulates malate transportation/accumulation via interaction with MdMYB1/10; MdNAC52: MD15G1079400) regulates anthocyanin/PA; heat shock transcription factor B4 (MD15G1080700); stress-response protein MdWRKY7 (MD15G1078200)
    1553.5–54.9MdC3H (MD15G1436500) involved in chlorogenic acid biosynthesis; MdEIN3 (MD15G1441000) involved in regulating ethylene synthesis and anthocyanin accumulation
    168.9–10.9SAUR-like auxin-responsive protein (MD16G1124300) and MdNAC83: (MD16G1125800) associated with fruit ripening; MD16G1140800 regulates proanthocyanidin; MdPAE genes (MD16G1132100, MD16G1140500) regulates ethylene production.
     | Show Table
    DownLoad: CSV

    The ethylene-responsive factor 13 (MdERF13: MD02G1213600) resided within the significant region (21.9–23.5 Mb) on Chr2, while the ascorbate peroxidase 3 (MdAPX3: MD02G1127800, MD02G1132200) resided in the prominent region between 9.45 and 11.28 Mb) (Fig. 3, Table 1). There were several genomic regions showing association with FBD on Chr3. The first region (8.2–9.5 Mb) flanked MdAPX1 (MD03G1108200, MD03G1108300), while MdERF3 (MD03G1194300) and heat shock protein 70 (HSP70: MD03G1201800, D03G1201700) resided within the prominent peak region (25.5–27.5 Mb). Another significant region (36.5–37.5 Mb) at the bottom on Chr3 harboured ethylene response sensor 1 (MdERS1: MD03G1292200), a flavonoid-biosynthesis related protein MdMYB12 (MD03G1297100), HSP DnaJ and pectin methyl esterase (PME) inhibitor proteins (Fig. 3, Table 1).

    The significant G' region (27.6–29.6 Mb) on Chr4 flanked the genes MD04G1188000 and MD04G1188400 involved in the biosynthesis of hydroxycinnamoyl CoA shikimate/quinate hydroxycinnamoyl transferase (HCT/HQT), while the adjacent region (22.1–24.3 Mb) harboured a cluster of UDP-glucosyltransferase (UGT) proteins and HSP (Table 1). Another region between 11.0 and 13.0 Mb encompassed phenylalanine and lignin biosynthesis gene MYB85 (MD04G1080600)[17]. The only significant region associated with FBD on Chr6 spanned between 29.8 and 31.7 Mb, which included a SNP earlier reported associated with FBD[4]. This region also harboured several MYB proteins (MdMYB86: MD06G1167200; MdMYB98: MD06G1172900; MdMYB66: MD06G1174200) and a large cluster of cytochrome P450 proteins (Table 1).

    A significant FBD-associated region (13.9–15.9 Mb) on Chr7 encompassed the HSP 70B (MD07G1116300), while another significant region between 23.6 Mb and 25.2 Mb flanked the gene coding for phenylalanine ammonia lyase 2 (MdPAL2: MD07G1172700), a drought stress gene NGA1 (MD07G1162400), DNAJ HSP, and stress-response protein (MdNAC69: MD07G1163700, MD07G1164000) (Fig. 3, Table 1). A large significant region spanning between 7.9 and 11.8 Mb on Chr9 encompassed the gene MD09G1110500 putatively involved in ascorbate oxidase (AO), HSP (HSP70: MD09G1137300; HSP89.1: MD09G1122200), an ethylene-forming enzyme (MD09G1114800), and MdNAC42 (MD09G1147500, MD09G1147600). Another significant genomic region on Chr9 was between 14.9 and 16.5 Mb, which harbours the MYB domain protein MdMYB66 (MD09G1183800) (Table 1).

    A sharp G' peak region (1.7–3.0 Mb) on Chr11 associated with FBD encompassed ethylene insensitive 3 (MdEIN3: MD11G1022400) along with a cluster of UGT proteins, WD-40, and bHLHL proteins (Table 1). A significant region between 38.6 and 40.6 Mb at the bottom of Chr11 was dominated by clusters of senescence-related genes and cytochrome P450 enzymes. This genomic region also flanked a chalcone-flavanone isomerase (CHI) family protein (MD11G1273600) and a bHLH protein (MdbHLH3: MD11G1286900, MDP0000225680), along with the heat shock transcription factor A6B (HSFA6B: MD11G1278900) (Table 1, Supplemental Table S3).

    The region (2.2–3.3 Mb) associated with FBD on Chr12 flanked genes for the HSP 70 and 70-1, NAC proteins, ethylene and auxin-responsive proteins (Table 1). An adjacent significant region (7.2–8.4 Mb) flanked DNAJ HSP, along with bZIP, bHLH and WD-40 repeat-like proteins (Table 1; Supplemental Table S3). There was a large genomic region on Chr13 showing a significant association with FBD. In this region, the first G' peak (27.5–30.5 Mb) harboured MD13G1257800, which regulates polyphenol 4-coumarate: CoA ligase (4CL) synthesis. The second G' peak region (32.6–34.7 Mb) corresponded to an earlier mapped QTL for flesh browning in white-fleshed apples[18].

    The significant region (25.4–27.5 Mb) on Chr14 harboured various genes for proteins with different putative functions, such as AP2 proteins, WD-40 repeat family proteins, bHLH proteins, MdNAC83 (MD14G1150900), MdWRKY45 (MD14G1154500), chalcone synthase (CHS) family proteins (MD14G1160800 and MD14G1160900), cytochrome P450 proteins, and several MYB domain proteins (MdMYB86: MD14G1172900; MdMYB98: MD14G1179000; MdMYB66: MD14G1180700, MD14G1181000, MD14G1180900) (Fig. 3, Table 1, Supplemental Table S3). A sharp G' peak (15.2–17.2 Mb) on Chr14 did not meet the significance threshold corresponding to the FBD QTL in white-fleshed apples[18].

    The upper 1.5 Mb region on Chr15 associated with FBD encompassed several transcription factor families, including WD-40 repeats, bHLH, bZIP, ethylene synthesis proteins, and dihydroflavonol reductase (DFR) protein (MD15G1024100) (Fig. 3, Table 1; Supplemental Table S3). Another significant region on Chr15 (4.6−6.8 Mb) harboured HSF B4, MdMYB73 (which interacts with MdMYB1/10 to modulate malate transportation) and MdNAC52 (MD15G1079400), which regulates anthocyanin and PA synthesis by directly regulating MdLAR[19]. A significantly associated region at the bottom of Chr15 (53.5–54.9 Mb) harboured MdNAC35 (MD15G1444700), MdC3H (MD15G1436500) involved in the production of p-coumarate 3-hydroxylase (C3H) enzyme, which plays a role in chlorogenic acid biosynthesis, and MdEIN3 (MD15G1441000).

    The significant region between 8.9 and 10.9 Mb on Chr16 harboured a gene for SAUR-like auxin-responsive protein (MD16G1124300) and MdNAC83 (MD16G1125800), both of which have been reportedly associated with apple fruit ripening[20]. This region also encompassed MD16G1140800, which regulates PA[11], and a gene for the pectin acetyl esterase protein MdPAE10 (MD16G1132100) involved in ethylene production and shelf-life[21]. A sharp G' peak region (20.4–22.4 Mb) in the middle of Chr16 flanked MdERF1B (MD16G1216900) and MdMYB15 (MD16G1218000 and MD16G1218900), involved in altering anthocyanin and PA concentrations[10] (Fig. 3, Table 1, Supplemental Table S3).

    WCI-associated genomic regions, and candidate genes screening

    The average sequencing depth of the SNP loci (~160,000) retained for marker-trait association was similar for the two WCI pools (41 vs 44). The significant regions were located on Chrs 2, 4, 6, 7, 10, 15 and 16 (Fig. 4). The genomic intervals within ±1.0 Mb of the significant G' peaks, and the putative candidate genes within those intervals, are listed in Table 2. Additional genomic regions, which did not meet the significance threshold but displayed distinct G' peaks, were also identified across most chromosomes ( Supplemental Fig. S2, Supplemental Table S3). A significant region on Chr4 encompassed chalcone synthase (CHS) genes, along with an ERF (MD04G1009000) involved in regulating PA biosynthesis[11].

    Figure 4.  G’-statistics across the linkage groups (LG) showing significant association with the weighted cortex intensity (WCI) in apple. The horizontal red lines indicate the significance threshold. The putative candidate genes (refer to Table 2) underpinning various G' peaks are also shown.
    DownLoad: Full-Size Img PowerPoint
    Table 2.  A list of the genomic regions significantly associated with the weighted cortex intensity (WCI) in apples. Putative candidate genes residing within these regions are also listed using GDDH13v1.1 reference genome assembly.
    ChrGenomic region (Mb)Putative genes functions
    211.2–13.2UDP-glucosyltransferase (UGT) proteins (MD02G1153000, MD02G1153100, MD02G1153200, MD02G1153300; MD02G1153400; MD02G1153500; MD02G1153700; MD02G1153800; MD02G1153900);
    40–1.2Chalcone synthase (CHS) genes (MD04G1003000; MD04G1003300 and MD04G1003400); DNAJ heat shock protein (MD04G1003500); MdERF (MD04G1009000) involved in regulating PA biosynthesis.
    69.5–11.0Ubiquitin protein (MD06G1061100); stress-response WRKY protein MdWRKY21 (MD06G1062800);
    612.1–13.6pectin methylesterase inhibitor superfamily protein (MdPME: MD06G1064700); phenylalanine and lignin biosynthesis gene (MdMYB85)
    616.0–17.6MD06G1071600 (MDP0000360447) involved in leucoanthocyanidin dioxygenase (LDOX) synthesis
    624.0–26.0UDP-glycosyltransferase proteins (MD06G1103300, MD06G1103400, MD06G1103500 and MD06G1103600); auxin response factor 9 (MdARF9: MD06G1111100) and heat shock protein 70 (Hsp 70; MD06G1113000).
    74.1–5.6Ethylene insensitive 3 protein (MdEIN3: MD07G1053500 and MD07G1053800) involved in proanthocyanidins (PA) biosynthesis; ubiquitin-specific protease (MD07G1051000, MD07G1051100, MD07G1051200, MD07G1051500, MD07G1051300, MD07G1051700 and MD07G1051800).
    1015.9–18.3bHLH proteins (MD10G1098900, MD10G1104300, and MD10G1104600); UGT protein (MD10G1101200); UGT 74D1 (MD10G1110800), UGT 74F1 (MD10G1111100), UGT 74F2 (MD10G1111000).
    1035.6–38.6Stress-response WRKY proteins (MdWRKY28: MD10G1266400; MdWRKY65: MD10G1275800). ethylene responsive factors (MdERF2: MD10G1286300; MdERF4: MD10G1290400, and MdERF12: MD10G1290900) and a NAC domain protein (MdNAC73: MD10G1288300); polyphenol oxidase (PPO) genes ((MD10G1298200; MD10G1298300; MD10G1298400; MD10G1298500; MD10G1298700; MD10G1299100; MD10G1299300; MD10G1299400).
    1524.8–26.8heat shock factor 4 (MD15G1283700), drought-stress response WRKY protein 7 (MdWRKY7: MD15G1287300), MdMYB73 (MD15G1288600) involved in ubiquitination and malate synthesis
    1531.7–34.2MYB domain protein 93 (MdMYB93: MD15G1323500) regulates flavonoids and suberin accumulation (Legay et al. 2016); ubiquitin -specific protease 3 (MD15G1318500),
    161.5–3.4MYB domain proteins (MD16G1029400) regulates anthocyanin; senescence-associated gene 12 (MD16G1031600); ethylene response factor proteins (MdERF118: MD16G1043500, and MD16G1047700 (MdRAV1); MdMYB62 ( MD16G1040800) flavonol regulation; malate transporter MdMa2 (MD16G1045000: MDP0000244249), Ubiquitin-like superfamily protein (MD16G1036000),
    165.4–7.4MdMYB88 (MD16G1076100) regulates phenylpropanoid synthesis and ABA-mediated anthocyanin biosynthesis; MdMYB66 (MD16G1093200) regulates triterpene biosynthesis, MdWRKY72 (MD16G1077700) mediates ultraviolet B-induced anthocyanin synthesis.
     | Show Table
    DownLoad: CSV

    There were distinct G' peaks within a large genomic region (spanning between 9.0 and 18.0 Mb) significantly associated with WCI on Chr6 (Fig. 4). The G' peak region 12.1–13.6 Mb encompassed MdMYB85 (MD06G1064300), and a gene for a pectin methyl esterase

    (PME) inhibitor protein (MD06G1064700); while the region between 16.0 and 17.6 Mb flanked the genes involved in leucoanthocyanidin dioxygenase (LDOX) synthesis (Table 2). The auxin response factor 9 (MdARF9: MD06G1111100) and the HSP70 (MD06G1113000) resided in the significant region between 24.0 and 26.0 Mb on Chr6 (Table 2).

    The genomic region (4.1–5.6 Mb) with significant association with WCI on Chr7 flanked MdEIN3 (MD07G1053500, MD07G1053800), which plays an important role in the regulatory network of PA biosynthesis[11]. A significant region (35.6–38.7 Mb) on Chr10 encompassed gene clusters for bHLH and WRKY proteins along with an ethylene repressor factor (MdERF2: MD10G1286300) (Fig. 4, Table 2; Supplemental Table S3). This region also harboured MdERF4 (MD10G1290400), MdERF12 (MD10G1290900) and NAC domain genes (MdNAC73: MD10G1288300), which have been reported to be associated with fruit ripening[20]. In addition, there was a cluster of polyphenol oxidase (PPO) genes residing in this region (Table 2).

    On Chr15, the significant WCI-associated region spanning between 31.7 Mb and 34.2 Mb encompassed several TF families, including MdMYB93 (MD15G1323500) (Fig. 4, Table 2, Supplemental Table S3). A distinguished G' peak region (24.9–26.9 Mb) on Chr15 harboured genes for WD-40 repeat-like proteins, bHLH proteins, redox responsive transcription factor 1 (MD15G1283200), cytochrome P450 proteins, HSF4 (MD15G1283700), and MdWRKY7 (MD15G1287300). The MdMYB73 (MD15G1288600) residing in this region has been shown to interact with MdMYB1/10 and regulates several functions, including cold-stress response, ubiquitination and malate synthesis[22, 23].

    The significant 1.5–3.4 Mb region on Chr16 flanked an anthocyanin repressor MYB protein (MD16G1029400), senescence-associated gene 12 (MD16G1031600), malate transporter MdMa2 (MD16G1045000: MDP0000244249), and two ethylene response factor genes (MdERF118: MD16G1043500; MdRAV1: MD16G1047700), which interact in retaining flesh firmness[24]. The MdMYB62 (MD16G1040800) gene residing in this region is phylogenetically linked to MdMYB8 (MD06G1217200), which plays a major role in flavonoid biosynthesis[25]. Another significant region (5.4–7.4 Mb) on Chr16 harboured several bHLH genes, including MdMYB88 (MD16G1076100) involved in phenylpropanoid synthesis resulting in drought resistance[26] and ABA-mediated anthocyanin production[27]. The MdMYB66 (MD16G1093200) gene, which is involved in triterpene biosynthesis, and the MdWRKY72 (MD16G1077700) gene involved in anthocyanin synthesis, were also present in this region (Table 2, Fig. 4).

    Homologous regions associated with FBD and WCI

    The genomic regions tagged by XP-GWAS were enriched for regulatory functions, and several of these have a paralog (Table 3). For example, paralogues for a trio of genes (MdMYB86, MdMYB98 and Cytochrome 450) resided in the FBD-associated genomic regions on Chr6 and Chr14. Several pairs of genes resided together in the FBD-associated paralog regions; for example, MdNAC90 and germin-like protein 10 on Chr3 and Chr11; WRKY55 and WRKY70 on Chr1 and Chr7; and ERF1 and ERF5 on Chr4 and Chr6. The ethylene response factor MdERF1B had a paralog in the FBD-associated region on Chr13 and Chr16. Interestingly, three copies of MdNAC83 (Chrs 14, 16 17) and four copies of MdMYB66 (Chrs 4, 6, 9 and 14) resided in the FBD-associated regions. A pair of genes (RAV1 and MYB62) resided together in the WCI-associated paralog regions on Chrs13 and 16 (Table 3).

    Table 3.  A list of homologues genes/genomic regions associated with the internal flesh browning disorder (IFBD) and red flesh (WCI) in apples. Putative candidate genes residing within these regions are also listed using GDDH13v1.1 reference genome assembly.
    TraitChr (genomic region: Mb)Gene namePredicted gene IDPutative function
    IFBDChr6 (29.8–31.7 Mb)MdMYB66MD06G1174200Suberin/triterpene deposition
    IFBDChr14 (25.4–27.5 Mb)MdMYB66MD14G1180700, MD14G1181000, MD14G1180900Suberin/triterpene deposition
    IFBDChr4 (7.0–8.1 Mb)MdMYB66MD04G1060200Suberin/triterpene deposition
    IFBDChr9 (14.9–16.5 Mb)MdMYB66MD09G1183800Suberin/triterpene deposition
    IFBDChr6 (29.8–31.7 Mb)MdMYB86MD06G1167200Anthocyanin regulation
    IFBDChr14 (25.4–27.5 Mb)MdMYB86MD14G1172900Anthocyanin regulation
    IFBDChr6 (29.8–31.7 Mb)MdMYB98MD06G1172900Drought stress response
    IFBDChr14 (25.4–27.5 Mb)MdMYB98MD14G1179000Drought stress response
    IFBDChr6 (29.8–31.7 Mb)Cytochrome P450MD06G1162600, MD06G1162700, MD06G1162800, MD06G1163100, MD06G1163300, MD06G1163400, MD06G1163500, MD06G1163600, MD06G1163800, MD06G1164000, MD06G1164100, MD06G1164300, MD06G1164400, MD06G1164500, MD06G1164700Flavonoid and triterpenic metabolism
    IFBDChr14 (25.4–27.5 Mb)Cytochrome P450MD14G1169000, MD14G1169200, MD14G1169600, MD14G1169700Flavonoid and triterpenic metabolism
    IFBDChr3 (14.7–16.6 Mb)MdNAC90MD03G1148500Senescence-related
    IFBDChr11 (16.5–18.5 Mb)MdNAC90MD11G11679000Senescence-related
    IFBDChr14 (25.4–27.5 Mb)MdNAC83MD14G1150900Senescence/ripening-related
    IFBDChr16 (8.9–10.9 Mb)MdNAC83MD16G1125800Senescence/ripening-related
    IFBDChr17 (0–0. 7 Mb)MdNAC83MD17G1010300Senescence/ripening-related
    IFBDChr3 (14.7–16.6 Mb)Germin-like protein 10MD03G1148000Polyphenol oxidase
    IFBDChr11 (16.5–18.5 Mb)Germin-like protein 10MD11G1167000, MD11G1167100, MD11G1167400, MD11G1169200Polyphenol oxidase
    IFBDChr1 (26.9–28.5 Mb)MdWRKY55MD01G1168500Drought stress response
    IFBDChr7 (30.3–31.9 Mb)MdWRKY55MD07G1234600Drought stress response
    IFBDChr1 (26.9–28.5 Mb)MdWRKY70MD01G1168600Drought stress response
    IFBDChr7 (30.3–31.9 Mb)MdWRKY70MD07G1234700Drought stress response
    IFBDChr4 (7.0–8.1 Mb)MdERF1MD04G1058000Ethylene responsive factor
    IFBDChr6 (5.5–7.3 Mb)MdERF1MD06G1051800Ethylene responsive factor
    IFBDChr4 (7.0–8.1 Mb)MdERF5MD04G1058200Ethylene responsive factor
    IFBDChr6 (5.5–7.3 Mb)MdERF5MD06G1051900Ethylene responsive factor
    IFBDChr13 (18.0–20.0 Mb)MdERF1BMD13G1213100Ethylene response factor 1
    IFBDChr16 (20.4–22.4 Mb)MdERF1BMD16G1216900Ethylene response factor 1
    WCIChr13 (2.6–4.4 Mb)MdRAV1MD13G1046100Ethylene responsive factor
    WCIChr16 (1.5–3.4 Mb)MdRAV1MD16G1047700Ethylene responsive factor
    WCIChr13 (2.6–4.4 Mb)MdMYB62MD13G1039900Flavonol biosynthesis
    WCIChr16 (1.5–3.4 Mb)MdMYB62MD16G1040800Flavonol biosynthesis
    IFBDChr9 (7.9–11.8 Mb)MdNAC42MD09G1147500, MD09G1147600Anthocyanin accumulation
    WCIChr17 (11.4–12.4 Mb)MdNAC42MD17G1134400Anthocyanin accumulation
    IFBDChr9 (7.9–11.8 Mb)HSP 70MD09G1137300Heat stress response
    WCIChr17 (11.4–12.4 Mb)HSP 70MD17G1127600Heat stress response
    IFBDChr4 (11.0–13.0 Mb)MdMYB85MD04G1080600Phenylalanine and lignin biosynthesis
    WCIChr6 (12.1–13.6 Mb)MdMYB85MD06G1064300Phenylalanine and lignin biosynthesis
    IFBDChr15 (13.2–14.2 Mb)MdEBF1MD15G1171800Ethylene inhibition
    WCIChr8 (15.2–17.2 Mb)MdEBF1MD08G1150200Ethylene inhibition
    IFBDChr15 (53.5–54.9 Mb)MdEIN3MD15G1441000Ethylene insensitive 3 protein
    WCIChr8 (30.3–31.6 Mb)MdEIN3MD08G1245800Ethylene insensitive 3 protein
    IFBDChr11 (1.7–3.0 Mb)MdEIN3MD11G1022400Ethylene insensitive 3 protein
    WCIChr7 (4.1–5.6 Mb)MdEIN3MD07G1053500, MD07G1053800Ethylene insensitive 3 protein
    IFBDChr15 (4.6–6.8 Mb)MdMYB73MD15G1076600, MD15G1088000Cold-stress response & malate accumulation
    WCIChr15 (24.8–26.8 Mb)MdMYB73MD15G1288600Cold-stress response & malate accumulation
    IFBDChr15 (4.6–6.8M b)MdWRKY7MD15G1078200Anthocyanin accumulation
    WCIChr15 (24.8–26.8 Mb)MdWRKY7MD15G1287300Anthocyanin accumulation
    IFBDChr15 (43.5–45.5 Mb)MdMYB93MD15G1369700Flavonoid & suberin accumulation
    WCIChr15 (31.7–34.2 Mb)MdMYB93MD15G1323500Flavonoid & suberin accumulation
     | Show Table
    DownLoad: CSV

    Paralogs of several regulatory functions were also found in the regions associated with either FBD or WCI; for example, a significant region (7.9–11.8 Mb) harbouring a gene trio (MdNAC42, HSP70 and HSP89.1) on Chr9 was associated with FBD, but the paralogs of this trio also resided in a distinct G' region associated with WCI on Chr17 (Table 3, Supplemental Table S3). Some other examples included MdMYB85 (Chr4 for FBD, and Chr6 for WCI), MdEBF1 (Chr8 for WCI, and Chr15 for FBD), MdEIN3 (Chr8 for WCI, and Chr15 for FBD), and EIN3 (Chr7 for WCI, and Chr11 for FBD). A pair of genes (MdMYB73 and MdWRKY7) resided together in the separate regions associated with FBD (4.6–6.8 Mb) and WCI (24.8–26.8 Mb) within Chr15 (Table 3).

    Differentially expressed genes in R6:MdMYB10 apple fruit

    The red pigmentation in fruit flesh differed amongst transgenic lines, with fruit from line A10 presenting the most deeply pigmented tissues (Supplemental Fig. S3), while those from lines A2 and A4 were similar in having a lower intensity of pigmentation. No pigmentation was observed in the flesh of control 'Royal Gala' (RG) fruit. RNAseq analysis of a representative transgenic line (A2) compared with RG revealed that a total of 1,379 genes were differentially expressed (log 2-fold), with 658 genes upregulated and 721 genes downregulated. This list was then assessed for commonality with the genomic regions and candidate genes from the XP-GWAS.

    Genes that contained mis-sense SNPs (which would result in a change in predicted protein) and that were also differentially expressed between A2 red-fleshed transgenic line and white-fleshed 'Royal Gala' (control) apples included anthocyanin-related flavanone 3-hydroxylase, chalcone synthase and dihydroflavonol 4-reductase (Table 2 & 4). Many more upstream DNA variants (e.g. in potential promoter-controlling elements) were seen in this group of differentially expressed genes (DEGs) that were also in regions underlying WCI or FBD. Genes encoding enzymes that may be involved in FBD, such as a Rho GTPase activating protein, peroxidase, lipoxygenase 1, and ethylene-forming enzyme (ACO4), were DEGs and showed mis-sense SNPs, including a potential stop codon in the Rho GTPase-activating protein (Table 4). This intersection between DNA change and differential expression warrants further research to evaluate the functions of these genes.

    Table 4.  List of candidate genes associated with WCI or FBD in the XP-GWA and R6:MdMYB10 apple datasets.
    Mutations in GWAS apple populationPredicted gene functionExpression in R6 and ‘Royal Gala’ apples
    GeneMissense
    SNPs    		SNP
    Stop    		Upstream
    variants    		DatasetLocusAnnotation and TAIR IDAverage
    RPKM
    R6 flesh    		Average RPKM
    WT flesh    		log2Fold
    Change
    MD02G11322001FBD Supplemental Table S3Chr02:9450615-11289014flavanone 3-hydroxylase
    (F3H, TT6, F3'H) AT3G51240    		387752.47
    MD02G11336003FBD Supplemental Table S3Chr02:9450615-11289014fatty acid desaturase 5
    (FAD5) AT3G15850    		18109.98
    MD02G11537001WCI Table 2Chr02:11278254-13234556UDP-Glycosyltransferase, lignin related AT2G1856010385341.03
    MD03G10592002FBD Supplemental Table S3Chr03:3177994-4910866Peroxidase AT5G053402514.21
    MD03G1143300422FBD Table 1Chr03:14797661-16611554bZIP transcription factor (DPBF2, AtbZIP67) AT3G44460114852-2.80
    MD03G1147700121FBD Table 1Chr03:14797661-16611554Rho GTPase activating protein AT5G22400205891.29
    MD04G1003400331WCI Table 2Chr04:1-1284894Chalcone synthase (CHS, TT4) AT5G1393014433122.34
    MD04G120410063FBD Table 1Chr04:27629536-29612282lipoxygenase 1 (LOX1, ATLOX1) AT1G550205462681.09
    MD06G1160700129FBD Table 1Chr06:29862341-31737341peptide met sulfoxide reductase AT4G251302121254562.04
    MD06G116140014FBD Table 1Chr06:29862341-31737341Pectin lyase-like protein AT5G63180614129530-2.22
    MD07G1240700126FBD Supplemental Table S3Chr07:30357332-31979025Fe superoxide dismutase 2 AT5G511001051826-3.99
    MD07G13069006FBD Supplemental Table S3Chr07:34607570-36531467UDP-glucosyl transferase 78D2 AT5G170506771112.78
    MD08G12491002WCI Supplemental Table S3Chr08:30486569-31607516HSP20-like chaperone (ATHSP22.0) AT4G102503186456.12
    MD09G111400032FBD Table 1Chr09:7999212-11883202fatty acid desaturase 5 (FAD5) AT3G158507608.68
    MD09G11468001FBD Table 1Chr09:7999212-11883202PHYTOENE SYNTHASE (PSY) AT5G1723038456160901.32
    MD10G132810045FBD Supplemental Table S3Chr10:40235258-41736791ethylene-forming enzyme (ACO4) AT1G050108769033910331.21
    MD15G10236001FBD Table 1Chr15:1-1487288jasmonic acid carboxyl methyltransferase (JMT) AT1G19640472111552.07
    MD15G10241008FBD Table 1Chr15:1-1487288dihydroflavonol 4-reductase (DFR, TT3, M318) AT5G428008723091.63
    MD17G11334004WCI Supplemental Table S3Chr17:11422073-12433463PHYTOENE SYNTHASE (PSY) AT5G172306001721.87
    MD17G12606002FBD Supplemental Table S3Chr17:31098955-32776972dehydroascorbate reductase 1 (DHAR3) AT5G167101867−1.78
     | Show Table
    DownLoad: CSV

    DISCUSSION

    Apple Type 1 red flesh is a complex trait influenced by interplay between several TF families

    Several gene families (e.g. chalcone synthase, UGT, anthocyanin synthase, HSP, PAL, ERF, WRKY proteins, ABA, and bZIP transcription factors) residing in WCI-associated genomic regions have been reported to be associated with RF in apple[7, 12, 28]. Several of these genes are implicated in stress responses, suggesting that flavonoids and anthocyanin biosynthesis could also be associated with stress (e.g., drought, water loss) tolerance of RF apples[7].

    MdMYB73, which may play a role in the cold-stress response[23, 29], resided in the WCI-associated regions. Ethylene is among the various modulators of environmental stresses induced by factors such as drought and cold temperatures[9, 28]. Several ERF proteins (e.g. MdERF1B, MdERF3) interact with the promoters of MYB domain proteins to regulate anthocyanin and proanthocyanidin (PA) accumulation in apple[10, 11, 30]. The WCI-associated genomic regions in our study flanked several ERFs previously reported to have roles in anthocyanin and PA regulation. MdERF4 (MD10G1290400), one of the three ERFs residing in the WCI-associated region on Chr10, has high phylogenetic similarity with MdERF38, which interacts with MdMYB1 to regulate drought-related anthocyanin biosynthesis[31]. Some ERFs (e.g. MdERF2: MD10G1286300; MdERF4: MD10G1290400, and MdERF12: MD10G1290900) and MdNAC73 residing in the WCI-associated regions have been reported to be associated with fruit maturation[20] – suggesting an interaction between ethylene production and RF colour[12].

    The WCI-associated region on Chr15 included a gene MdMYB73 (MD15G1288600) involved in malate acid synthesis. Previous studies have reported that MdMYB1 regulates both anthocyanin and malate accumulation, and perhaps MdBT2 regulates MdMYB73-mediated anthocyanin accumulation[22, 23]. The malate transporter gene MdMa2 and the MdMYB7 gene, involved in the regulation of anthocyanin and flavanols in RF apple[32], resided together in the WCI-associated upper region of LG16. Taken together, these results support the hypothesis that there is interplay between anthocyanin and malate accumulation in the RF apple. There was a cluster of ubiquitin-specific proteases underpinning the WCI-associated regions on Chrs 7 and 15, which is supported by earlier reports suggesting that ubiquitin-specific proteases respond to auxin and might suppress anthocyanin biosynthesis proteins[33]. The auxin response factor 9 (MdARF9: MD06G1111100), which has also been shown to suppress anthocyanin biosynthesis in RF callus samples[33], also resided in the WCI-associated genomic region on LG7.

    Multiple transcription factors play key roles in FBD

    Seedlings in both the low- and high-FBD pools carried the MdMYB10 gene, which suggests that MdMYB10 itself is not the causal factor of FBD in RF apples. Long-term cold storage generally results in senescence-related flesh breakdown, and several transcription factor genes (e.g. MYB, WRKY, NAC, ERF, cytochrome P450, and HSP) have been shown to express differentially during long-term cold storage[9]. The FBD-associated genomic regions in our study harboured several ERFs, suggesting that ethylene synthesis proteins may be contributing to cell wall disassembly, allowing PPO enzymes to come into contact with phenolic compounds and potentially leading to FBD symptoms[18].

    Pectin methyl esterase (PME) genes resided in the FBD-associated significant regions on Chrs 3, 13 and 17. Volatile generation and senescence degradation have been suggested to be bio-markers of FBD, and the expression levels of methyl esters were found to be associated with FBD and senescence in 'Fuji' apples after cold storage[34]. The MdPME2 gene has also been reported to be associated with apple flesh firmness and mealiness[21, 35]. The co-location of genes encoding cell wall-degrading enzymes (MdPME) and QTLs for FBD has been reported in apple[18, 36].

    The clusters of cytochrome P450 enzymes and senescence-related genes resided in some of the FBD-associated regions in this study. There are no earlier reports of the involvement of P450 enzymes in apple FBD expression, but some genes related to cytochrome P450 were found to be upregulated during litchi fruit senescence[37]. Pericarp browning in litchi is mainly attributable to the degradation of anthocyanin, and the ABA-initiated oxidation of phenolic compounds by PPO[37]. Flavonoids are among the major polyphenols in RF apples[5], and cytochrome P450 is part of the regulatory mechanism for flavonoid metabolism[38]. Association of FBD with the genomic region harbouring P450 would suggest its role in enhanced polyphenol synthesis causing FBD. The co-occurrence of a senescence-related gene MdNAC90 (MD03G1148500; MD11G11679000) and the PPO regulator germin-like proteins (in the FBD-associated paralog regions on Chrs 3 and 11) lends support to an interplay between senescence and oxidation of phenolics and anthocyanins.

    A significant region on Chr9 encompassed a cluster of UGT proteins, which play a role in the regulation of flavonoids and phenolic compounds, as well as converting phloretin to phloridzin[39]. Cytochrome P450, which resided within several FBD-associated regions, is reportedly involved in flavonoid metabolism, such as chlorogenic acid (CGA) acid and phloridzin, which have also been positively associated with suberin production and cell wall disassembly[40, 41].

    Reactive oxygen species (ROS) play an important role in regulating physiological processes in plants, such as senescence[37]. Legay et al.[42] suggested that MdMYB93 (MD15G1369700), found residing in the FBD-associated regions in this study, plays a critical role in remobilisation of flavonoid/phenolic compounds, which can be utilised for detoxification of ROS in the case of oxidative stress. However, flavonoid biosynthesis has also been linked with suberin production causing cuticle cracks in apples[42, 43], and the development of cuticular cracks could accelerate flesh browning as a result of an enhanced oxidative process[44].

    Several HSP (e.g. HSP70, HSP70-1, HSP60, HSP89.1, and HSP DNAJ) and HSF (e.g. HSF4, HSFB4, HSFA6B) resided in the FBD-associated genomic regions. Ferguson et al.[45] showed that, during summer, apple flesh temperature could reach as high as 43 °C, and that an increase in the expression of HSP in apples was associated with high daily flesh temperatures, suggesting a role of HSP to counter heat stress. HSPs have been reported to interact with AP2/ERFs and to play a role in flavonoid biosynthesis and drought tolerance in apple[46]. Heat stress affects lignin accumulation and its substrate, O-phenols, and has been reported to play role in enzymatic browning[47]. Additionally, HSF that regulate HSP expression have also been reported to be regulated by cold stress to generate heat-induced cold tolerance in banana[48]. HSF1 was shown to transcriptionally regulate the promoters of HSP to enhance chilling tolerance in loquat fruit[49]. Activity of the enzymes (PAL, C4H, 4CL) of the phenylpropanoid pathway was positively correlated with loquat fruit lignification, whilst suppression of their expression by heat shock treatment and low-temperature conditioning significantly reduced fruit lignification[50].

    Wang et al.[7] showed that ascorbate peroxidase (APX) was among the genes that were upregulated in RF apple compared with in white-fleshed apples. Several genes (e.g. MdAPX1, MdAPX3, MdAPX4, and MdDHAR1) involved in ascorbate synthesis resided in the genomic regions associated with FBD. Co-localisation of MdDHAR and ascorbic acid (AsA) synthesis genes in the FBD-linked genomic regions have been reported[51], suggesting that the low AsA content increases fruit susceptibility to FBD[52]. It has been shown that the expression of MdMYB1 and MdDHAR genes was strongly correlated in RF apples, and that AO and APX were upregulated by anthocyanin regulatory genes[31].

    There were several genes associated with CGA biosynthesis residing in the FBD-associated regions on various linkage groups, including the genes MYB19 (MD07G1268000) and MdC3H (MD15G1436500). Higher concentrations of CGA were reported in transgenic apple lines carrying MdMYB10[12], suggesting a role of CGA metabolism in the expression of FBD[53]. Interestingly, some of the FBD-associated regions (e.g. Chrs 9, 11, 13, 14 and 17) reported here in RF apples coincide with those reported earlier for FBD in white-fleshed apples[18, 36, 53], suggesting some common underlying genetic mechanisms.

    The nexus between WCI and FBD

    As discussed above, ERFs have been reported to be involved in the accumulation of anthocyanin and PA biosynthesis, while ethylene synthesis proteins also contribute to cell membrane breakdown, allowing the PPO enzyme to come into contact with phenolic compounds, potentially leading to FBD symptoms[12, 18]. We observed clusters of anthocyanin biosynthesis proteins (bHLH), ERFs (MdERF2, MdERF4, MdERF12) and PPO genes together in the genomic region associated with WCI on Chr10. The co-occurrence of these gene families perhaps facilitates potential interactions that contribute to the genetic correlation between WCI and FBD. We noted that genes involved in flavonoid regulation and ethylene synthesis occurred together in the FBD-associated regions on several chromosomes (e.g. MD16G1140800 and MdPAE10: MD16G1132100; MdERF1B: MD16G1216900 and MdMYB15: MD16G1218000, MD16G1218900) – suggesting these genes could be in linkage disequilibrium and this would contribute to the expression of WCI and FBD.

    MYB7 (MD16G1029400) resided in the WCI-associated upper region of Chr16, and the expression level of MYB7 was shown to be correlated with that of LAR1 in peach fruit[54]. The MdLAR1 protein (MDP0000376284), which is located about 1.2 Mb upstream of MYB7 (MD16G1029400), was reported to be associated with WCI and FBD in apple[3]. Mellidou et al.[36] reported that 4CL (MD13G1257800 in the FBD-linked region on Chr13), which catalyses the last step of the phenylpropanoid pathway, leading either to lignin or to flavonoids, was upregulated in browning-affected flesh tissues. The gene MdMYB85, involved in the regulation of flavonoid and lignin biosynthesis, resided in the WCI-associated region on Chr6 and FBD-associated paralogous region on Chr4. Metabolic interactions between anthocyanin and lignin biosynthesis have been reported for apple[55] and strawberry[56], while flesh lignification and internal browning during low-temperature storage in a red-fleshed loquat cultivar was shown to be modulated by the interplay between ERF39 and MYB8[57].

    The co-localisation of MdMYB66 and cytochrome P450 proteins, along with the anthocyanin regulatory protein MdMYB86, in paralogous FBD-associated genomic regions on Chr6 and Chr14 suggests that these genes interact as a 'hub' contributing to the WCI-FBD genetic link. The paralogs of some other genes were found to be residing in the regions associated with either WCI or FBD. For example, MdNAC42 and HSP70 co-localised in the FBD-associated region on Chr9, but this same pair of genes also resided in the most prominent region associated with WCI on Chr17. Similarly, paralogs of MdEIN3 resided in the WCI-associated region on Chr7 (MD07G1053500, MD07G1053800) and FBD-associated region on Chr11 (MD11G1022400). Interestingly, the paralogs of MdMYB73 and MdWRKY7 co-localised in the WCI- (24.6–26.8 Mb) and FBD-associated (4.6–6.8 Mb) regions on Chr15. The WCI-associated region at the bottom of Chr11 hosted a cluster of senescence-related genes, along with the anthocyanin biosynthesis gene MdbHLH3 (MD11G1286900), suggesting they might interact in the genetic nexus between FBD and WCI.

    FBD in RF apples can be caused by senescence, injury via extreme temperature exposure (chilling or heat), or enzymatic (cut fruit) reaction. Genes reported to be connected to all three factors were located in various FBD-associated genomic regions in this study. Postharvest strategies that both delay senescence and limit exposure to low temperatures may be needed to manage FBD. We also hypothesise that high ascorbic acid content could help to minimise expression of FBD in Type-1 RF cultivars. The adverse genetic correlation between WCI and FBD appears to arise from dual and/or interactive roles of several transcription factors, which would pose challenges for designing a conventional marker-assisted selection strategy. The use of bivariate genomic BLUP to estimate breeding values to simultaneously improve adversely correlated polygenic traits (e.g. WCI and FBD), could be an alternative approach[58, 59].

    MATERIALS AND METHODS

    Generation of red-fleshed seedling populations

    A population of 900 apple seedlings composed of 24 full-sib families was generated in 2011 by selected crossings between six red-leaved pollen parents and six white-fleshed female parents. All six pollen parents inherited their red-leaf phenotype from the same great-grandparent 'Redfield'[1]. Each pollen parent was involved in four crosses, and the female parents were involved in three to six crosses each. Foliage colour of young seedlings is a phenotypic marker for Type 1 RF apple. The main purpose of this trial was to understand the flesh colour variation and FBD in the Type 1 RF seedlings, so only the seedlings with red foliage (i.e. carrying MdMYB10) were kept for this trial. The number of seedlings per family varied from 10 to 95. The seedlings were grafted onto 'M9' rootstock and were planted in duplicate at the Plant & Food Research orchard in Hawke's Bay, New Zealand (39°39′ S, 176°53′ E) in 2015.

    Phenotyping for RF and FBD was conducted over two consecutive fruiting seasons (2017 and 2018). Fruit were harvested once, when judged mature, based on a change in skin background colour from green to yellow, and when the starch pattern index (SPI) was between 1.0 and 2.0 (on a scale of 0 to 7). In each season, six fruit were harvested from each plant and stored for 70 d at 0.5 °C, followed by 7 d at 20 °C before fruit evaluation. Fruit were cut in half across the equator and the proportion of the cortex area (PRA) that was red in colour, and the intensity of the red colour (RI) (= 1 (low) to 9 (high)) was scored. A weighted cortical intensity (WCI) was then calculated (PRA × RI) as an estimation of the amount of red pigment in the fruit. The proportion of the cortex area showing symptoms of FBD was also recorded. WCI and FBD were averaged over all fruit for a particular seedling. Fruit were also assessed for the following eating quality traits on a 1 (lowest) to 9 (highest) scale: firmness, crispness, juiciness, sweetness, sourness and astringency, to understand the genetic correlations of eating-quality traits with WCI and FBD.

    Generation of R6:MdMYB10 apple fruit

    The binary vector pSA277-R6:MYB10 was transferred into Agrobacterium tumefaciens strain LAB4404 by electroporation. Transgenic 'Royal Gala' plants were generated by Agrobacterium-mediated transformation of leaf pieces, using a method previously reported[5]. Wild-type 'Royal Gala' and three independent transgenic lines (A2, A4 and A10) of R6:MdMYB10 were grown under glasshouse conditions in full potting mix with natural light. The resulting fruit were assessed for flesh colour phenotypes at harvest (around 135 d after full bloom). Fruit peel and cortex from three biological replicates were collected and frozen in liquid nitrogen, with each replicate compiled from five pooled mature fruit for each transgenic line or wild-type control.

    Real-Time quantitative RT-qPCR analysis

    Total RNA of 36 samples (3 R6:MYB10 lines and 1 wild type control, 3 time points, 3 biological replicates) was extracted, using Spectrum Plant Total RNA Kit (SIGMA). Removal of genomic DNA contamination and first-strand cDNA synthesis were carried out using the mixture of oligo (dT) and random primers according to the manufacturer's instructions (QuantiTect Reverse Transcription Kit, Qiagen). Real-time qPCR DNA amplification and analysis was carried out using the LightCycler 480 Real-Time PCR System (Roche), with LightCycler 480 software version 1.5. The LightCycler 480 SYBR Green I Master Mix (Roche) was used following the manufacturer's method. The qPCR conditions were 5 min at 95 °C, followed by 45 cycles of 5 s at 95 °C, 5 s at 60°C, and 10 s at 72 °C, followed by 65 °C to 95 °C melting curve detection. The qPCR efficiency of each gene was obtained by analyzing the standard curve of a cDNA serial dilution of that gene. The expression was normalized to Malus × domestica elongation factor 1-alpha MdEF1α (XM_008367439) due to its consistent transcript levels throughout samples, with crossing threshold values changing by less than 2.

    Phenotypic data analysis

    Individual fruit measurements were first averaged for each seedling. As the phenotyping was repeated over two years, we used a mixed linear model (MLM) accounting for this 'permanent environmental effect', as previously described[58]. Pedigree-based additive genetic relationships among seedlings were taken into account for estimation of genetic parameters using ASReml software[60]. Product-moment correlations between best linear unbiased predictions (BLUP) of breeding values of all seedlings for different traits were used as estimates of genetic correlation among traits.

    Identification of extreme phenotypes, pool-sequencing, and QC of sequence data

    A selective DNA pooling procedure was adopted to construct DNA pools. A high-pool and a low-pool were constructed separately for the two traits (WCI and FBD). Genomic DNA was extracted from the leaves of selected seedlings, and quantified by fluorimetry using the picogreen reagent (Cat#P11496, Thermo). The low and high pools consisted of 35 seedlings each, and normalised amounts (~300 ng) of DNA from individuals were pooled. The pools were dried down with DNA Stable reagent (Cat#93021001, Biomatrica) in a centrifugal evaporator and shipped for sequencing. Each DNA pool was sequenced using paired-end 125 bp reads on the Illumina HiSeq 2500 platform. The quality of raw sequence reads was checked with FastQC/0.11.2 and MultiQC/1.2. Based on the quality control reports, the reads were aligned to the published apple reference genome GDDH13 v1.1[61] using the program bowtie2/2.3.4.3[62] with trimming from both ends before alignments and aligning in full read length ("-5 6 -3 5 –end-to-end"). The mapping results were marked for duplicate alignments, sorted, compressed and indexed with samtools/1.12[63]. Based on the alignment of binary alignment map (BAM) files of the high and low pools, single nucleotide polymorphism (SNP) identification was performed using samtools/1.12 ('samtools mpileup') and bcftools/1.12 ('bcftools call –mv')[63, 64]. Variant sites with missing genotype in any of the pools, or having the same genotype between the pools, were discarded. To minimise the influence of sequencing quality on association analysis, the identified SNPs were further filtered according to the following criteria: 1) a Phred-scaled quality score > 20; and 2) the read depth in each pool was neither < 35, nor > 500.

    Association test using DNA pools

    The allele frequencies between each pair of bulk DNAs (low versus high WCI; low versus high FBD) were compared at each SNP locus. Differences in the allele frequencies between the low and high pool were expected to be negligible for unlinked SNP markers, but allele frequency differences would be larger for SNPs closely linked to the underlying quantitative trait loci (QTLs) contributing to the extreme phenotypes. A nonparametric test (G-statistic = 2 × Σni ln(ni/nexp), where ni (i = 1 to 4) represented counts of reference and alternate alleles at a particular SNP generated from sequencing of the low and high pool, and nexp was the expected allele count assuming no allele frequency divergence between the two DNA pools[65].

    We then calculated a modified statistic (G'), which took into account read count variation caused by sampling of segregants as well as variability inherent in short-read sequencing of pooled samples[65]. Using R package QTLseqr[66], firstly a G-statistic was calculated for each SNP marker, and then a weighted average using Nadaraya-Watson kernel was obtained to yield a G' statistic for a sliding genomic window of 2 Mb size. The Nadaraya-Watson method weights neighbouring markers' G-values by their distance from the focal SNP so that closer SNPs receive higher weights. The 95th percentile value of G' was used as a threshold to identify significant hotspots and to identify the putative candidate genes residing within the ±1 Mb region around the G' peak. For comparison purposes, a standard two-sided Z-test[14] was also performed to determine the significance of allele frequency differences at SNP loci between the pools for each trait.

    Candidate gene identification

    The GDDH gene models intersecting with the XP-GWA hotspots were pulled out with bedtools/2.30.0 ("bedtools intersect -wo -nonamecheck"). The selected genes were further blasted to TAIR10 ("-evalue 1e-5") and the annotated functions from Arabidopsis genes with the best blast score, the highest % identity, and the longest aligned length, were used. Then the expressions of genes located in the GWA hotspots were extracted from the RNAseq analysis of the R6:MdMYB10 representative transgenic line, and the log 2-fold change between the R6:MdMYB10 and 'Royal Gala' apples were calculated.

    ACKNOWLEDGMENTS

    This research was funded in 2017/18 by the Strategic Science Investment Fund of the New Zealand Ministry of Business, Innovation and Employment (MBIE) and from 2019 by the Plant & Food Research Technology Development – Pipfruit programme. We thank our colleague Jason Johnston for providing some pictures of the flesh browning disorder in red-fleshed apples. Richard Volz and Jason Johnston provided constructive comments and suggestions on the manuscript.

  • The authors declare that they have no conflict of interest.

  • [1]

    Rapport DJ, Costanza R, McMichael AJ. 1998. Assessing ecosystem health. Trends in Ecology & Evolution 13:397−402

    doi: 10.1016/S0169-5347(98)01449-9

    CrossRef   Google Scholar

    [2]

    Chazdon RL. 2008. Beyond deforestation: restoring forest and ecosystem services on degraded lands. Science 320:1458−60

    doi: 10.1126/science.1155365

    CrossRef   Google Scholar

    [3]

    Henderson JP, Grissino-Mayer HD. 2009. Climate-tree growth relationships of longleaf pine (Pinus Palustris Mill.) in the Southeastern Coastal Plain, USA. Dendrochronologia 27:31−43

    doi: 10.1016/j.dendro.2008.08.001

    CrossRef   Google Scholar

    [4]

    Pessin LJ. 1938. Effects of soil moisture on the rate of growth of longleaf and slash pine seedlings. Plant Physiology 13:179−89

    doi: 10.1104/pp.13.1.179

    CrossRef   Google Scholar

    [5]

    Peet RK, Allard DJ. 1993. Longleaf pine vegetation of the southern Atlantic and eastern Gulf coast regions: a preliminary classification. Proceedings of the Tall Timbers Fire Ecology Conference, Tallahassee, Florida 1993. 18: 45−81. http://labs.bio.unc.edu/Peet/pubs/TTFEC-1993.pdf
    [6]

    Patterson TW, Knapp PA. 2016. Observations on a rare old-growth montane longleaf pine forest in central North Carolina, USA. Natural Areas Journal 36:153−61

    doi: 10.3375/043.036.0206

    CrossRef   Google Scholar

    [7]

    Chapman HH. 2013. Is the longleaf type a climax? Fire Ecology 9:2−7

    doi: 10.1007/BF03400625

    CrossRef   Google Scholar

    [8]

    Harcombe PA, Glitzenstein JS, Knox RG, Orzell SL, Bridges EL. 1993. Vegetation of longleaf pine region of the west gulf coastal plain. Proceedings of the Tall Timbers Fire Ecology Conference. Tallahassee, Florida 1993, 18: 83−104. https://talltimbers.org/wp-content/uploads/2018/09/83-Harcombeetal1993_op.pdf
    [9]

    Landers JL, Van Lear DH, Boyer WD. 1995. The longleaf pine forests of the southeast: requiem or renaissance. Journal of Forestry 93:39−44

    doi: 10.1093/jof/93.11.38

    CrossRef   Google Scholar

    [10]

    Drewa PB, Platt WJ, Moser EB. 2002. Community structure along elevation gradients in headwater regions of longleaf pine savannas. Plant Ecology 160:61−78

    doi: 10.1023/A:1015875828742

    CrossRef   Google Scholar

    [11]

    Oswalt SN, Smith WB, Miles PD, Pugh SA. 2012. Forest Resources of the United States, 2012: a technical document supporting the Forest Service 2010 update of the RPA Assessment. USDA Forest Service General Technical Report WO-92. https://www.fs.usda.gov/research/treesearch/47322
    [12]

    Christensen NL. 1993. The effects of fire on nutrient cycles in longleaf pine ecosystems. Proceedings of the Tall Timbers Fire Ecology Conference. Tallahassee, Florida 1993, 18: 205−14. https://talltimbers.org/wp-content/uploads/2018/09/205-Christensen1993_op.pdf
    [13]

    Outcalt KW. 1993. Southern pines performance on sandhill sites in Georgia and South Carolina. Southern Journal of Applied Forestry 17:100−02

    doi: 10.1093/sjaf/17.2.100

    CrossRef   Google Scholar

    [14]

    Stambaugh MC, Sparks JC, Abadir ER. 2014. Historical pyrogeography of Texas, USA. Fire Ecology 10:72−89

    doi: 10.4996/fireecology.1003072

    CrossRef   Google Scholar

    [15]

    Brockway DG, Lewis CE. 1997. Long-term effects of dormant-season prescribed fire on plant community diversity, structure and productivity in a longleaf pine wiregrass ecosystem. Forest Ecology and Management 96:167−83

    doi: 10.1016/S0378-1127(96)03939-4

    CrossRef   Google Scholar

    [16]

    Rosenberg N, Trajtenberg MA. 2001. A general-purpose technology at work: the Corliss steam engine in the late19th century US. The Journal of Economic History 64:61−99

    Google Scholar

    [17]

    Frost CC. 1993. Four centuries of changing landscape patterns in the longleaf pine ecosystem. Proceedings of the Tall Timbers Fire Ecology Conference. Tallahassee, Florida. 18: 17−23. https://talltimbers.org/wp-content/uploads/2018/09/17-Frost1993_op.pdf
    [18]

    Eviner VT, Hawkes CV. 2008. Embracing variability in the application of plant–soil interactions to the restoration of communities and ecosystems. Restoration Ecology 16:713−29

    doi: 10.1111/j.1526-100X.2008.00482.x

    CrossRef   Google Scholar

    [19]

    Matson GC. 1916. The Catahoula Sandstone. Report. Professional Paper 98, USGS published series, Mississippi. pp 209−26. https://doi.org/10.3133/pp98M
    [20]

    United States Climate Data. August 24, 2017. https://www.usclimatedata.com/climate/united-states/us
    [21]

    USDA Forest Service. 1929. Volume, yield, and stand tables for second-growth southern pines. Handbook 50. U.S. Department of Agriculture, Washington, DC. 202 pp. https://doi.org/10.5962/bhl.title.65546
    [22]

    Klute A. 1986. Methods of soil analysis: part 1physical and mineralogical methods. Madison, Wisconsin: American Society of Agronomy, Inc. Soil Science Society of America, Inc. xxviii, 1188 pp. http://doi.org/10.2136/sssabookser5.1.2ed
    [23]

    Soil Science Division Staff. 1993. Soil survey manual. United States Department of Agriculture Handbook No. 18. 639 pp. https://www.nrcs.usda.gov/sites/default/files/2022-09/SSM-intro.pdf
    [24]

    United States Department of Agriculture-Natural Resource Conservation Services (USDA- NRCS). 1992. Soil-woodland correlation field data sheet (ECS-005) for longleaf pine on selected Letney, Stringtown, and Tehran soils. Retrieved from https://esi.sc.egov.usda.gov/ESIForestland/.
    [25]

    Saxton KE, Rawls WJ. 2006. Soil water characte1istic estimates by texture and organic matter for hydrologic solutions. Soil Science Society America Journal 70:1569−78

    doi: 10.2136/sssaj2005.0117

    CrossRef   Google Scholar

    [26]

    Bolioluzzi EC, Tessier D, Rheinheimer DS, Julien JL. 2006. The cation exchange capacity of a sandy soil in southern Brazil: an estimation of permanent and pH-dependent charges. European Journal of Soil Science 57:356−64

    doi: 10.1111/j.1365-2389.2005.00746.x

    CrossRef   Google Scholar

    [27]

    McWhorter IV. 2005. Effects of fire exclusion on the longleaf pine ecosystem of upland island wilderness. Thesis. Stephen F. Austin State University, Nacogdoches, Texas.
    [28]

    Farrish KW, Doolittle IA, Gambel EE. 1990. Loamy substrata and forest productivity of sandy glacial drift soils in Michigan. Canadian Journal of Soil Science 70:181−87

    doi: 10.4141/cjss90-020

    CrossRef   Google Scholar

    [29]

    Hunter LR, Gibson AR. 1984. Predicting Pinus radiata site index from environmental variables. New Zeeland Journal Forestry Science 14:53−64

    Google Scholar

    [30]

    Corona P, Scotti R, Tarchiani N. 1998. Relationship between environmental factors and site index in Douglas-fir plantations in central Italy. Forest Ecology and Management 110:195−207

    doi: 10.1016/S0378-1127(98)00281-3

    CrossRef   Google Scholar

    [31]

    Coile TS. 1935. Relation of site index for shortleaf pine to certain physical prope1iies of the soil. Journal of Forestry 33:726−30

    doi: 10.1093/jof/33.8.726

    CrossRef   Google Scholar

    [32]

    Dexter AR. 2004. Soil physical quality: Part I. Theory, effects of soil texture, density, and organic matter, and effects on root growth. Geoderma201−14

    doi: 10.1016/j.geoderma.2003.09.004

    CrossRef   Google Scholar

    [33]

    Fox TR, Lee Allen H, Albaugh TJ, Rubilar R, Carlson CA. 2007. Tree nutrition and forest fertilization of pine plantations in the southern United States. Southern Journal of Applied Forestry 31:5−11

    doi: 10.1093/sjaf/31.1.5

    CrossRef   Google Scholar

    [34]

    Mead DJ, Gadgil RL. 1978. Fertilizer use in established radiate pine stands in New Zealand. New Zeeland Journal Forestry Science 8:105−34

    Google Scholar

  • Cite this article

    Oswald BP, Svehla R, Farrish KW. 2024. Soil parameters affecting longleaf pine (Pinus palustris) site quality in east Texas. Forestry Research 4: e002 doi: 10.48130/forres-0023-0031
    Oswald BP, Svehla R, Farrish KW. 2024. Soil parameters affecting longleaf pine (Pinus palustris) site quality in east Texas. Forestry Research 4: e002 doi: 10.48130/forres-0023-0031
    shu

Tables(7)

Download PDF

Article Metrics

Article views(3420) PDF downloads(468)

Other Articles By Authors

ARTICLE   Open Access    

Soil parameters affecting longleaf pine (Pinus palustris) site quality in east Texas

  • 1.

    Arthur Temple College of Forestry and Agriculture, Stephen F. Austin State University, Nacogdoches, Texas 75965, USA

  • 2.

    Texas Parks and Wildlife Department, Nacogdoches, Texas 75965, USA

  • Received: 25 August 2023
  • Revised: 12 December 2023
  • Accepted: 19 December 2023
  • Published online: 12 January 2024
Forestry Research  4,  Article number: e002  (2024)  |  Cite this article

Abstract: The decline since European colonization in longleaf pine (Pinus palustris Mill.) within its range in the southeastern United States, attributed to factors including both site conversion and fire exclusion has spurred interest in the re-establishment of the species. Land that originally supported longleaf pine in the southeastern United States has often been converted for agricultural use, loblolly pine (Pinus taeda Mill.) plantations, and urban development. Longleaf pine was found on a wide range of soil properties due to frequent fires which kept many competing species suppressed; fire has often been excluded due to human health, safety, and liability concerns. Longleaf pine ecosystem restoration efforts might be best focused on soils that have characteristics that naturally restrain herbaceous and hardwood competition. Properties of three soil series in east Texas that historically or are currently supporting longleaf pine ecosystems were evaluated. Analysis of Variance, Principal Component Analysis, and regression techniques were used to compare soil properties; while all three soils historically supported longleaf pine, they vary in texture, depth to argillic horizons, nutrient availability, available water capacity, and other parameters which are likely related to site quality, as measured by site index. Longleaf pine site index is influenced by depth to E and the first argillic B horizons, B horizon texture and nutrients. B horizon physical and chemical variables appear to be the most influential for longleaf pine site index on these sites, and should be considered when evaluating potential sites for longleaf pine restoration efforts.

    HTML

    Introduction

    • Many ecosystems have been degraded through exploitation of their natural resources, or land-use conversion to agricultural and urban use[1], and restoration is often challenging due to modifications of soils, introduction of exotic invasive species, and lack of adequate resources to adequately conduct the restoration. Site selection is an important step in ecosystem restoration because the original ecosystems may have been greatly altered due to human activities[2]. The longleaf pine (Pinus palustris Mill.) ecosystems of the southeastern United States are no exception to this degradation. Prior to European settlement, longleaf pine ecosystems occupied vast areas of the southern Atlantic and Gulf Coastal Plain regions of the United States, with approximately 30 million hectares extending between east Texas to Virginia, and stretching as far south as Florida, covering several climatic, physiographic, and many soil types[3−6].

      Longleaf pine was found in a wide range of ecosystems and sites from excessively drained sandhills to poorly drained flatwoods[7−11]. Longleaf pine was most competitive on the sandier, well-drained sites across the region; however, a relatively frequent low intensity fire return interval, every two to eight years, set by native peoples or from lightning, is regarded as a key factor in historically reducing hardwood and shrub encroachment on most sites where it was found[12−15]. During the logging and naval stores industry boom of the 1920s, old growth longleaf pine was quickly reduced[16]; by the mid-1930s, only 10% of the old growth longleaf pine forest remained in east Texas and west Louisiana, but was mostly secondary growth[17]. This dramatic decline led conservation groups and government agencies to begin conserving the remaining longleaf stands, and also to initiate longleaf pine ecosystem restoration on sites where the ecosystem once existed. However, many challenges exist that hinder this process.

      One cause of longleaf pine ecosystem restoration failure is the inadequate consideration of soil suitability for longleaf pine. Soil type can affect the vegetation present, while vegetation can affect the condition of the soils[18]. Due to these challenges, restoration efforts hypothetically should focus in areas that fit site specific soil/site parameters that support longleaf pine ecosystem restoration with the least management inputs. The objective of this study was to evaluate select soil properties on three soil mapping units (series) currently supporting longleaf pine stands in east Texas and relate these properties to longleaf pine site index.

    Materials and methods

      Study area

    • This study was conducted in the Western Gulf Region of the native longleaf pine range in eastern Texas within portions of the Angelina National Forest (31°2'52.3" N, 94°21'48.96" W) and Sabine National Forest (31°10'56.21" N, 93°43'34.68" W), United States Forest Service Forests and Grasslands of Texas (USA). Sites contained longleaf pine ecosystems before and after the logging boom, and are considered optimum reference sites for possible longleaf pine restoration. All are located on the Catahoula geologic formation, stretching from east Texas to the Mississippi River, that consisted of sandstone, ranging from a few meters to approximately 18 m thick[19]. As with most of these national forests, recurring prescribed fires on a 3−5 year interval have been used to maintain the site conditions and reduce fuel loading. All sites contained an overstory of longleaf pine, with minimal mid-story or longleaf pine advanced regeneration, and with a variety of herbaceous species and woody plants dominating the understory. The climate for east Texas is described as humid and subtropical, with mild winters with mean low temperatures in January between 2.8 °C to 3.9 °C, with summer temperatures reaching 33.3 °C to 34.4 °C for mean highs in August. Annual rainfall ranges from 1,240 to 1,510 mm with a relatively long growing season[20].

      Sampling locations were located on three different soil series mapping units, exhibiting different soil characteristics on well-drained or excessively-drained soils and ranged in depth and texture to the argillic B horizon: Letney Series (loamy, siliceous, semi-active, thennic Arenic Paleudults), Tehran Series (loamy, siliceous, semiactive, thermic Grossarenic Paleudults), and the Stringtown Series (fine-loamy, siliceous, semiactive, thennic Typic Hapludults).

      Plot establishment

    • Ten, 50 m radius plots were established within each soil series across the two national forests, for a total of 30 plots. Prior to selection, each potential plot was randomly located on relatively pure soil map units determined by soil profile assessments at five points, one point in the center and four in each cardinal direction, 50 m from the center point; verification and identification of the soil series was accomplished using a bucket auger. Any of the points that failed to be consistent with the range of characteristics for the given soil map unit for the site were rejected.

      Field sampling

    • A 10 Basal Area Factor prism was used to determine basal area at each plot center. Site index trees were chosen from the trees recorded with the prism by selecting the six closest trees to plot center that were either dominant or codominant and free of wounds. If six were not recorded by the prism, the nearest suitable trees still within the plot were measured. Annual growth rings from the six trees were quantified from a tree core extracted at DBH (Diameter at Breast Height) to determine age. A laser range finder was used to estimate total height to the nearest 3.05 cm. Total age and height were used in the site index curves developed for longleaf pine[21]. Soil samples were taken using a bucket auger at plot center to correlate with the collected longleaf pine data. Soil samples were taken from the first three horizons (A, E, and the first argillic B) while individual horizon depths were measured to a depth of 150 cm.

      Laboratory methods

    • Soil textural (sand, silt, and clay) analyses were conducted using the Bouyoucos method[22] from the A, E, and B horizon samples. For coarse textured soils, 100 g of oven-dried soil was used, while 50 g was used for medium and fine textured soils; each sample was mixed with 100 ml of sodium hexametaphosphate, left for 12 h in deionized water, and then agitated for 15 min. Hydrometer readings were then made at 40 s and at 2 h to obtain total suspended solids. Samples were then poured into a series of sieves dividing the sample into the five sand particle sizes and clay plus silt[23]: very coarse sand, coarse sand, medium sand, fine sand, and very fine sand with the range in sizes being 1−2, 0.5−1, 0.25−0.5, 0.10−0.25, 0.05−0.10, and < 0.05 mm, respectively. The samples were placed in a forced-draft drying oven at 105 °C until a constant weight was reached, then dry-sieved using a Ro-Tap® Shaker utilizing the same size classifications. Soil samples were dried and weighed prior to sieving and each sand fraction was weighed post sieving.

      Bulk density was measured following standard procedures[22] adjacent to plot center where the other soil samples were collected using a core sampler with 48.25 mm diameter rings, and samples oven-dried at 105 °C until constant weight was achieved and weighed prior and after drying. Field capacity and wilting coefficient were measured using a soil pressure plate apparatus and chambers. Field moist samples were soaked in water for 24 h prior to being placed under the pressure plates at both −31 and −1,500 kPa. Subsamples were weighed moist and then oven-dried at 105 °C to constant weight and then reweighed.

      Standard lab methods using an ICP Thermoscientific Analyzer were performed to obtain phosphorus, potassium, calcium, magnesium, nitrogen, organic carbon, and ammonium at the Stephen F. Austin State University Plant, Soil, and Water Laboratory. To obtain pH, a one to two ratio of soil to water using 12.5 g of soil and 25 ml of deionized water method was determined using a pH probe. Electrical conductivity was taken following the completion of the pH using the same prepared sample using an E.C. meter.

      Statistical analysis

    • Analysis of variance (ANOVA) using Proc GLM (General Linear Model) procedure in SAS was used to determine significant differences (p = 0.05). If differences were found among variables, Tukey's mean separation test was then used. Because large set of variables inherently have some correlations, principal component analysis (PCA) was used to summarize all of the variables into unrelated variables (PC1, PC2 ...), and important or significant PCAs were selected to perform regression. The number of principal components evaluated was determined by using randomization in PC-ORD. The top 10 composite variables from each significant PCA were selected and used in step wise regression to determine which variables most influenced longleaf pine site index.

    Results

    • The official descriptions for all three series were: are they are deep, well drained to excessively drained, with some variations in texture, color, and depth of each horizon[24]. Depth to the first argillic horizon ranged from 23 to 49 cm (mean = 42.5 cm) in the Stringtown series, from 55 to 88 cm (mean = 67.1 cm) in the Letney series, and from 101 to 155 cm (mean = 111 cm) in the Tehran series. The greatest difference is depth to the first argillic (Btl) horizon: Stringtown < 50 cm, Letney 50 to 100 cm, and Tehran Bt1 > 100 cm.

      Longleaf pine site indices

    • ANOVA indicated significant differences for longleaf pine site index (Table 1). Mean site indices for Letney and Stringtown soils were within the USDA Natural Resources Conservation Service (NRCS) range of site indices, but was below for Tehran soils (Table 2).

      Table 1.  Means, standard deviations, and coefficient of variations for site index (base age 50) for natural longleaf pine stands on three soil series in east Texas.

      Soil seriesnSite index
      (m)      		Standard deviationCoefficient of variation
      Stringtown1022.2a2.3510.579
      Letney1022.6a1.285.564
      Tehran1020.0b1.607.980
      n = number of plots. Same letter within a column indicates no significant difference (p = 0.05).

      Table 2.  Mean, low and high site index values (base age 50) by USDA-NRCS for Stringtown, Letney, and Tehran soils.

      Soil seriesMean site index (m)Low site index (m)High site index (m)
      Stringtown24.520.726.5
      Letney24.821.332.0
      Tehran26.224.130.8
      n = number of plots.

      Soil physical parameters

    • Within the unweighted soil physical parameters, 12 were significantly different (Table 3). Both depth of A and depth to E on Tehran soils were significantly deeper than Stringtown. As expected, depth to B was significantly different, with Tehran being the deepest and Stringtown being the shallowest. Depth of E was also found to be significantly different, with Tehran being greater than both Stringtown and Letney; depth of B was also significantly different, with Stringtown approximately 73 cm thicker than Tehran, and 31 cm thicker than Letney. Wilting coefficient of the A horizon showed significant differences between Stringtown and Letney soils, with 50% more water held in the Stringtown series (Table 4). B horizon wilting coefficient was significantly greater in Stringtown than the Tehran soils.

      Table 3.  Significant (p = 0.05) soil physical parameters not weighted by horizon thickness, means, and p-values.

      HorizonVariableStringtownLetneyTehranp-value
      AThickness (cm)14.80a19.23ab25.35b0.008
      WC (g·cm−3)0.09a0.06b0.07ab0.026
      MS (%)28.74a31.77ab40.59b0.049
      EDepth to E (cm)14.80a20.53ab25.15b0.010
      Thickness (cm)24.20a49.17b86.45c<0.001
      BDepth to B38.90a70.40b111.80c<0.001
      Thickness (cm)111.10a79.60b38.80c<0.001
      MS (%)20.99a23.21a36.14b0.003
      Silt + Clay (%)46.11a35.46ab29.60b0.014
      Sand (%)64.92a72.58ab78.27b0.004
      Clay (%)26.95a18.45ab13.53a0.013
      Same letter within a row indicates no significant difference (p = 0.05). WC = Wilting Coefficient, MS = Medium Sand.

      Table 4.  Significant (p = 0.05) soil physical parameters weighted by horizon thickness with p-values.

      HorizonVariableStringtownLetneyTehranp-value
      AWC (g·cm−3)1.291.211.680.090
      OM (g·cm−3)0.04a0.05ab0.07b0.005
      EFC (g·cm−3)3.32a6.25a18.94b0.006
      AWC (g·cm−3)2.134.4614.260.019
      OM (g·cm−3)0.06a0.11a0.19b<0.001
      BFC (g·cm−3)36.05a22.24b10.59b0.001
      WC (g·cm−3)26.70a12.14b2.99b<0.001
      AWC (g·cm−3)0.32a0.23ab0.13b0.012
      Same letter within a row indicates no significant difference (p = 0.05). FC = Field Capacity; WC = Wilting Coefficient, AWC = Available Water Capacity, OM = Organic Matter.

      Medium sand in the A and B horizons had the highest percent by weight in Tehran soils over the other soils. Medium sand in the B horizon and wilting coefficient of the B horizon were inversely correlated; as medium sand increased, wilting coefficient decreased. As the depth to the first argillic B horizon increased, both total silt + clay and total clay in the B horizon decreased.

      Six physical variables weighted by horizon thickness were determined to be significantly different by soil series (Table 4). Field capacity in the E horizon was higher in Tehran soils than the others. Stringtown soils were significantly different from Letney and Tehran soils for field capacity and wilting coefficients weighted by thickness of the B horizon, and Stringtown soils held more moisture at field capacity and at wilting coefficient in the B horizon than Letney and Tehran. A and E horizon organic matter content was highest in Tehran. Organic matter content in the B horizon had the opposite trend, where Stringtown soils were significantly greater than Tehran soils.

      Soil chemical parameters

    • Of the 36 soil chemical parameters not weighted by horizon thickness, exchangeable Ca in the A horizon was the only parameters found to be significantly different; Ca concentration in the A horizon in the Letney soils was significantly higher than in the other two soils.

      Weighted by horizon thickness, 17 variables were significantly different (Table 5). Ca weighted by E horizon thickness was not significantly different, but were in the A and B horizons. Organic C in the A horizon was greater in Tehran than in Stringtown soils; and in the E horizon was greater than in Stringtown and Letney soils. The B horizon had the opposite effect, as Stringtown soils contained more organic C than Tehran. Overall, Stringtown contained more total N than Tehran soils, while in the E horizon Tehran soils had more total N; Stringtown had more total N in the B horizon than Letney soils, which had more than Tehran soils. Tehran had more NH4 in the E horizon, but Stringtown had more in the B horizon than Tehran soils. Tehran had more P in the A and E horizons than Stringtown soils, and more K in the E horizon than Stringtown; Stringtown and Letney soils contained more K in the B horizon than Tehran soils. Stringtown soils contained more Mg in the B horizon than Tehran, and Stringtown soils contained more S in the B horizon than Tehran.

      Table 5.  Significant (p = 0.05) mean chemical parameters (mg·Kg−1) by horizon thickness by soil series.

      HorizonVariableStringtownLetneyTehranp-value
      ATotal N19.36a24.60ab34.64b0.0034
      P0.03a0.10b0.09ab0.0308
      K0.26a0.53ab0.60b0.0345
      Ca2.41a6.394.35ab0.0444
      C181.17a269.76ab361.33b0.0054
      ETotal N42.84a82.11b152.10c<0.0001
      NH40.10a0.15a0.36b<0.0001
      P0.05a0.09ab0.14b0.0042
      K0.87a1.30ab1.70b0.0461
      C292.42a534.65a959.73b<0.0001
      BTotal N217.14a164.65b8.74c<0.0001
      NH40.48a0.32ab0.22b0.0105
      K4.75a4.86a1.51b0.0062
      Ca67.87a65.92a14.35b0.0026
      Mg21.21a15.04ab2.55b0.0056
      S2.62a1.54ab0.50b0.0254
      B0.010.010.000.0753
      C1577.93a1125.85ab669.62b0.0118
      Same letter within a row indicates no significant difference.

      Generally, Stringtown had higher concentrations of nutrients in the B horizon than Tehran soils, although Tehran had higher concentrations in the A and E horizons. Within the A horizon, clay content was highest in the Letney soils which would provide a higher cation exchange capacity. K and Ca within the A horizon which were higher in Tehran and Letney soils; Stringtown averaged lower silt and clay in the A horizon resulting in lower quantities of those nutrients within the A horizon. Total N was highest in the A horizon in the Tehran which also contained the most organic C.

      Soil profile nutrients were weighted by horizon depth and then summed for the entire 150 cm soil profile; Ca, Mg, and S were significantly different (Table 6). Stringtown and Letney soils contained more Ca than Tehran, and Stringtown soils contained more total Mg and S than Tehran. Soils with argillic B horizons closer to the surface (Stringtown and Letney) tended to have higher total available nutrient contents than Tehran. Total amounts of Ca, Mg, and S were found to be greatest in the Stringtown soils; Stringtown had the thickest B horizon relative to the 150 cm profile depth, and also had the highest amounts of silt and clay.

      Table 6.  Significant mean (g) chemical parameters within the 150 cm soil profiles with means (g) by soil series.

      VariableStringtownLetneyTehranp-value
      Ca79.85a86.93a32.90b0.0040
      Mg23.62a18.01ab5.53b0.0048
      S3.11a1.86ab0.880.0174
      Same letter within a row indicates no significant difference (p = 0.05).

      Principal Component Analysis ordination

    • Five variable combinations accounted for approximately 62% of the variation (Table 7) using principal component analysis. PC1 (21% of the variance) was strongly driven by depth to the B horizon, thickness of the B and E horizons, percent silt and clay in the B horizon, total wilting coefficient of the B horizon and entire profile, and total organic matter in the E horizon. PC2 (15% of the variance) was driven by percent medium sand, total sand and silt in the A horizon as well as percent medium sand, total sand, and silt in the E horizon. PC3 (10% of the variance) was driven by field capacity and available water capacity of the A and B horizons, total potential field capacity and available water capacity of the A horizon, total potential available water capacity of the B horizon, and total potential available water capacity for the profile. PC4 (7% of the variance) was driven by field capacity, wilting coefficient, and available water capacity of the E horizon and total field capacity of the entire profile, while PC5 (7% of the variance) was driven by percentage of very coarse sand, coarse and medium sand in the A horizon, percentage of very coarse sand and medium sand in the E horizon, and percentage of very coarse sand, coarse sand, and total clay in the B horizon.

      Table 7.  Results with p-values from each of the first 10 principal components from 999 randomizations to determine significant components based on relationship to the maximum theoretical eigenvalue vs the true eigenvalue for all physical variables, chemical variables and physical and chemical variables combined with associated % variance.

      AxisEigenvalueMaximum Eigenvalue% of VariationCumulative variationp-value
      Physical parameters
      113.097.46720.77920.7790.001
      29.6835.82915.37136.1500.001
      36.5855.18710.45246.6020.001
      44.8924.8957.76554.3670.002
      54.5194.5407.17361.5400.002
      63.2434.0755.14766.6871.000
      72.9993.7514.76071.4471.000
      82.5255.5324.00875.4551.000
      92.2123.2943.51178.9661.000
      102.0823.0503.30582.2711.000
      Chemical parameters
      116.2167.43923.50123.5010.001
      211.5605.98616.75340.2540.001
      39.7575.48814.14054.3940.001
      45.8225.1128.43862.8320.001
      54.2484.7306.15668.9880.2.92
      63.2694.3784.73873.7261.000
      72.6494.0743.83877.5641.000
      82.4503.7393.55181.1151.000
      91.8903.5003.50083.8541.000
      101.7323.4012.51086.3641.000
      Combined parameters
      125.10410.93919.01819.0180.001
      215.5019.35311.74330.7620.001
      313.8218.66810.47041.2320.001
      411.7918.1978.93350.1650.001
      58.0317.6676.08456.2490.001
      67.5957.2865.75462.0030.001
      76.9896.8325.29567.2980.001
      85.6846.4714.30671.6040.983
      95.0906.2003.85675.4601.000
      104.0865.9333.09678.5561.000

      Four significant PCA’s accounted for approximately 63% of the variation among the soil chemical variables (Table 7). PC1 (24% of variance) were concentrations of K, Ca, Mg, and boron in the B horizon, as well as total Mg and boron weighted by depth of the B horizon, and total K, Ca, Mg, and S weighted by depth of the 150 cm soil profiles. PC2 (17% of the variance) was driven by concentration of K, Ca, Mg, S, and Boron in the E horizon, total K, Ca, Mg, S, and B weighted by depth in the E horizon. PC3 (14% of the variance) was driven by total C, P, K, Ca, and Mg weighted by depth in the A horizon, as well as total grams of P weighted by depth of E horizon and total NH4+ and total N weighted by depth of the B horizon. PC4 (8% of the variance) was driven by total C and P within the entire profile, and total N, P, and C in the B horizon.

      Seven variables accounted for 67% of the variation (Table 7) when the physical and chemical variables were combined for analysis. PC1 (19% of the variance) was driven by depth to B, thickness of E and B, wilting coefficient of the B horizon, percentage of silt and clay in the B, total potential wilting point of the B, total wilting point in the profile, organic matter in the E horizon, Mg in the B, total N, K, Ca, Mg, S, and B in the B horizon, and total K, Ca, Mg, and S. PC2 (12% of the variance) was driven by field capacity and available water capacity in the A horizon, field capacity and total field capacity of the B horizon, total field capacity in the profile, P, K, Ca, and Mg in the A horizon, NH4 in the E horizon, NH4 in the B horizon, total P, K, Ca, Mg and S in the A horizon, and total NH4 in the profile. PC3 (10% of the variance) was driven by Ca, Mg, and B in the E horizon, Mg, S, and B in the B horizon, total grams of P in the A horizon, and total Mg in the E, and total B in the B horizon. PC4 (9% of the variance) was driven by percent coarse sand, very fine sand, silt and clay, and total sand in the A horizon, percent coarse sand, fine sand, very fine sand, total sand, and total silt in the E horizon, total organic matter in the profile, P in the B horizon, and total C, P, and B in the profile. PC5 (6% of the variance) was driven by concentrations of K, Ca, Mg, S, and B of the E horizon and concentration of C in the B horizon. PC6 (6% of the variance) was driven by bulk density and organic matter in the E and B horizons, concentration of B in the A horizon, concentration of C in the E horizon, concentration of P and C in the B horizon, and total grams of NH4, and B within the profile, while PC7 (5% of the variance) was driven by clay, wilting point, and total potential wilting point of the E horizon, concentration of Ca in the A horizon, and total Ca in the A horizon.

      Regression for site index

    • Seven variables were the most significant soil physical factors affecting longleaf pine: depth to B, thickness of the E and B horizons, percent silt and clay in the B horizon, wilting coefficient of for the B horizon, wilting coefficient of the profile, and percent organic matter in the E horizon.

      The best two-variable model (1) included depth to the B horizon and total wilting coefficient of the B horizon (R2 of 0.3984):

      Siteindex=88.71063−(Depth(cm)toBhorizon∗0.19074)−(TotalBhorizonwiltingpotential∗0.26955) (1)

      Ten soil chemical variables correlated most with site index of longleaf pine: total K, Ca, Mg, and S in the profile, total Mg and B in the B horizon, and concentrations of K, Ca, Mg, S, and B in the B horizon. Only a one variable model (2) best fit the site index (R2 of 0.2026):

      SiteIndex=66.93652+(TotalCa(mg)inProfile∗0.05947) (2)

      Combining all variables, the variables most correlated to site index were depth to the B horizon, wilting coefficient of the B horizon, percent silt and clay and depth weight wilting coefficient of the B horizon, the profile weight wilting point of the whole profile, organic matter of the E horizon, concentration of Mg in the B horizon, total N, K, Ca, Mg, S, and boron in the B horizon, and total K, Ca, Mg, and S in the B horizon.

      Using step-wise regression, the top variables that affect longleaf pine site index were total N and S in the B horizon, concentration of Mg in the B horizon, total Mg and S in the profile, and wilting coefficient weighted by horizon thickness in the B horizon. These six-variables proved to be the best model (R2 = 0.6668). Regression Eqn (3) for site index using these six variables was:

      SiteIndex=64.98+(TotalN(mg)inB0.05119)+(TotalMg(mg)inprofile1.66002)+(TotalS(mg)intheBhorizon5.87648)−(concentrationofMg(mg⋅cm−3)inBhorizon0.22445)−(TotalS(mg)inprofile5.25599)−(TotalwiltingpotentialinBhorizon0.53062) (3)
    Discussion

      Soil physical and chemical parameters

    • The wilting coefficient was affected by the amount of silt and clay within the profile; as silt and clay decreased, so did the wilting coefficient. Within the unweighted soil physical parameters, field capacity and available water capacity should be affected by this texture correlation; however, it was not found in this study. In fact, available water capacity was highest in the deeper sand soils, suggesting that the pressure plate method used in this study may not have produced reasonable results.

      The A horizon, as expected, contained more organic matter than either of the other two horizons (Table 3). The E horizon is characterized as where leaching of humus, silt and clay, and various ions occur, while the B horizon is where accumulation of humus, silt and clay, and various ions occur. In this study the B horizon did contain higher percentages of silt and clay than did the A and E horizons. Within all soils, as depth to the first argillic B horizon increased, the percentage of silt and clay decreased in the B horizon. Conversely, sand increased as depth to the first argillic horizon increased.

      Similar to a previous study[25], the wilting coefficient was influenced by the proportion of B horizon in the entire profile, which had higher wilting coefficients in all soils. Texture and B horizon thickness played a big role due to the inherent ability of fine textured soils to hold more water at the wilting coefficient. However, neither field capacity or available water capacity were statistically different between soils, likely due to the pressure plate system retaining more water than it should have. The B horizon total potential field capacity was significantly highest in the Stringtown soils, likely a result of the increase in silt and clay content in the thicker B horizon in those soils.

      The higher percentage of clay may have reduced Ca leaching to the lower profile. The presence of finer texture soils increases cation exchange capacity (CEC), retaining cations. It is unclear why Stringtown had lower concentrations of Ca than Letney and Tehran, or why concentrations of other nutrients were not significantly different. The E horizon total potential field capacity was highest in the Tehran soils, again indicating that the pressure plate method did not produce reasonable results. Finer texture soils have higher CEC, which can result in the presence or ability to hold more cations[26].

      As depth to the B horizon increased, clay content decreased, as did water holding capacity, available water capacity, and wilting coefficient. Texture is inherently related to the amount of water a soil can hold. All three of the soil series in this study had A, E, and B horizons. In this study, texture did not prove to be as important as expected.

      Many studies have looked at the relationship between the depth of horizons and site index, with varying results. While no correlation between site index and depth to the first argillic horizon B for longleaf pine was found in east Texas[27], a negative correlation between site index and depth to finer textured layers in sandy soils for white oak (Quercus alba) was found in Michigan[28], but the influence began to wane at depths greater than 1.5 m. Site index of radiata pine (Pinus radiata) increased with increasing depth of the topsoil[29], but evaluated total soil depths in Douglas-fir did not show any significance with soil depths ranging from 50 cm to 100 cm[30]. No single physical variable had a well-defined correlation with site index for shortleaf pine[31]; however, depth to the B horizon and texture of that horizon was found to be a good indicator of shortleaf pine on the soils studied, is similar to what our study discovered. Higher soil organic matter resulted in increased site quality for most soils[25, 32] and soil texture was not the only parameter affecting water holding capacity: other factors included OM and soil bulk density along with gravel content and salinity, which also affected water availability.

      Ca, N, P, K, and Mg did not affect site index when added to radiata pine[29]. This partially conflicts with our results, as Ca was the only nutrient that was significant. Growth is often limited in many forests in the southern United States by nutrient availability, promoting fertilization in silvicultural practices. Nitrogen and phosphorous are often considered the most common nutrients that limit growth in southern pine forests[33]. Similarly, nutrients had a positive correlation with site index in radiata pine[29], but they did not specify a given depth at which the nutrients were most effective. Our study found that nutrient levels in the B horizon had a strong correlation with site index increase in longleaf pine as did total amounts of nutrients in the profile, which were usually correlated to soils with higher silt and clay concentrations and shallower B horizons. Boron has been shown to positively correlate with growth in pine[34]; our study also showed boron being positively correlated to longleaf pine site index, but did not differ between soils. Principal component analysis found depth to each horizon as well as certain nutrients within the B and E horizons as important. Higher nutrient levels correlated positively with site index[29], but they did not specify if this relationship was by horizon or total in the profile. Our study showed low total nutrients within all three soils, and as depth to the first argillic B horizon increased, nutrients available decreased. Stringtown appeared to have greater amounts of nutrients within the B than the Tehran in most situations with Letney soils intermediate in nutrient availability.

      Total N, Mg, and S in the B horizon had a positive effect on longleaf pine site index, while concentrations of Mg and S had a negative impact. Indirectly, depth to the B horizon had a correlation to these values as the thickness of the B impacts total available nutrients in that horizon. In addition, as B horizon depth increased, the percent clay within the horizon decreased. The presence of sand in the A and E horizon resulted in a lower cation exchange capacity, allowing for nutrients to leach through these horizons while the B horizon had an increase in clay.

      While all of the three soils found on sites in this study historically supported longleaf pine, these results highlight how site index for this species might be driven by variables other than a sandy, well-drained A horizon. The importance of the B horizon depth and the lower ability of the B horizon to allow water to drain was an important variable found in the PCA analysis. Two of the soils in this study had site indices within the NRCS range for longleaf pine, but both were lower than the mean for those soils. The Tehran soil was slightly below the minimum site index for that soil. It could be on that soils, the NRCS underestimated the importance of the depth of the B horizon. In addition, longleaf pine historically on those soils may have benefited more from the short-interval fire frequency than on the other two soils.

    Conclusions

    • Soil physical parameters in the A and E horizons did not appear to greatly influence site index for longleaf pine on these soils in east Texas. However, the depth to B and wilting coefficient of the B influenced site index of longleaf pine on these three soils, which suggest that water availability may play the largest role in affecting site index on these deep, coarse textured soils. Soil chemical parameters in the A and E horizons did not appear significant; however, soil chemical parameters in the first argillic B horizon, as well as nutrient availability in the whole profile did. Soil variables in the B horizon affect site index for longleaf pine the most, while some variables within the whole 150 cm profile also had an effect. This is likely due to the effect of the weighted by horizon thickness of the B horizon had on the total profile because of clay content of the horizon providing for higher available water content, and nutrient storage. Some A horizon parameters showed some slight effect on longleaf pine site index, but this could possibly be due to the amount of organic matter within the A horizon. Each model highlighting different variables reflects the complexity of the interaction of soil variables with site index. Productive forests tend to have soils with favorable physical properties that enhance biological functions. Separating and choosing the most significant of these soil variables can be challenging due to the inherent complexity and interactions among many of them.

      Future studies should look at rooting depth within each of these three soils as well as the effect of soils with drainage classes that are known to hold more water. Studies should also consider planting on these three-soil series using the same treatments to how these soil variables affect longleaf pine regeneration.

      For those making management decisions on locations with the greatest potential success for longleaf pine establishment, soils with similar A and B horizon characteristics may have the greatest success, along with the use of recurring prescribed fire.

    Author contributions

    • The authors confirm contribution to the paper as follows: study design and manuscript preparation: Oswald BP, Farrish KW, Svehla R; data collection and data analysis: Svehla R. All authors reviewed the results and approved the final version of the manuscript.

    Data availability

    • The datasets generated during and/or analyzed during the current study are available from the corresponding author on reasonable request.

      Acknowledgments

      • The Division of Environmental Science in the Arthur Temple College of Forestry and Agriculture at Stephen F. Austin State University provided financial support for this project. Our sincere appreciation is extended to Megan McCombs, Cassity Aguilar and Jacalyn Jones for their assistance in the field and with the lab analysis.

      Conflict of interest

      • The authors declare that they have no conflict of interest.

      Rights and permissions
      • Copyright: © 2024 by the author(s). Published by Maximum Academic Press, Fayetteville, GA. This article is an open access article distributed under Creative Commons Attribution License (CC BY 4.0), visit https://creativecommons.org/licenses/by/4.0/.
    Table (7) References (34)
  • About this article
    Cite this article
    Oswald BP, Svehla R, Farrish KW. 2024. Soil parameters affecting longleaf pine (Pinus palustris) site quality in east Texas. Forestry Research 4: e002 doi: 10.48130/forres-0023-0031
    Oswald BP, Svehla R, Farrish KW. 2024. Soil parameters affecting longleaf pine (Pinus palustris) site quality in east Texas. Forestry Research 4: e002 doi: 10.48130/forres-0023-0031
    shu

Catalog

  • HTML

  • Introduction
  • Materials and methods

  • Study area

  • Plot establishment

  • Field sampling

  • Laboratory methods

  • Statistical analysis

  • Results

  • Longleaf pine site indices

  • Soil physical parameters

  • Soil chemical parameters

  • Principal Component Analysis ordination

  • Regression for site index

  • Discussion

  • Soil physical and chemical parameters

  • Conclusions
  • Author contributions
  • Data availability
  • Acknowledgments
  • Conflict of interest
  • Rights and permissions
  • About this article
logo

Contents

Policies

Services

Partnerships

© 2025 Maximum Academic Press, except Open Access articles. MAP is a member of ALPSP, COPE, CrossRef, cOAlition S, OASPA and STM.

Privacy Policy  |  Site Map  |  Terms and Conditions

/

DownLoad:  Full-Size Img  PowerPoint
Return
Return