Foundational hypotheses addressing plant-insect codiversification and plant defense theory typically assume a macroevolutionary pattern whereby closely related plants have similar chemical profiles. However, numerous studies have documented variation in the degree of phytochemical trait lability, raising the possibility that phytochemical evolution is more nuanced than initially assumed. We utilize proton nuclear magnetic resonance (¹H NMR) data, chemical classification, and double digest restriction-site associated DNA sequencing (ddRADseq) to resolve evolutionary relationships and characterize the evolution of secondary chemistry in the Neotropical plant clade Radula (Piper; Piperaceae). Sequencing data substantially improved phylogenetic resolution relative to past studies, and spectroscopic characterization revealed the presence of 35 metabolite classes. Metabolite classes displayed phylogenetic signal, whereas the crude ¹H NMR spectra featured little evidence of phylogenetic signal in multivariate tests of chemical resonances. Evolutionary correlations were detected in two pairs of compound classes (flavonoids with chalcones; p-alkenyl phenols with kavalactones), where the gain or loss of a class was dependent on the other’s state. Overall, the evolution of secondary chemistry in Radula is characterized by strong phylogenetic signal of traditional compound classes and weak phylogenetic signal of specialized chemical motifs, consistent with both classic evolutionary hypotheses and recent examinations of phytochemical evolution in young lineages.

Study system and sample collection

            For phylogenetic and chemical analyses, we collected leaf material from 71 individuals representing 65 Neotropical Piper species from the following clades: Churumayu (N = 3), Hemipodium (N = 1), Isophyllon (N = 5), Macrostachys (N = 4), Peltobryon (N = 2), Pothomorphe (N = 1), Radula (N = 44), and Schilleria (N = 5). This study complied with all local and national regulations/guidelines, and vouchers for all collections were deposited in herbaria in the country of origin as stipulated in the permit documents (Table S1). Brazilian collections were made under permit No. 15780-6 from the Sistema de Autorização e Informação em Biodiversidade (SISBIO). Costa Rican collections were made under the permits R-054-2018-OT-CONAGEBIO and R-055-2018-OT-CONAGEBIO from the Ministerio del Ambiente y Energía (MINAE). Collections from Ecuador were conducted under the permit 03-IC-FAU/FLO-DNP/MA granted by the Ministerio del Ambiente. Collections from Panamá were covered by the permit SE/AP-15-13 from the Autoridad Nacional del Ambiente (ANAM). Finally, Peruvian collections were covered by the permit 288-2015-SERFOR-DGGSPFFS granted by the Servicio Nacional Forestal de Fauna Silvestre (SERFOR). All collections were identified by E.J.T. in the field, and confirmed with vouchers in the herbarium using regional keys, where available, comparison with type specimens, and experience with the genus. For chemical profiling and DNA sequencing, we collected the youngest, fully expanded leaves and dried them immediately with silica gel. While drying on silica gel may not inhibit enzymatic activity and could limit our analyses to relatively stable molecules, this is not an issue for the phylogenetic analyses described below. Collections were only made from mature individuals in the field. Vouchers were pressed, dried, and deposited in one or more herbaria for future reference and species verification (Table S1). To investigate the evolution of phytochemistry at a relatively shallow evolutionary scale, we conducted the majority of our sampling within Radula³⁴.

Phylogenetic analyses

            Genome-wide polymorphism data was generated for 71 individuals for phylogenetic analyses. Either the same accession sampled for chemical analysis, or an individual from the same population as the one sampled, were sequenced with a genotyping-by-sequencing approach⁶³ that is analogous to ddRADseq⁶⁴. Briefly, genomic DNA was digested with two restriction enzymes, EcoRI and MseI. Sample-specific barcoded oligos containing Illumina adaptors were annealed to the EcoRI cut sites, and oligos containing the alternative Illumina adaptor were annealed to the MseI cut sites. Fragments were PCR amplified and pooled for sequencing. The library was size-selected for fragments between 350 - 450 base pairs (bp) with the Pippin Prep System (Sage Sciences, Beverly, MA), and sequenced on two lanes of an Illumina HiSeq 2500 at the University of Texas Genome Sequencing and Analysis Facility (Austin, TX). Single-end, 100 bp, raw sequence data were filtered for contaminants (E. coli, PhiX, Illumina adaptors or primers) and low quality reads using bowtie2_db⁶⁵ and a pipeline of bash and perl scripts (https://github.com/ncgr/tapioca). We used custom perl scripts to demultiplex our reads by individual and trim barcodes and restriction site-associated bases.

            Assembly and initial filtering was conducted with ipyRAD v.0.7.30⁶⁶. ipyRAD was specifically designed to assemble ddRADseq data for phylogenetic applications, permits customization of clustering and filtering, and allows for indel variation among samples⁶⁶. Because a suitable Piper genome was not available at the time of analysis, we generated a de novo consensus reference of sampled genomic regions with ipyRAD. Briefly, nucleotide sites with phred quality scores lower than 33 were treated as missing data. Sequences were clustered within individuals according to an 85% similarity threshold with vsearch⁶⁷ and aligned with muscle⁶⁸ to produce stacks of highly similar ddRADseq reads (hereafter, ddRADseq loci). The sequencing error rate and heterozygosity were jointly estimated for all ddRADseq loci with a depth >6, and these parameters informed statistical base calls according to a binomial model. Consensus sequences for each individual in the assembly were clustered once more, this time across individuals, and discarded if possessing >8 indels (max_Indels_locus), >50% heterozygous sites (max_shared_Hs_locus), or >20% variable sites (max_SNPs_locus). To reduce the amount of missing data in our alignment matrix, ddRADseq loci were retained if they were present in at least 50 of 71 samples. The nexus file of concatenated consensus sequences for each individual, including invariant sites, were used as input for the Bayesian phylogenetic methods described below. The nexus alignment as well as complete information on additional parameter settings for this analysis are archived at Dryad (https://doi.org/10.5061/dryad.j6q573nc7).

            To resolve patterns of diversification and to provide a foundation for investigating variation in patterns of phytochemical evolution, we estimated a rooted, calibrated tree according to a relaxed clock model in RevBayes v.1.0.12⁶⁹, which provides the ability to specify custom phylogenetic models for improved flexibility compared with other Bayesian approaches. The prior distribution on node ages was defined by a birth-death process in which the hyper priors on speciation and extinction rates were exponentially distributed with λ = 10. We relaxed the assumption of a global molecular clock by allowing each branch-rate variable to be drawn from a lognormal distribution. After comparing the relative fits of JC, HKY, GTR, and GTR+Gamma nucleotide substitution models with Bayes factors, we modeled DNA sequence evolution according to the best-fit HKY model. Eight independent MCMC chains were run for 100,000 generations with a burn-in of 1,000 generations and sampled every 10 generations. Chains were visually assessed for convergence with Tracer v.1.7.1⁷⁰ and numerically assessed with effective sample sizes (ESS), the Gelman−Rubin convergence diagnostic⁷¹, and by comparing the posterior probabilities of clades sampled between MCMC chains. The maximum clade credibility (MCC) tree provided the ultrametric fixed tree topology and relative node ages for phylogenetic comparative methods described below.