DOI: 10.5061/dryad.j6q573nc7
title: Phytochemistry reflects different evolutionary history in traditional classes versus specialized structural motifs

description: 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 (1H 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 1H 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.
description: <b>Study system and
 sample collection</b>             For
 phylogenetic and chemical analyses, we collected leaf material from 71
 individuals representing 65 Neotropical <i>Piper </i>species
 from the following clades: Churumayu (<i>N</i> = 3),
 Hemipodium (<i>N</i> = 1), Isophyllon (<i>N</i> =
 5), Macrostachys (<i>N</i> = 4), Peltobryon
 (<i>N</i> = 2), Pothomorphe (<i>N</i> = 1), Radula
 (<i>N</i> = 44), and Schilleria (<i>N</i> = 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<sup>34</sup>.   <b>Phylogenetic
 analyses</b>             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<sup>63</sup> that is
 analogous to ddRADseq<sup>64</sup>. Briefly, genomic DNA was
 digested with two restriction enzymes, <i>Eco</i>RI and
 <i>Mse</i>I. Sample-specific barcoded oligos containing
 Illumina adaptors were annealed to the <i>Eco</i>RI cut sites,
 and oligos containing the alternative Illumina adaptor were annealed to
 the <i>Mse</i>I 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
 (<i>E. coli</i>, <i>Phi</i>X, Illumina adaptors or
 primers) and low quality reads using bowtie2_db<sup>65</sup>
 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<sup>66</sup>. ipyRAD was specifically designed
 to assemble ddRADseq data for phylogenetic applications, permits
 customization of clustering and filtering, and allows for indel variation
 among samples<sup>66</sup>. Because a suitable <i>Piper
 </i>genome was not available at the time of analysis, we generated a
 <i>de novo </i>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<sup>67</sup> and aligned with
 muscle<sup>68</sup> 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
 &gt;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 &gt;8 indels (max_Indels_locus), &gt;50%
 heterozygous sites (max_shared_Hs_locus), or &gt;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<sup>69</sup>, 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
 <i>λ </i>= 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<sup>70</sup> and
 numerically assessed with effective sample sizes (ESS), the Gelman−Rubin
 convergence diagnostic<sup>71</sup>, 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.