Gomphrena meyeniana is an extremely variable species from the Andean highlands, which has attracted the attention of many botanists because it is the world’s highest-elevation C4 eudicot and because of its wide morphological variability. It has the typical high-Andean plant morphology, with small leaves tightly clustered on a thick rootstock. The large range of morphological variation within this species, coupled with the varying opinions on the existence of several species or infra-specific taxa, and the lack of molecular information has made the clarification of the G. meyeniana complex a challenge. Our approach was to perform a broad spectrum molecular sampling to identify its phylogenetic position within Gomphrena genus and to perform a multivariate analysis to objectively differentiate taxa based on morphological characters. The ITS and trnL-F regions were analyzed individually and combined using Bayesian inference and maximum parsimony methods. To analyze the morphological characters we performed a clustering method (partitioning around medoids with the Gower’s dissimilarity algorithm). The molecular analyses supported the monophyly of the G. meyeniana complex but did not support segregation into varieties. The morphological analyses gathered the infra-specific taxa into only three varieties that can be easily distinguished through three simple characters: the presence of leaves on the flowering shoot, the habit of the flowering shoot, and the pilosity of the tepals. The varieties of G. meyeniana accepted here are var. meyeniana, var. caulescens, and var. flaccida. A dichotomous key to identify the infra-specific taxa is here presented and illustrated. The varieties tucumanensis, and conwayi were synonymized with var. caulescens, and var. meyeniana respectively.

Molecular Data: Total genomic DNA was extracted from silica-dried leaves using a modified CTAB protocol (Doyle and Doyle 1987). For herbarium specimens, DNA was isolated using the DNeasy plant mini kit (Qiagen, Hilden, Germany) following the manufacturer’s recommendations. ITS1 and ITS2 were amplified in one segment using the primers ITS 28S Fwd and ITS 18S Rev of Kadereit et al. (2005). The chloroplast trnL-F region was amplified in one or two fragments using primers c, d, e, and f described by Taberlet et al. (1991). The PCR amplification for nuclear and cpDNA regions was performed as described by Kadereit et al. (2005) and Bena et al. (2017), respectively. The PCR products were purified with Gel Band Purification Kit (PB-L, Buenos Aires, Argentina), following the manufacturer’s recommendations. Sequencing reactions were performed by Macrogen, Inc. (Seoul, Korea). Low-quality regions were trimmed from the chromatogram files, and reverse and forward sequences for each specimen voucher were assembled and edited using the software MEGA v.6 (Tamura et al. 2013). Sequences for each marker were then aligned using ClustalX v2 (Larkin et al. 2007) with the default options. The alignment was checked and improved manually when necessary using the program MEGA v.6 (Tamura et al. 2013). The aligned data matrix is available from the Dryad Digital Repository (Bena et al. 2019).

Morphological Data: The association between morphological variables was studied by calculating the Pearson correlation coefficient between attributes in the Hmisc package of R (R Core Team 2014). We discarded the highly correlated attributes (correlation value > 0.6) Fig. S3.