DArTseq molecular markers associated with the piping leaf margin phenotype in pineapple (Ananas comosus L.)

This study sought to understand the genetic basis of the piping leaf margin phenotype in pineapple. To achieve this aim, a genome-wide association study (GWAS) using mixed linear regression and logistic regression analysis was conducted on three pineapple diversity panels including seedling populations segregating for spiny, spiny-tip and piping leaf margins. This study identified single nucleotide polymorphism (SNP) markers associated with the piping and spiny-tip leaf margin phenotypes. A broad quantitative trait locus (QTL) positioned on chromosome 23 between positions 240,475 and 2,369,197 bp was the most highly associated with piping leaf margin in all analyses. Major candidate genes proposed are a Zinc finger protein 2, a Zinc finger protein 3, a WUSCHEL-related homeobox 2, a WUSCHEL-related homeobox 1 and a Zinc finger protein CONSTANS-like. Some other genes of a lower association, linked or nearby genes of interest, are also considered potentially involved to varying degrees. All candidate genes are known to be involved in aspects of stem cell maintenance, cell proliferation, epidermal cell differentiation, organogenesis, leaf polarity, cell wall modification or hormone signalling. It is possible each plays a role in either differentiation or morphological aspects of the spiny-tip and piping leaf margin phenotypes. It is expected the relative role of each associated gene might vary with genetic background.


Introduction
Pineapple exhibits five leaf margin types and all, except complete spiny, are considered the product of domestication (Coppens d'Eeckenbrugge and Sanewski 2011). Genotypes exhibiting this phenotype are typically found in the north-western regions of Brazil (Coppens d'Eeckenbrugge et al. 2018a), the country of origin of Ananas. The piping margin typically exhibits a thin white border along the leaf margins and is devoid of any spines ( Fig. 1A and B). Leaves with a spiny margin are considered the natural state as all wild genotypes possess this phenotype and it is reasonable to assume piping represents a mutation selected at some stage in the estimated 6,000 to 10,000 years of domestication of pineapple (Coppens d'Eeckenbrugge et al. 2018b).
The piping leaf margin does not appear to offer any obvious evolutionary advantage and hence its widespread existence in the western Amazon is likely all attributable to mankind. Because the trait exists in several different genotypes which are considered to have arisen from clonal selection and sexual recombination (Chen et al. 2019), the piping phenotype most likely has been selected multiple times.
The first obvious sign of the piping phenotype occurs about 8 cm behind the tip and extends to the base of the leaf where it emerges from the meristematic zone. The demarcation or boundary between adaxial and abaxial 'domains' appears clear but a narrow strip of abaxial surface has been folded onto the adaxial surface to present a new 'false' margin (Collins 1960). The piping strip will vary in width in different genotypes from almost non-existent to a strip approximately 2 mm wide. The width can also vary on different leaves of the same plant or on different margins of the same leaf possibly revealing an environmental effect on gene penetrance, meristematic cell development or leaf cell turgor. The fold can be exaggerated and obvious or, in some cases, only present as a shallow valley. There are however no published studies describing the botanical nature of piping leaf margin.
Leaf margin spines do not usually occur with piping suggesting an additional genetic mechanism. According to Collins (1960), the piping gene is epistatic to the spiny and spiny-tip alleles. He goes further to say, plants homozygous for piping have a more pronounced piping margin, suggesting a dosage effect and plants homozygous recessive for the piping alleles will display as either spiny or spiny-tip. The piping phenotype appears dominant over spiny and spiny-tip phenotypes. The piping leaf margin appears very stable and there are few reports of reversion to spiny, unlike the spinytip phenotype which frequently reverts to spiny and displays considerable variability in gene penetrance (Sanewski 2020).
This study examined the botanical nature of piping leaf margin and identified SNP markers associated with spiny-tip and piping leaf margin phenotypes. The use of three diversity panels, independently analysed, was a novel approach that enabled testing between pairs of phenotype states; viz, piping and spiny, piping and spiny-tip and spiny and spinytip. Putative candidate genes associated with the piping leaf margin, all located on chromosome 23, are discussed.

Methods
In this study, the use of the term 'piping' is meant to encompass smooth margins with or without folding and/or an obvious white edge. Plants were phenotyped by visually inspecting the leaf margins and carefully feeling the margin to the very tip for the presence of small spines. Those without any spines anywhere on the leaf were considered as piping, those with small spines at the tip as spiny-tip, and plants with spines completely over the margins were classed as spiny.
Histology 10 µm thick transverse leaf margin sections were prepared from cultivars, 'Tapiricanga' (piping) and 'Smooth Cayenne' (spiny-tip), using a freezing microtome to examine the differences between a piping margin and a margin with spine suppression but no piping. The sections were fixed in 4% formalin and mounted in Mowiol ® solution for brightfield microscopy.

Panel Descriptions for GWAS
Three diversity panels were used. Diversity panel one was a sub-set of 90 plants from a population of seedlings exhibiting piping (33 plants) or completely spiny leaf margins (57 plants). This panel was intended to reveal the genetic difference between spiny and piping margin. The population Fig. 1 A and B. The piping leaf margin typically exhibits a thin white border along the leaf margin (1A) and is devoid of spines even at the tip. On some leaves of some genotypes, it might also be characterised by an upward fold in the margin (1B) was developed using the parents that included spiny, spinytip and piping genotypes. Sixty-two of these plants were of the parent combination 'PRI 59-656' (pollen) and 'MD-2' (seed). It is considered that the piping margin trait in all progenies used was inherited from the piping margin 'Smooth Guatemala' component and the spiny trait was derived from 'MD-2'.
Panel two comprised 516 genotypes with 12 full sib families and some of their parental varieties. There were 440 plants with piping leaf margins and 76 with spiny-tip. The population had all derived from varieties from the Pineapple Research Institute (PRI) of Hawaiian breeding program as for panel one. Panel two differed from panel one in that it represented a larger number of recombinant events and spiny leaf margin plants were not included. This panel was intended to reveal the genetic difference between spiny-tip and piping margin genomes.
Diversity panel three was an amalgamation of panel two seedlings (508 seedlings) and a reciprocal bi-parental population (367 seedlings) ('MD2' × 'Cayenne', and 'Cayenne' × 'MD2') of seedlings segregating for spiny and spinytip margin. Panel three also included 76 parental genotypes and several unrelated genotypes including wild accessions. Panel three was the most diverse of the populations used. There were 492 plants with piping margins, 351 with spinytip and 108 with spiny margins to give a total of 951 plants.

Molecular Marker Development
Freeze-dried leaf base samples (white tissue) were supplied to Diversity Arrays Pty Ltd for DArTseq marker development. Marker development followed the procedure described in Sanewski et al. (2017). The methylation-sensitive restriction enzymes PstI and MseI were used to avoid highly repetitive regions that bias GWAS, and to reduce the genome complexity. The marker sequences were trimmed to 69 bp.
Subsequent to the sequencing of panel one genotypes, a revised version of this genome sequence became available and was used for panels two and three. Marker sequences for panel two and panel three genotypes were referenced against the revised genome of the spiny-tip pineapple cultivar 'Smooth Cayenne' clone 153 (Acomosus_321_v7.1, assembly GCA_902162155.1 at European Nucleotide Archive; https:// www. ebi. ac. uk/ ena/ brows er/ home). This version of the genome is slightly different from the previous version. Chromosomes 24 and 25 have been amalgamated as a larger chromosome 25. Chromosome one was divided and the later portion re-labelled as chromosome 24. This was done by the authors of that work to correct for an error in the original assembly (Chen et al. 2019). These two reference genomes should be identical for chromosomes two to 23 and hence the different panel results are comparable where markers of interest are positioned only on chromosomes two to 23.
Before analysis, the SNP data sets were filtered on reproducibility (> 0.99), call rate (> 0.75) and MAF (> 0.025) before imputation using the LD KNNi approach in TASSEL 5.2.44 (Bradbury et al. 2007;Money et al. 2015). The accuracy of imputation was tested by masking 10% of SNPs and comparing the imputed file with the original.

Principal Component Analysis (PCoA)
Principal component covariates were generated in TASSEL 5.2.54 (Bradbury et al. 2007). The number of covariates for GWAS was determined based on either the Bayesian Information Criterion (BIC) approach (Schwarz 1978) performed in GAPIT, or the QQ and Scree plots. The smallest number of covariates that produced a QQ plot that was neither inflated or under-inflated and had the smallest P values and was consistent with the Scree plot (visualized in TASSEL but not reproduced) was used in the GWAS.

Association Analyses
The SNP marker data was analysed by genome-wide association study methods (GWAS) using mixed linear models (MLM) with structure correction. The multi-locus mixed linear models (MMLM) of Segura et al. (2012), FarmCPU ) and BLINK (Huang et al. 2019) and the single loci ECMLM model, all in the GAPIT R package version 3 (Zhang et al. 2010;Tang et al. 2016) were used. Kinship was estimated using the Van Raden method in GAPIT (Van Raden 2008).
The phenotype data was binary and coded as 0 and 1 for panels one and two and 0, 1 and 2 for panel three. The Bonferroni-adjusted P value of P ≤ 0.01 was used as the level of significance.
The R package Bayesian Genomic Linear Regression (BGLR) was used, through the graphic user interface, Intelligent Prediction and Association Tool (iPat v1.3) (Chen and Zhang 2017) to test SNP markers in panels one, two and three. The Bayes A model was used with a probit link to optimise analysis of binary or ordinal data through logistic regression. The algorithm was run with 50,000 iterations and the first 3,000 discarded. The Van Raden kinship matrix and principal components were used to adjust for relatedness. The marker effects posterior values were converted to absolute values and the Z-score calculated by dividing the absolute marker effect by its standard deviation. The P value was then calculated by the following equation according to Campbell et al. (2019), P = 2*(1-pnorm(Z-Score)). A Bonferroni-adjusted P ≤ 0.05 was used as the cut-off for significance.
A mixed model logistic regression was also performed using the R program Generalized Linear Mixed Model Association Test (GMMAT v1.1.2) (Chen et al. 2016) for panels one and two data. This model uses a Bernoulli distribution and logit function with covariate adjustment by principal components as a fixed effect and kinship (Van Raden) as a random effect to adjust for population structure. Logistic models such as GMMAT that can include fixed and random effects to account for population structure are considered more suitable for binary traits than linear mixed models especially if some populations have more or less 'affected' status rather than 'control' status. In such situations there is an increased risk of type I errors with mixed linear models as they assume error terms are normally distributed (Shenstone et al. 2018). The model 10 algorithm in GMMAT was used for to perform score tests on marker associations. The data were treated as binomial.
Logistic regression was also performed in gPLINK using a logit function (Purcell et al. 2007) for panel one only.

Linkage Analysis
Linkage disequilibrium (LD) for the linkage decay curve was estimated for the entire genome for panel one data using the software package TASSEL 5.2.44 using a sliding window of 100 SNPs. The r 2 of SNP marker pair frequencies and their distance was used to generate a linkage decay curve. A Loess curve was fitted using the package SigmaPlot 14. Linkage disequilibrium decay distance was calculated as the point on the x axis corresponding to where the Loess curve passed 50% of the maximum r 2 (Noble et al. 2018).
Marker linkages were also analysed for individual chromosomes for panel one SNP markers using Haploview 4.2 (Barrett et al. 2005). In the case of Haploview, the Hardy-Weinberg p value cut-off was set at zero.

Nearby Genes
Genes within the LD decay distance were identified using the annotated 'Smooth Cayenne' v3 genome for panel one genotypes and v7 genome for panels two and three in the software package 'Geneious R10'. This differential use of genome versions was because these different genomes were used as the reference for marker positioning in those marker datasets. Likely candidate genes were selected based on proximity and supporting published reports of gene/protein biological function.

Histology
The section for spiny-tip cv. 'Smooth Cayenne' (Fig. S1) shows the co-ordinated adaxial and abaxial layers. The abaxial layer is characterized by a corrugated epidermal surface which appears slightly thinner than the adaxial surface. The transition between the adaxial and abaxial layers dissects the vascular bundles such that the phloem bundle and its sclerenchyma sheath are in the abaxial layer, and the xylem bundle and its sclerenchyma sheath are in the adaxial layer. Alternating between the phloem-xylem bundles, but only in the abaxial layer, are aeration canals. Both the abaxial and adaxial layers have a line of fibre bundles relatively close to the epidermal layer. Both the adaxial and abaxial layers converge to the leaf margin.
The piping leaf margin section shows the abaxial piping 'bulge' consists of one to two repeats of the phloem vascular bundle and fibre strands in two different directions ( Fig. S2A and B). The piping section therefore appears to be one or two repeats of the abaxial layer at the margin sandwiched on top of each other before the terminating epidermal cell layer.

Molecular Marker Datasets
The original seedling SNP dataset for panel one comprised 3,574 markers with 18% heterozygosity and an average minor allele frequency (MAF) of 0.23. The average PIC was 0.31, one ratio 0.47 and call rate 0.85. There were 13.5% null alleles before imputation. Imputation accuracy was 91% and there were 2,141 SNP markers with an average MAF of 0.23 and 0.18 heterozygosity after imputation. 77 genotypes remained for association analyses. The original SNP dataset for panel two comprised 516 genotypes and 40,686 markers and contained 4.8% null alleles. After imputation, with an accuracy of 88%, 516 genotypes and 16,393 markers remained with heterozygosity of 0.32, average MAF of 0.29 and average PIC of 0.18. After filtering for MAF ≤ 0.05 and imputation (accuracy 0.88), the panel three dataset included 16,393 SNP markers for 952 genotypes and had an average MAF of 0.29 and heterozygosity of 0.33.

Principal Component Analysis (PCoA)
Based on the PCoA analysis, zero to four principal components (PCs) were used as covariates for panel one depending on the algorithm, up to nine for panel two depending on the algorithm and three for panel three. For panel one, the first four PCs accounted for 27% of the genotypic variance. For panel two, the first 9 PCs accounted for 37% and for panel three, the first three PCs accounted for 18%.

Marker Density
Density for SNP markers for chromosome 23 using the panel three dataset after filtering and imputation is shown in Fig. 2. Figure 2 showed that marker density was high on chromosome 23 at the starting end up to mid-way around position 4,504,232 bp. Chromosome 23 appears to be telocentric.

Linkage Disequilibrium
LD decay distance for the entire genome using panel one SNPs indicated an effective decay distance of 1,222 kb for the entire genome (Fig. S3). This LD decay distance was used as a guide for selecting candidate genes in all analyses. For panel one, most of the significant markers were intergenic and linkages were generally not strong between associated SNPs.

Association Analyses
Panel one comprised piping and spiny phenotypes only. The highest associated SNP marker (#4714214) was the most significantly associated in all panel one mixed linear model (MLM) analyses (Table 1). Genes close to this position include a Zinc-finger protein 2 (ZFP2), Zinc-finger protein 3 (ZFP3), Ras-related protein Rab-8B, and WUSCHELrelated homeobox 2 (WOX2) genes (Table S1). The WOX2 and ZFP2 are considered the main candidates associated with piping and spiny leaf phenotypes in this panel of genotypes. SNP marker #4714214 was homozygous for the reference allele (spiny-tip 'Smooth Cayenne') in spiny genotypes and homozygous for the SNP or heterozygous in piping genotypes. SNP marker #4713909 was homozygous for the SNP allele or heterozygous in spiny genotypes and homozygous for the reference allele in piping genotypes. Several other markers on chromosomes one and two were also significantly associated but most were intergenic and not matched to genes. Several other genes of interest were however positioned close to these markers. The Manhattan and QQ plots for GWAS of panel one SNPs are shown in Figs. 3A and 4A. Figure 5 shows the significantly associated SNP markers by ECMLM and gPlink for the main QTL on chromosome 23. SNP marker #4714214, at position 798,190 bp on chromosome 23, was the most highly associated SNP marker by GMMAT and gPlink. SNP marker #4713909 was the only significantly associated SNP marker using Bayes A in BGLR (Table 1). This intergenic marker is positioned on chromosome 23 at 762,438 bp and was close to the most significantly associated SNP marker (#4714214). Both #4714214 and #4713909 are positioned close to a WOX2 (Aco015463.1) and ZFP2 (Aco015474.1) (Table S1). Linkage data for the most highly associated SNP markers in panel one (spiny and piping margins) is shown in table S2. The highest associated SNP marker by BGLR, #4713909, was close to the WOX2 and ZFP2 but also highly linked to #100067496 which was also close to the WOX2 and ZFP2. There were no other strong linkages with associated SNP markers.
Panel two comprised spiny-tip and piping margin plants. SNP marker #54316030 was the most significantly associated in all panel two GWAS analyses and logistic regression by BGLR and GMMAT (Table 2). This marker was positioned on chromosome 23 at 1,068,662 bp, 9.3 kb from a WUSCHEL-related homeobox 1 (WOX1) (Aco007021.1) and 31.6 kb from a Zinc finger CONSTANS-like protein  (Table S3). These are considered the best candidates for the piping and spiny-tip genotype data set studied. This marker was homozygous or heterozygous for the polymorphism in piping genotypes. Spiny-tip genotypes were homozygous for the reference allele (spiny-tip cv. 'Smooth Cayenne'). Marker #28878068 positioned at 671,978 bp on chromosome 23 was highly associated by MLM only. This marker was unmatched but positioned 6.9 kb from a Myb domain protein 97 (Aco015454.1), 28 kb from an O-fucosyltransferase family protein (Aco015458.1) and 66.8 kb from a WUSCHEL-related homeobox 2 (WOX2) (Aco015463.1). Piping genotypes were homozygous for the polymorphism and null or homozygous for the reference allele (spiny-tip cv. 'Smooth Cayenne') in spiny-tip cultivars. Other markers of relatively lower significance were also demonstrated by some algorithms on chromosome 5, 6, 12, and 13 and un-positioned. The Manhattan and QQ plots for GWAS of panel two SNP data are shown in Figs. 3B and 4B.
Panel three comprised spiny, spiny-tip and piping margin phenotypes. Marker # 54316030 at position 1,068,662 on chromosome 23 was the most significantly associated in three MLM algorithms (Table 3). This marker was only 9.3 kb from a WOX1 (Aco007021.1) and 31.6 kb from a Zinc finger CONSTANS-like protein (Aco007020.1) (Table S4). These are the same genes proposed for panel two. Marker # 4723993 on chromosome six at 14,330,763 bp was also highly associated by three    (Table S4).

Discussion
There do not appear to be any global disturbances in leaf polarity. The development of a piping margin appears consistent with a cessation of the adaxial layer while the abaxial layer is repeated an additional one to two times at the margin thus creating a thicker composite abaxial layer, an imbalance and upward folding. Piping leaf margin appears to be a lack of co-ordination between development of the abaxial and adaxial tissues at the leaf margin rather than changes in adaxial-abaxial patterning.

Diversity Panel One
Panel one comprised a population segregating for spiny and piping margin and as such is considered the definitive data set for examining polymorphisms associated with the selection of piping-leafed variants from the presumed original spiny genotypes. The assumption at the commencement of this study is that the piping phenotype arose from the spiny phenotype. Discussion of candidate genes within this context is therefore based primarily on the findings from panel one.

WUSCHEL-Related Homeobox 2 (WOX2)
The WOX2 (Aco015463.1) gene is 59 kb from the highest associated SNP by all four MLM models and two logistic regressions and hence considered a highly likely principal causal variant, along with the ZFP2. There are ten WOX genes in pineapple, spread across six chromosomes with three on chromosome 23 (Rahman et al. 2017). Research on WOX genes in other species supports a potential role in leaf margin morphology. Rice plants carrying a mutation in a WOX3 called LSY1, exhibited curled leaf margins suggesting greater cell growth or proliferation in the abaxial leaf tissues compared with the adaxial side. In another example, a WUS ortholog in rice is expressed mostly in leaf margins (Somssich et al. 2016). Also, overexpression of a WOX3 in rice results in adaxial curling of leaf margins. Similarly, a WOX1 double mutant in rice exhibits reduced leaf blades and thickened leaf margins (Honda et al. 2018). WOX1 promotes cell proliferation (Tadege 2016) suggesting a loss or reduced-function polymorphism would reduce cell proliferation.

Zinger Finger Proteins (ZFPs)
The highest associated SNP marker was 1.8 kb from a ZFP2 (Aco015474.1). ZFPs form a large family of transcription factors with roles in plant development and stress response and appear conserved across plant evolution (Li et al. 2013). As in the case of WOX genes, ZFPs have also been implicated in leaf morphology in other species. Superman (SUP), one of the better known ZFPs affecting morphology, suppresses cell division at a specific stage and is associated with flower organogenesis and plant morphogenesis (Jiang et al. 2008). SUP overexpression in Nicotiana tabacum produces a percentage of plants with curled-up leaf margins and thickened roots (Nibau et al. 2011). In plants with moderate expression of SUP and curly leaves, the organisation of the adaxial/ abaxial mesophyll layers was disrupted with more cells and cells of variable shape throughout. Overexpression of a ZFP in rice, OsZHD1, induced abaxial curling due to an increase in the number and abnormal arrangement of bulliform cells (Xu et al. 2014). Other ZFPs, KNUCKLES and JAGGED, have also been shown to be involved in boundary maintenance and/ or cell proliferation. In both Arabidopsis and tomato, the KNUCKLES transcription factor indirectly represses WUSCHEL to allow a temporary inhibition of meristematic activity (Bollier et al. 2018). Ubiquitin E3 ligases are also ZFP family genes. Plant U-box E3 ligases have been implicated in plant development and morphology. As an example, a U-box E3 ligase mutation in barley has been associated with semi-dwarfing, wavy leaf margins and erect growth (Braumann et al. 2018). Cullin4-based E3 ligases are somehow involved in aspects of plant morphology given that mutation of the gene in Arabidopsis causes stunted growth and aberrant leaf shapes (Chen and Hellmann 2013).

Ras-related Protein Rab-8B
The Ras-related Rab-8B gene (Aco015473.1) is positioned only 4.0 kb from the highest associated SNP marker in panel one. Again, data in other species suggests a potential role in leaf margin phenotype in pineapple. The Rab GTPase family of proteins plays an important role in cell wall synthesis and modification through trafficking of cell wall polymers to specific target membranes. Cell wall synthesis and modification such as might occur during cell expansion, are integral to plant development and possibly architecture (Speth et al. 2009;Lycell 2008). A study in Arabidopsis looked at the membrane where RabA localizes, and in conjunction with down-regulation manipulation, demonstrated a corresponding change in cell shape (Kirchhelle et al. 2016). Following on from this it was postulated by Rahni and Birnbaum (2016) that Rab proteins might mediate local wall stiffness thus affecting cell shape and consequently tissue and organ morphology. This potential role in specific cell edge mechanics, rather than that of the entire cell wall, might have implications in pineapple piping leaf margin where a weakness in one edge might allow bending or folding, given growth and turgor pressure in surrounding tissues, especially if cell proliferation is increased in some surrounding tissues.
Although not studied here, it might also explain the position and direction of spines.

Diversity Panel Two
Panel two results represent a difference between spiny-tip and piping leaf margin phenotypes. It is known that the piping margin is dominant over spiny-tip (Collins 1960). The results here (table two, table S3) suggest the piping margin plants might carry a polymorphism in either a WOX1, WOX2 and/ or ZFP that the spiny-tip plants do not. There is no difference between the genotypes represented in panel two with respect to the spiny-tip polymorphisms putatively associated in a previous study (Sanewski 2020). This suggests that both piping and spiny-tip genotypes in panel two carry the spiny-tip polymorphism. Because piping margin is dominant over spiny-tip, piping margin genotypes as derived in this diversity panel must have arisen from spinytip genotypes, not spiny genotypes.

Diversity Panel Three
Because panel three contains all three phenotypes, piping, spiny-tip and spiny, it is expected that the SNPs associated with both piping margins and spiny-tip will be identified as a difference between them and/or spiny margins. Accordingly, two main QTLs were identified (Table 3), both identical to those identified previously in panels one and two. Of some interest is the fact that three different WOX genes were found close to markers associated with piping margins. These were all positioned on chromosome 23 and some distance apart. WOX1 is positioned at 1,057,435-1,059,380 bp, WOX2 is at 738,739-739,215 bp and WOX5 is at 159,566-160,802 bp. There were no linkages observed between markers near these genes and, as such, it is assumed they are complimentary or independently functional. If we assume these different WOX genes are each independently associated with the piping margin phenotype, then they are likely independent selection events. The piping margin phenotype has therefore likely been selected multiple times rather than once. It is possible however, given the panels used here, that some of this selection has occurred in the post-Columbian period.
In the case of the largest and most diverse set of genotypes used here, panel three, SNP marker #54316030 was the most highly associated. This marker was close to a WOX1 and ZFP CONSTANS. When referenced against the spiny-tip 'Smooth Cayenne', piping genotypes were homozygous or heterozygous for a polymorphism. Spinytip and spiny genotypes were homozygous for the 'Smooth Cayenne' variant. Leitão (2018), using a small F1 population, estimated the piping locus to be between 1,177,327 and 1,976,804 bp on chromosome 23. That study used an unrelated population and different mapping strategy but was in very good agreement thus supporting the current study. The QTL identified in that study was close to the second QTL found in the current study which was positioned close to a WUSCHEL related homeobox 1 (Aco007021.1) and Zinc finger CON-STANS-like protein (Aco007020.1). Nashima et al. (2022) used a fine mapping strategy to identify approximately the same QTL on chromosome 23 between positions 2,100,182 and 2,262,401 and nominated one of the same candidate WOX genes (Aco007021.1) with supporting RNAi data but described the gene as a WOX3 not a WOX1 as annotated in the reference genome (Ming et al. 2015).
In a previous study, Sanewski (2020) identified a QTL on chromosome six associated with the spiny-tip leaf margin.
The most likely causal gene was a ZEP, a precursor to ABA, although other genes including a MATE efflux were also implicated. A much larger diversity panel used in the current study comprised of many segregating populations and unrelated genotypes, and using additional analytical methods, but identified the same broad QTL for the spiny-tip phenotype thus validating the earlier study to some extent.
No major SNP markers associated with a piping margin were positioned on chromosome six near the spiny-tip loci in the current study. The piping margin phenotype is however known to suppress other leaf margin phenotypes. This suggests the piping leaf genes function upstream of the causal gene on chromosome six, as far as determination of leaf margin is concerned, or cause a morphological change that prevents spine development as well as induce marginal rolling in the case of piping specifically. What is evident is that piping does not usually occur with spines or spiny-tip. The putative WOX and ZFP polymorphisms are therefore upstream or dominant over the putative ABA pathway polymorphisms regarding leaf morphology.
There were no significantly associated markers identified in panel two genotypes (spiny-tip and piping) on chromosome six where the spiny-tip QTL has previously been shown to be positioned. This suggests there is no difference in the spiny-tip allele between spiny-tip and piping margin plants. This then suggests most of the piping margin plants might carry the same polymorphism for spiny-tip and therefore are descendants or closely related to a spiny-tip cultivar such as 'Smooth Cayenne'. As previously stated, it is believed the piping margin genetics incorporated into the genotypes represented in this study originated from 'Monte Lirio'. This suggests 'Monte Lirio' is related to a spiny-tip cultivar like 'Smooth Cayenne'.

Conclusion
All MLM and logistic regression models were in good agreement regarding the broad QTL and in most cases specific markers most significantly associated with piping leaf margin. Positions within 0.6-1.1 mb on chromosome 23 contained the most significantly associated loci using three logistic regression models, four mixed linear models and three diversity panels. In all three diversity panels, a WOX and ZFP were positioned close to the most significantly associated markers. While other significantly associated markers on other chromosomes were demonstrated, it is likely they are of less importance in most genotypes studied. It appears there are two discreet QTLs associated depending on the diversity panel. One QTL encompasses a WOX2 and ZFP2, the other a WOX1 and ZFP CONSTANS. It is highly probable a WOX and/ or a ZFP are the main gene(s) controlling the piping leaf margin in modern pineapple varieties that originate from the Hawaiian PRI pineapple breeding program. The only piping leaf margin variety used in that program was 'Smooth Guatemala' (syn. 'Monte Lirio') (PRI breeding records) and it is evident that this is the source of the piping phenotype in the current study.
The relatively high level of heterozygosity in pineapple, coupled with strong selection pressure for this trait by plant breeders and early domesticators has likely led to the accumulation of genes which play a role. This raises the possibility that other piping or smooth leaf genotypes not studied in detail here might possess a slightly different genetic determinism as geographical separation and sexual recombination has likely played a role in their development. Most polymorphisms exhibiting a highly significant association for the piping phenotype in this study are all positioned on chromosome 23.
Homologues of several of the proposed candidate genes are known to have functions in important plant processes other than plant morphogenesis in other species. Further study is needed to better understand the role of 'piping leaf genes' in not only that phenotype but also other plant functions important to horticulture.