In plants, it has been reported that high temperature has negative effects in various physiological processes such as photosynthesis, primary and secondary metabolism, water relations, or lipid metabolism (Xu et al., 2006). Specifically, it is known that heat stress generates damage to the cell membrane, overproduction of reactive oxygen species (ROS), senescence, inhibition of photosynthesis, and cell death. However, plants have designed specific protection mechanisms to minimize and repair the damage caused by high temperatures in order to cope with heat waves that they will inevitably face during their lifespan.

When a plant perceives environmental stress, multiple signaling cascades are activated. These include the intricate crosstalk between the different plant hormones and other signaling pathways involving kinases and phosphatases, calcium, ROS, and lipids. Previous studies in P. radiata showed salicylic acid (SA) and abscisic acid (ABA) as crucial factors in the initial response to heat stress (Escandón et al., 2016), probably due to the rush of the plant to quickly regulate stomatal closure as observed in other species (Acharya and Assmann, 2009). Other phytohormones such as cytokinins (CKs) and indolacetic acid (IAA) seem to be more important for the acclimation and recovery of the plant (Skalák et al., 2016). Nevertheless, how these hormones interact and trigger changes in the plant metabolome to deal with heat stress is still unknown.

Environmental stress results in a reorganization of the metabolism in order to assure homeostasis, which is often accomplished by maintaining essential metabolism and synthesizing metabolites with stress-protective and signaling properties (Fernández de Simón et al., 2017). Changes in secondary metabolism are usually triggered by environmental stresses (Routaboul et al., 2012; Fernández de Simón et al., 2017) such as extreme temperatures, salinity, water availability, or high light intensity, being involved in most plant adaptive responses in combination with signaling mechanisms regulated by plant hormones. Phenolic compounds are one of the most important classes of secondary metabolites in plants as they play important roles in the response to high-temperature stress. Specifically, high-temperature stress promotes the production of phenolic compounds such as flavonoids, phenylpropanoids, anthocyanins, and lignins which are related to the suppression of stress-induced oxidation of most cell molecules (Wahid et al., 2007).

Lipid metabolism is also altered by high-temperature stress as a basic mechanism to control membrane fluidity, cell signaling, and movement of substances (Falcone et al., 2004). According to Grover et al. (2000) and Umesha (2005), the accumulation of highly saturated fatty acids might confer tolerance to hightemperature stress by means of reducing structural membrane fluidity when this is increased by environmental warmth.

Metabolomics profile analyses can be considered an important diagnostic tool for verifying the physiological responses of plant species to environmental changes and to understand the mechanisms behind the complex biological response (Gibbons et al., 2015; Zhang et al., 2015; Meijón et al., 2016), providing a snapshot of the physiological status of the plant in response to an environmental stress (Gibbons et al., 2015; Patel and Ahmed, 2015; Zhang et al., 2015). The feasibility of metabolome analysis for biomarker discovery relies on the assumption that metabolites are important players in biological systems (Monteiro et al., 2013) and stress situations cause drastic changes of metabolomics pathways, which are not new concepts (Shulaev et al., 2008; de Leonardis et al., 2015; Pascual et al., 2017). Currently, using MS-based platforms and combining different analytical technologies, it is possible to increase metabolome coverage. The use of GC-MS technique allows measuring the most of the primary metabolites, while LC-MS provide a better coverage of large hydrophobic metabolites predominant in secondary metabolisms (Doerfler et al., 2013). Both techniques together can assure a better metabolome, allowing the elicitation of a complete view of metabolic dynamics involved in the heat response of P. radiata. However, these kinds of studies also require a system biological approach using bioinformatics tools to understand their implications for cell function and to attach the missing connections between molecules and plant physiology (Bruggeman and Westerhoff, 2007; Meijón et al., 2016).

In this study, P. radiata plants were exposed to high temperatures aiming to mimic future scenarios of increased warmth. Metabolome and physiological data were analyzed comprehensively using a multivariable approach combining both sets of data. The analysis revealed the dynamic behavior of the metabolic and signaling transduction pathways, as well as connections among the pathways. Classical physiological measurements combined with cutting-edge technologies such as mass spectrometry-based analytical procedures for characterizing the variations in the metabolome revealed key metabolites related to high-temperature response in P. radiata. These key metabolites have a possible use as biomarkers in P. radiata and other species, due to the similarity of metabolites and basic metabolic pathways between very different species, while proteins, genes, and mRNAs are diversified from one species to another (Peng et al., 2015). Altogether, this work provides a deep knowledge of the response and acclimation process to heat stress in P. radiata, as well as the selection of possible universal thermotolerance biomarkers or crucial metabolites which can be further considered by breeders and forest managers.

# MATERIALS AND METHODS

#### Plant Material and Experimental Design

The assay was conducted in a climate chamber under controlled conditions (Fitoclima 1200, Aralab). One-year-old P. radiata seedlings (plant size about 33 ± 4 cm) in 1-dm3 pots (blond peat:vermiculite, 1:1) were kept under a photoperiod of 16 h

(400 µmol m−2 .s −1 ) at 25◦C and 50% relative humidity (RH), and 15◦C and 60 % RH during the night period. The plants had been previously acclimated over a 1-month period inside the climate chamber, being watered with nutritive solution (NPK, 5:8:10).

Control plants (C) were collected before starting the heat exposure and were maintained at 25◦C for the duration of the trial. Heat exposure treatment began with an increasing temperature gradient from 15 to 40◦C over 5 h and maintained for 6 h. This experimental procedure was repeated for 5 days. Sampling was performed at: 3 h after 40◦C was reached on day 1 (T1/2) and at the end of the 6-h heat exposure on day 1 (T1), day 2 (T2), day 3 (T3), and day 5 (T5). Plants were watered every day to 80% FC (full capacity) and weekly fertilized with a nutritive solution (NPK, 5:8:10). Plants of each treatment were allowed to recover for 1 month under the control conditions. Recovered plants (R) represent an intermediate exposure time (T3 recovered plants) because there was no significant difference (at the morpho-physiological level) between recovered plants exposed to the different heat exposures. This experimental design aimed to cover the entire stress sensing–response–adaption process, increasing the density of analysis at short-term (T1/2, T1, T2), longer exposures (T5), and recovered plants (R), complementing previous short-term response (C, T1, and T3) analysis (Escandón et al., 2017).

Mature needles from each plant (16 plants/exposure) were sampled, cleaned with a moistened cloth, and immediately frozen in liquid nitrogen until metabolites were extracted. Pools of 3 plants for each biological replicate were performed. Data of physiological measurements – electrolyte leakage (EL), relative water content (RWC), maximum quantum efficiency of PSII (Fv/Fm), quantum yield of photosystem II photochemistry (ϕPSII), malondialdehyde content (MDA), proline content, starch content, total soluble sugars (TSS), chlorophyll a content (Chla), chlorophyll b content (Chlb), and carotenoid content (Carot) – and fitohormone data – SA, indol-3-acetic acid (IAA), ABA, zeatin riboside (ZR), dihydrozeatin riboside (DHZR), gibberellin 7 (GA7), jasmonic acid (JA), gibberellin 9 (GA9), isopentenyl adenosine (iPA), isopentenyl adenine (iP), and castasterone (BK) – were taken from Escandón et al. (2016) for multivariable and integrative analysis.

#### Metabolite Extraction

Metabolite extraction was performed according to Valledor et al. (2014a) using 100 mg of needle fresh weight. Briefly, samples (C, T1/2, T1, T2, T3, T5, R) were ground in liquid nitrogen and 600 µL of cold (4◦C) metabolite extraction solution – Methanol:Chloroform:H2O (2.5:1:0.5) – was immediately added to each tube. Then, samples were centrifuged at 20,000 g for 4 min at 4◦C and the supernatant transferred to new tubes. Finally, 800 µL of Chloroform:water (1:1) were added and the tubes vortexed and again centrifuged at 20,000 g for 4 min at 4 ◦C. Two layers formed: an upper aqueous layer, containing the polar metabolites and a lower organic layer, containing the nonpolar. Both fractions were transferred to new tubes and dried in a speed vac.

# Polar Metabolite Identification and Quantitation Using LC-Orbitrap-MS Analysis

The polar fraction of each sample was analyzed twice on an LC-Orbitrap-MS, first in positive ion mode and then in negative. A Dionex Ultimate 3000 (Thermo Fisher Scientific, United States) UHPLC was used and a LC-Orbitrap LTQ XL-MS system (controlled by Xcalibur version 2.2, Thermo Fisher Corporation) was run according to the procedure described in Escandón et al. (2017). The resolution and sensitivity of the Orbitrap were controlled by the injection of a mixed standard after the analysis of each batch, and the resolution was also checked with the aid of lock masses (phthalates). Blanks were also analyzed during the sequence.

LC-Orbitrap-MS raw data were processed and compared using MZmine software version 2.10 (Pluskal et al., 2010). MS1 spectra were filtered establishing a noise threshold at 5.5E03 and minimum peak height at 6E03 with a minimum time peak of 0.15 min. Peaks were smoothed and deconvoluted by using a local minimum search algorithm (98% chromatographic threshold, minimum retention range 5 min, minimum relative height of 90%, and minimum ratio top/edge of 1.2). Chromatograms were aligned using the RANSAC algorithm with a tolerance of 5 ppm of and 1.0 min retention time. Normalized peak areas were used for quantification, and their values were log transformed before statistical analyses (Supplementary Table S1a).

The individual peaks were identified following different approaches; the first step was performed against an in-house library (>100 compounds) and manual annotation considering m/z and retention times. In the second step, masses were assigned using the KEGG, PubChem, METLIN, MassBank, HMDB, and Plantcyc databases as reported by Escandón et al. (2017) with built-in MZmine plugins with a 5 ppm threshold and considering as "identified" beyond doubt those metabolites that were defined after the comparison to our standard compound library or by a matching of MS/ MS to the small number of plant compounds for which their MS/MS is available in public databases (Meijón et al., 2016); and as "tentatively assigned" those with molecular ions with exact masses corresponding to identified metabolites in databases. Metabolite identification against our library was confirmed by retention time (RT), mass, isotopic pattern, and ring double bound parameters. Supplementary Data S1 includes the detailed interpretation of experimental MS/MS spectra which support our tentative identifications of the candidate metabolites that were not identified beyond doubt in the first term.

### Non-polar Metabolite Identification and Quantitation Analysis

Non-polar metabolites were derivatized with 295 µL tertmethyl-Butyl-Ether (MTBE) and 5 µL of trimethylsulfonium hydroxide (TMSH) for 30 min at room temperature. The tubes were centrifuged for 3 min at 20,000 g to remove insoluble particles before transferring the supernatants to GC-microvials. GC-MS measurements were carried out following a previously developed procedure (Valledor et al., 2014b) on a triple quad instrument (TSQ Quantum GC; Thermo, United States). The

mass spectrometer was operated in electron-impact (EI) mode at 70 eV in a scan range of m/z 40–600. Metabolites were identified based on their mass spectral characteristics and GC retention times through comparison with the retention times of reference compounds in an in-house reference library and the current version of Golm Metabolome Database (Hummel et al., 2007) using LC-Quant software (Supplementary Table S1b).

#### Quantitative Real-Time PCR of Selected Genes

RNA was extracted from 100 mg of needle fresh weight as described by Valledor et al. (2014a). cDNA was obtained from 1,000 ng of RNA using the RevertAid kit (Thermo Scientific, United States) and random hexamers as primers following the manufacturer's instructions. Later, qPCR reactions were performed in a CFX Connect Real Time PCR machine (Bio-Rad) with SsoAdvanced Universal SYBR Green Supermix (Bio-Rad, United States); three biological and two analytical replicates were performed for each treatment.

ACTIN (ACT), RIBOSOMAL PROTEIN 18S, GLYCERALDEH YDE 3-212 PHOSPHATE DEHYDROGENASE (GAPDH), and UBIQUITINE (UBI) genes were tested as endogenous control employing geNorm following the criteria of Hellemans et al. (2007). ACT and UBI were the most stable and consequently selected as endogenous genes. Normalized Relative Quantities (NRQ) and Standard Errors of RQ were determined according to Hellemans et al. (2007). Primers were designed using transcript sequences available in a P. radiata in-house database obtained from RNA-Seq data (unpublished results). Detailed information about the primers used for qPCR experiments is available in Supplementary Table S2.

#### Statistical and Bioinformatics Analysis

Five biological replicates were used for metabolites and physiological parameter statistical analysis. The procedures were conducted with the R programming language running under the open-source computer software R v2.15.2 (R Development Core Team, 2015) and RStudio (RStudio Team, 2015). Metabolome datasets were pre-processed following the recommendations of Valledor et al. (2014c). In brief, missed values were imputed using a k-nearest neighbors approach, and variables were filtered out if they were not present in all replicates of one treatment or in at least 45% of the analyzed samples. Data was transformed following a sample-centric approach followed by log transformation. Centered and scaled values (z-scores) were subjected to multivariate analysis and Heatmap clustering. The calculation of the number of common metabolites for all combinations of treatments and unique metabolites in a single treatment was performed using core functions of R.

Metabolomics pathways of each metabolite (Supplementary Table S3a) were searched against KEGG pathway maps (KEGG Mapper, Kanehisa et al., 2012) and p-values of each metabolomics pathways (Supplementary Table S3b) in MBROLE 2.0 (López-Ibánez et al., 2016). Heat mapping was carried out using the Manhattan distance method to group metabolites in different KEGG pathways with an MBROLE FDR correction of less than 0.05. Multivariate analysis of metabolites and physiological parameters (Supplementary Table S4) were conducted with mixOmics (Lê Cao et al., 2009) using Principal Component Analysis (PCA), Sparse Partial Least Squares (sPLS), and network analyses. The normalization of the datasets was performed before combining them. sPLS algorithm was used to find correlations between predictor (metabolites matrix) and response variables (physiological parameters) and its graphic representation in the network analysis. Network topology was defined after applying sPLS regression using the function network() of the mixOmics package and filtered (only edges equal or higher than |0.60| were maintained) in Cytoscape v.3.3.0 (Cline et al., 2007). Univariate analyses were conducted: one-way ANOVA, p < 0.05 for metabolites and physiological parameters and Student's t-test, p < 0.1 for qPCR analysis. Graphics were plotted employing ggplot2 (Wickham, 2009) and pheatmap (Kolde, 2015).

# RESULTS

# Characterization of the Metabolome During Heat Treatment

The combination of LC-Orbitrap-MS for polar and GC-MS for non-polar metabolites allowed the reliable quantification of a total of 2,287 ions based on the obtained m/z and retention times (Supplementary Table S1). These different ions can be considered different metabolites since each combination of m/z and retention time should be unique for each metabolite and spatial conformation. The fusion of a customized searching/identification algorithm based on in-house and public databases resulted in the unequivocal identification of 41 ions (identical matches to our compound library) and 747 ions that were tentatively assigned after comparing its very accurate mass against reference compound databases (Supplementary Table S1). The combination of both ionization modes (GC-MS and LC-Orbitrap-MS) gave a broad characterization of the pine metabolomes during heat-induced response, which covered most of the primary and secondary metabolism pathways.

Metabolome analysis showed that most of the metabolites are present in all treatments or at least in two different treatments (880 or 2,014 metabolites, respectively, "common" in **Figure 1** and Supplementary Table S5). Conversely, a total of 273 metabolites were identified only in one of the treatments (**Figure 1**), Control (C) being the sample that showed the highest number of characteristic metabolites (78), followed by T5 (44) and R (42). Otherwise, T5 and C were the treatments that shared a greater number of metabolites (106 metabolites; Supplementary Table S5).

The complexity of the metabolome data was reduced by focusing on the specific pathways in relation to the relative abundance of the metabolites identified. Heatmap-clustering analysis (**Figure 2** and Supplementary Figure S1), distinguished four different groups in relation to metabolic pathways identified: C, T1/2-T1, T2-T3, and T5-R group that is highly separated of C and shorter-term treatments.

On the other hand, the clustering based on the profile of KEGG pathways did not show any clear trend, each group presenting a different pathway profile across the sample times, e.g., glucosinate biosynthesis and flavone and flavonol biosynthesis increased their activities in treatments T1/2, T1, T2, or T3; conversely, pentose phosphate pathway and flavonoid biosynthesis showed the maximum in T5 and/or R. Other pathways, such as the biosynthesis of unsaturated fatty acids increased the activity in T2, T3, and T5 or only in the longerterm exposure (T3 and T5) as arginine and proline metabolism pathways. In contrast, pathways related with phenylalanine metabolism enhance their activity during the stress but decrease at T5 (at the levels of C and R). Lastly, riboflavin metabolism showed a substantial increase only in T5 and glutathione metabolism showed a fluctuating tendency, highlighting its maximum in T3.

These results seem to underlie two different types of heatresponse in P. radiata (shorter-term and longer-term response) where different pathways are necessary at a different time of stress.

#### Shorter-Term and Longer-Term Responses to Heat Stress Were Confirmed by Integrative Analysis of Metabolome and Physiological Datasets

To simplify the dimensionality of the results and integrate metabolome data with physiological datasets, sPLS (**Figure 3A** and Supplementary Table S6) and PCA (Supplementary Figure S2 and Supplementary Table S7) analyses were used. sPLS and PCA scores revealed a trajectory of the different sampling times by the combination of the two main components. Recovered (R) and long-term exposed (T5) plants were separated from the other treatments by the first component, which seems to be gathering the variance related to long-term stress adaption, while the second component, considering top correlated variables (Supplementary Figure S2 and Supplementary Table S6) is related to heat-response/tolerance. This component revealed differences between Control-Recovered and heat-treatments.

Altogether, multivariate results provided hints about two different responses to stress. First, an initial response to stress (shorter-term response), involves increased ABA and SA activities in samples T1/2, T1, T2, and T3 (**Figure 3A**, Supplementary Table S6a, Comp 1). Secondly, adaptive mechanisms, such as TSS and ZR, seem to be involved in the distinction between T5 and T3, as well as between the recovered and control plants according to sPLS analysis (**Figure 3A**, Supplementary Table S6a, Comp 1).

The interaction networks built from this analysis (**Figure 3B** and Supplementary Table S6) showed a complex correlation between different metabolites and hormones and physiological parameters. The important role of CKs was demonstrated during heat response, as well as the links of these hormones with different compounds of the metabolome and physiological parameters in response to high temperatures. Two main nodes in the constructed network were detected: ZR and iPA nodes. The ZR node was negatively correlated with numerous saturated and unsaturated fatty acids (C16:0 and C18:0 with their unsaturated forms), flavonoids (including kandelin A-1 or eujambolin), amino acids (such as L-proline and L-serine), and L-phenylalanine (key metabolite in several pathways, such as phenylalanine metabolism, phenylpropanoid biosynthesis and biosynthesis of plant hormones). The presence of L-proline in the network indicates the relevance of the accumulation of this amino acid in longer-term response (T3 and T5) to high temperatures in P. radiata, as it has already been showed in Escandón et al. (2016) and confirmed by the increased activity of its metabolism (**Figure 2**, arginine and proline metabolism). Additionally, ZR was positively correlated with fatty acids involved in sphingolipid metabolism, terpenoids (including abscisic-alcohol 11-glucoside) and other secondary metabolites: dihydrokaempferol (related to flavonoid biosynthesis) and cis-1,2-dihydro-3-ethylcatechol (implicated in the degradation of aromatic compounds). In this node, it is also important to note that the physiological parameter, EL, appeared negatively correlated by an unknown metabolite (p537) to IAA and DHZR which, in turn, were also positively linked to tetracosanoic acid (saturated fatty acid, C24:0). This confirms the relevance of these hormones (Du et al., 2013;

to KEGG pathway. C, control; R, recovered; T1/2, 3 h after 40◦C on day 1; T1, 6-h heat exposure on day 1; T2, day 2; T3, day 3; T5, day 5.

Cerný et al., 2014 ˇ ) and fatty acid signaling (Los and Murata, 2004) to repair membrane damage related to heat stress. On the other hand, iPA was positively correlated with fatty acids, such as tetradecanoic acid and (9Z)-octadecenoic acid, as well as metabolites involved in diterpenoid biosynthesis (sclareol); and this hormone was negatively correlated with compounds involved in flavonoid biosynthesis, such as dihydromyricetin.

Interestingly, the analysis of the network dynamics and particularly the quantitation of the represented variables showed a two-step response, with T3 as the transition point between shorter-term and longer-term responses (Supplementary Movie S1). This observation is consistent with the conducted sPLS and PCA analyses. This behavior can be considered an adaptive mechanism to cope with rapid environmental changes that occur daily (i.e., sun heat following a rainy period). In this case, plants require a mechanism to quickly overcome the first impact of stress; however, the physiology should return to ideal values after the removal of the stress factor

control that this node exerts over the interactions of other nodes in the network (higher control = lighter color). EL, electrolyte leakage; ZR, zeatin riboside; iPA, isopentenyl adenosine; IAA, indol-3-acetic acid; DHZR, dihydrozeatin riboside; BK, castasterone; GA7, gibberellin 7; C, control; R, recovered; T1/2, 3 h after 40◦C on day 1; T1 6-h heat exposure on day 1; T2, day 2; T3, day 3; T5, day 5.

in order to achieve an optimal energetic balance. A previous work (Escandón et al., 2016) showed that photosynthetic state was only slightly affected in the first impact to heat stress (T1/2, T1 and T2), recovered to control levels in T3, and even improved on T5, probably related to the beginning of the acclimation process; while lipids peroxidation analysis showed a slight accumulation of MDA in shorter-term exposures (T1/2, T1 and T2) prior to the activation of acclimation mechanisms. On the other hand, if the stress persists the plant must adapt to the new situation with the cost of reducing its growth and reproductive capacity compared to an ideal situation (Bradford and Hsiao, 1982).

#### Using Metabolomics to Explore Possible Thermotolerance Biomarkers

The metabolome analysis revealed the essential role of specific metabolites (**Figure 3**) and pathways (**Figure 2** and Supplementary Figure S1) in relation to high-temperature adaptation in P. radiata.

The accumulation of three different metabolites belonging to pathways with higher significant changes across hightemperature treatments is showed in **Figure 4**. One of these three key compounds is L-phenylalanine which participates in numerous pathways (Meyermans et al., 2000; Wittstock and Halkier, 2000; Teufel et al., 2010; Tzin and Galili, 2010; Yoo et al., 2013), such as phenylalanine metabolism, phenylpropanoid biosynthesis, phenylalanine, tyrosine and tryptophan biosynthesis, glucosinolate biosynthesis, and biosynthesis of plant hormones, most of them highly related/interconnected. Its accumulation profile reflected a level increase from T1/2 to T3 (**Figure 4A**), decreasing to control levels in T5 and reaching the lowest accumulation values in R. Phenylpropanoids and flavonoids play a key role in protecting plants against abiotic stress, largely by inhibiting the formation of ROS through a number of different mechanisms (Mierziak et al., 2014). The increment of L-phenylalanine in the shorter treatments may be related to the need to cope with the production of ROS. The decrease in L-phenylalanine in T5 and R could be related to the increased activity of PHENYLALANINE AMMONIA LYASE (PAL; crucial enzyme in phenylpropanoid pathway), which transforms L-phenylalanine into trans-cinnamic acid (Rivero et al., 2001). The activity increase of this enzyme is considered one of the most important ways of cell acclimation against stress in plants (Levine et al., 1994; Leyva et al., 1995; Rivero et al., 2001).

Other pathway highlighted by heatmap-clustering (**Figure 2** and Supplementary Figure S1) was the glutathione metabolism, which showed the highest accumulation in T3. Ascorbic acid is key in this pathway; however, it showed decreased levels in short-term treatments (**Figure 4B**), reaching the highest values in T5 and R. Ascorbic acid and glutathione are both antioxidants, which are crucial for plant defense against oxidative stress (Noctor and Foyer, 1998). The decrease of ascorbic acid in shorter treatments could be explained by their use in this defense against oxidative damage.

D-(−)-ribose (**Figure 4C**), which participates in pentose phosphate and ABC transporters pathways, showed a decline in shorter-term, recovering the control values in T5 and R. This may indicate that plants reduce their metabolism until T3, recovering the activity in T5 and R when the plants are adapted. day 3; T5, day 5.

fpls-09-00485 April 17, 2018 Time: 13:1 # 8

Although sugars play an important role against heat stress in many species (Wahid et al., 2007), in P. radiata it has been showed that the total amount of soluble sugars tends to decrease in the first moments of stress (Escandón et al., 2016). This may be because plants maintains growth patterns although they are driven by consuming carbohydrates reserves (Mitchell et al., 2013; Escandón et al., 2016).

Under high temperatures, plants alter lipid composition, causing membranes to become more fluid and thus interrupting membrane processes (Falcone et al., 2004). High-temperaturetolerant plants show an increased presence of saturated fatty acids to counteract increased fluidity during heat stress (Grover et al., 2000; Umesha, 2005). Levels of candidate metabolites related to biosynthesis of unsaturated fatty acids pathway are represented in **Figure 5** including both saturated fatty acids (**Figures 5A,D**) and unsaturated fatty acids (**Figures 5B,C,E,F**). Saturated fatty acids, hexadecanoic acid (C16:0, **Figure 5A**) and C18:0 (**Figure 5D**), shared an increasing tendency in short-term response, starting to fall in T3 (in the case of C16:0) or earlier in T2 (C18:0). In both, the lowest accumulation values were reached in R. This saturated fatty acids raise is consistent with the hypothesis by Grover et al. (2000), which indicates the plants increase saturation until seems to be already acclimated to the heat stress.

In the case of unsaturated fatty acids, they showed different patterns depending on the studied fatty acid. C16:2 and C18:2 showed a slight decrease in short-term response (**Figures 5B,E**), conversely, C16:3 and C18:3 (**Figures 5C,F**) presented increased levels in T1/2 and T1. In longer-term response, the most common tendency is the recovery at the control values (**Figures 5B,C,E**) even surpassing them, except for C18:3 where this only occurs in R (**Figure 5F**). This differential accumulation patterns seems to evidence the different roles of each unsaturated fatty acid during heat-stress response.

Flavonoids related pathways are also crucial elements according to heatmap-clustering analysis (**Figure 2** and Supplementary Figure S1). Bibliography indicate that flavonoids and anthocyanins are essential compounds to prevent and protect the plants against different abiotic and biotic stress (Dao et al., 2011; Zhang et al., 2012). Results showed that dihydromyricetin (**Figure 6A**), key metabolite in flavonoid biosynthesis pathway, continuously decreases its levels until T3, to then begin its accumulation, reaching the control values in R. On the contrary, flavonoids eujambolin and kandelin A-1 (**Figures 6B,C**), showed the lowest accumulation in T5 and R. Eujambolin presented a decreasing tendency in shorter treatments, showing a strong decrease in T5 and later in R. Conversely, kandelin A-1 showed an increase in shorter treatments in relation to control, displaying a strong decrease in T5 and R.

Dihydromyricetin (along with dihydroquercetin and dihydrokaempferol) is a dihydroflavonols implicated in the synthesis of anthocyanidins (Falcone Ferreyra et al., 2012). While, kandelin A-1 is a proanthocyanidin which are synthesized as oligomeric or polymeric end products of one of several branches of the flavonoid pathway, which shares the same upstream pathway with anthocyanins (He et al., 2008). Anthocyanins are usually accumulated during heat stress in vegetative tissues (Wahid and Ghazanfar, 2006) in order to decrease the transpirational losses caused by lower osmotic potential of the leaf (Chalker-Scott, 2002). However, the complexity of its metabolism and the high number of these compounds makes their study difficult.

Gene expression of PAL (**Figure 7A**), two DESATURASE (DES, **Figures 7B,C**) and a central gene involved in anthocyanins and proanthocyanins biosynthesis pathways, DIHYDROFLAVONOL 4-REDUCTASE (DFR, **Figure 7D**) (Ayabe et al., 2010; Katsu et al., 2017) were analyzed in order to confirm the possible role of these elements as biomarkers and validate the relevance of L-phenylalanine, fatty acid and flavonoid metabolism in the high-temperature response of P. radiata.

These data showed that PAL (**Figure 7A**) increases its expression in the first moments to stress (T1/2 and T1), returning to the control values in T2 and T3. In T5 and R

FIGURE 6 | Levels of key metabolites related to flavonoid biosynthesis across high-temperature treatments: Dihydromyricetin (A), eujambolin (B), and kandelin A-1 (C). Box plot representation of the LOG transform data (Supplementary Table S4). ID of the three metabolite, previously tentatively identified using public database, was validated by the interpretation of their MS/MS spectra (Supplementary Data S1). C, control; R, recovered; T1/2, 3 h after 40◦C on day 1; T1, 6-h heat exposure on day 1; T2, day 2; T3, day 3; T5, day 5.

its expression increases drastically, confirming the pattern of decrease of L-phenylalanine (**Figure 4A**) in these treatments by its consumption. In the case of DES (**Figures 7B,C**), both genes showed a reduction of the expression in shorter-term treatments: contig5799 presented the decrease in T1/2, T1 and T2, while contig04128 only in T2. This reduction is consistent with the

increase of saturated fatty acids identified (C16:0 and C18:0, **Figures 5A,D**). Finally, DFR (**Figure 7D**) showed an increase of expression across of stress time, reaching a peak in T5 and return to the control values in R. These results are according dihydromyricetin levels quantified (**Figure 6A**), which seems to be being consumed for the production of anthocyanins in response to high-temperature stress.

#### DISCUSSION

#### Metabolome Characterization: Dynamic of Metabolism Throughout High-Temperature Stress

The plant metabolome responds to an unfavorable environment in a dynamic way, being favored a characteristic type of metabolic pathways at every moment of stress. Depending on the timing, different compounds can be identified, as stress signal transduction molecules, stress metabolism by-products, or molecules that are part of the plant acclimation response (Shulaev et al., 2008). According to this, the results of this work showed that the response of P. radiata at high temperatures activates the synthesis of specific metabolites at each time of the stress (**Figure 1**) showing every sample-time a significant number of unique metabolites. Additionally, these results seem to reveal the beginning of the acclimation process in T5 that showed the greatest number of unique metabolites but also the greatest number of shared metabolites with the control plants (106 metabolites, Supplementary Table S4).

In-depth analysis of metabolite accumulation in relation to KEGG pathways (**Figure 2** and Supplementary Figure S1) allowed to confirm the high dynamism of the metabolome in relation to high-temperature response (Escandón et al., 2017). Although two main trends in the KEGG pathways have been identified by heatmap-clustering (**Figure 2** and Supplementary Figure S1), each pathway seem to have a specific role in a particular moment of the stress, seemingly following an orchestrated succession, particularly, in the case of flavonoids related pathways. Thus, phenylpropanoid biosynthesis seems to be activated in the first contact with stress (T1/2-T1), while in T2-T3 the most active route is flavone and flavonol biosynthesis, to finally increase the activity of the general pathway of flavonoid biosynthesis in T5 and R. The importance of the flavonoids across the heat stress is evident. According to bibliography the compounds in this group act as antioxidants, contributing to the adaptation to environmental changes such as cold, high temperatures or irradiation (Gould et al., 2002; Doerfler et al., 2013, 2014; Meijón et al., 2016).

Fatty acids also play a fundamental role in the response to stress in P. radiata. Biosynthesis of unsaturated fatty acids is a very active pathway in the longer exposure treatments (T2, T3 and T5), while saturated fatty acid biosynthesis is more relevant in the first moments of exposure to heat (T1/2 and T1). A reduction of unsaturated fatty acid content and increase in saturated fatty acids content has been positively associated with heat tolerance as they counteract the increase in membrane fluidity caused by high temperature (Grover et al., 2000; Larkindale and Huang, 2004). However, this hypothesis has not been fully confirmed since that unsaturated fatty acid are also known to be key elements in heat stress signaling (Königshofer et al., 2008) and to provide other essential characteristics to lipid membrane (Falcone et al., 2004), showing an increase in its accumulation at high temperatures.

In the case of P. radiata, biosynthesis of fatty acid is a key mechanism to overcome the fluidization of the membrane until T2 at the beginning of the stress, when the plants had overcome this fluidization. The global increment in saturated fatty acid would be consequence of the increase of the amount of a specific saturated fatty acid (e.g., hexadecanoic acid) or the accumulation of fatty acids with lower numbers of double bonds (the reduction of C18:3 to C18:2 or C18:1; Wang et al., 2017). From T2, plants activate the biosynthesis of unsaturated fatty acids which may be necessary to stabilize photosynthesis (Gombos et al., 1994) or activate signaling (e.g, lipid or calcium signaling) (Königshofer et al., 2008). These observations support the idea of Falcone et al. (2004) that the unsaturation level of lipid membranes also plays an important role in the plant's ability to tolerate high

temperatures, although other characteristics of plants membrane lipids are also likely to be important.

The study of the metabolome across heat stress and recovery, has given the possibility to deepen in the different dynamics of the metabolomic pathways in P. radiata, establishing the existence of key pathways in different moments of the stress response.

# Integrative Analysis of Physiological Response and Metabolome Dynamic

This work revealed that high temperature has a complex impact on cell function, suggesting that many and complex processes are involved in heat-resistance processes. However, the use of system-wide approach and integrative bioinformatics tools have allowed the understanding of the molecular basis, the identification of injury mediators, and the characterization of associated biomarkers. Thus, CKs, fatty acid metabolism and flavonoid and terpenoid biosynthesis were revealed as the most relevant pathways related to shorter-term and longer-term response clusters confirmed by multivariate analysis (**Figure 3A**), being ZR and iPA the key hormones that coordinate the different response processes (**Figure 3B**).

Some of the secondary metabolites with a higher loading in the shorter-term and longer-term response cluster included dihydrokaempferol, cis-1,2-dihydro-3-ethylcatechol and dihydromyricetin, which are all crucial elements in flavonoid metabolism. Flavonoids are a biologically and chemically diverse group widely represented in plants. Their diversity and multifunctionality demonstrate their importance in plants. Terpenoids are the other significant group of secondary metabolites identified in relation to the clusters, which are produced by a variety of plants and particularly in conifers (Zulak and Bohlmann, 2010). These volatile compounds are emitted by plants and play an important role in the interaction with their environment (Tholl, 2015). The best known and most studied group of terpenoids is the sesquiterpenoid plant hormone ABA, the central element in the plant stress response (McCourt et al., 2005; Raghavendra et al., 2010). Abscisic alcohol 11-glucoside, the glycosylated form of ABA, was also positively correlated to ZR as a key element in the network (**Figure 3B**). Both groups of secondary metabolites, flavonoids and terpenoids, are important in plant growth, development, and response against biotic and abiotic stress (Alder et al., 2012), as well as in plant adaptation to variable environmental conditions (Kliebenstein, 2004). Plasma membrane fluidity has been described to be an important temperature sensor in plants that appears to lie upstream of the unfolded protein response (Wu et al., 2012). Increased membrane fluidity appears to open calcium channels in the plasma membrane and the resultant inflow of calcium triggers signaling cascades, including an H2O2 burst (Königshofer et al., 2008), which activates the heat-stress response. In fact, fatty acids are the most prominent group of compounds revealed in the heat-stress response network, after CKs.

Integrated physiological and metabolome analysis across high-temperature treatments (Supplementary Movie S1) in P. radiata emphasizes the complex dynamics of the metabolome in response to heat stress and suggests the existence of a turning point (T3) at which P. radiata plants changed from an initial stress response (shorter-term response) to an acclimation one (longer-term response). The video highlights how the metabolites that are positively related to ZR increase drastically its accumulation in T5 and R. This is regulated by a complex interaction network that involves multiple pathways and groups of compounds where possible biomarkers of thermotolerance processes could be found.

### Evaluation of Proposed Thermotolerance Biomarkers in P. radiata

In plants, the concept of biomarker could be defined as "a characteristic that is objectively measured or evaluated as a predictor of plant performance" (Fernandez et al., 2016). The use of biomarkers originated from the field of medicine, but in plants in recent years, many authors have used metabolites as indicators for estimating plant performance under stress conditions (Quistian et al., 2011; Degenkolbe et al., 2013; Nam et al., 2015; Obata et al., 2015). One of main goals in this study is to find metabolic markers of high-temperature tolerance, potentially useful for P. radiata breeding programs.

L-phenylalanine, hexadecanoic acid (C16:0), and dihydromyricetin were confirmed as the three strongest biomarkers between the proposed candidates related to the results of this work. Their biological relevance, statistical strength in integrative analysis, and accumulation profile across the samples times analyzed validate their future use as thermotolerance biomarkers in P. radiata.

L-phenylalanine seems to be one of the clearest candidates. L-phenylalanine is crucial in numerous KEGG pathways which changed significantly during stress, particularly in shorter-term treatments (e.g., phenylalanine metabolism, phenylpropanoid biosynthesis, phenylalanine, tyrosine and tryptophan biosynthesis, glucosinolate biosynthesis, and biosynthesis of plant hormones). In addition, integrative analysis of metabolome and physiological measurements showed its relevance as a key compound in the network related to high-temperature response (**Figure 3B**). L-phenylalanine showed a direct correlation with the CKs associated with the trigger of response (negatively with ZR and positively with iPA). Old studies (Dedio and Clark, 1971; Deikman and Hammer, 1995) have already seen how the application of CKs stimulated flavonoids pathways in which L-phenylalanine is involved (flavonoids are synthesized by the phenylpropanoid metabolic pathway), such as the production of isoflavone and anthocyanins synthesis. More recently, Angelova et al. (2001) and Ali and Abbas (2003) have debated the impact of CKs on the accumulation of flavonol glycosides and Hamayun et al. (2015) found out how kinetin modulates isoflavone contents under salinity stress. Furthermore, PAL expression analysis (**Figure 7A**) confirms the beginning of the acclimation process in T5 when L-phenylalanine levels (**Figure 4A**) showed a significant decrease, validating the high value of this metabolite as a biomarker.

The second candidate biomarker was chosen within saturated fatty acids, given their importance in stabilizing membrane fluidity during heat stress (Grover et al., 2000; Larkindale and Huang, 2004; Wang et al., 2017). Hexadecanoic acid (also known as palmitic acid or C16:0) increases its accumulation in shorter

heat exposures, suffering a decrease in T5 and R, when the plant can be already acclimatized to stress (**Figure 5E**). Moreover, the low expression levels of both DES (**Figures 7B,C**) in shorter treatments underline the importance of saturated fatty acids, such as hexadecanoic acid, in the early response to high temperatures. This saturated fatty acid has already been studied by other authors such as Alfonso et al. (2001), Falcone et al. (2004), and Larkindale and Huang (2004) in relation to high-temperatures response. It has been even observed that the thermotolerance of an Arabidopsis mutant deficient in palmitic acid unsaturation is enhanced (Kunst et al., 1989). Furthermore, hexadecanoic is showed as an important element by the performed integrative analysis, appearing interconnected with the master regulator ZR in the network.

The last strong biomarker proposed is dihydromyricetin, also known as ampelopsin. Dihydromyricetin is a flavanonol included in anthocyanin biosynthesis inside flavonoid biosynthesis pathway. Flavonoids and anthocyanins have a photoprotective and antioxidant role (Dao et al., 2011; Zhang et al., 2012). Dihydromyricetin is consumed during anthocyanins synthesis, reducing its accumulation in shorter treatments and reaching the lowest levels in T3 (**Figure 6A**). This is validated by the increased expression of DFR across the exposure to high temperatures (**Figure 7D**), showing its maximum expression in T5 when dihydromyricetin begins to be accumulated again. Anthocyanins are produced under a variety of stresses such as UV-B (Valledor et al., 2012), low temperatures (Krol et al., 1995), or salinity (Wahid and Ghazanfar, 2006), although only a few studies have dealt with the effect of high temperatures on anthocyanin accumulation (Wahid and Close, 2007; Correia et al., 2014; de Leonardis et al., 2015). Dihydromyricetin was also prominent in multivariate and integrative analysis showing a negative correlation with iPA and BK. Cytokinin increases anthocyanin content and the transcript levels of PRODUCTION OF ANTHOCYANIN PIGMENT 1 (Das et al., 2012) which is according to the dropping levels of dihydromyricetin (anthocyanin precursor) observed in the first impact of the heat stress. Among many other biological functions, anthocyanins are considered the first line of defense against oxidative stress (Gould et al., 2002), scavenging oxygen radicals, and inhibiting lipid peroxidation (Chalker-Scott, 2002; Ling et al., 2007).

### CONCLUSION

This work shows that high temperatures induced a quick and dynamic change in the metabolome of P. radiata, in order to maintain homeostasis and facilitate survival. Integrative study of metabolome across high-temperature exposure and recovery plants allowed reaching a global view of molecular mechanism behind high-temperature response in P. radiata,

#### REFERENCES

Acharya, B. R., and Assmann, S. M. (2009). Hormone interactions in stomatal function. Plant Mol. Biol. 69, 451–462. doi: 10.1007/s11103-008- 9427-0

revealing complex interaction networks that involve CKs, fatty acid metabolism, and flavonoid and terpenoid biosynthesis being ZR and iPA the master regulators that trigger the global response. Additionally, novel potential thermotolerance biomarkers such as L-phenylalanine, hexadecanoic acid and dihydromyricetin have been proposed. However, these potential biomarkers need to be validated in further studies in P. radiata and their possible universality analyzed in other species.

### AUTHOR CONTRIBUTIONS

ME, MM, MC, and LV designed and performed the experimental work. ME processed the samples, integrated the datasets, and completed the metabolome analyses. LV and MM performed the mass spectrometry. GP realized the physiological analysis. ME and JP were involved in the preparation of samples. LV and JP helped with statistical analyses. ME and MM wrote the manuscript while MC, LV, JP, and GP supervised the manuscript. All authors read and approved the final manuscript.

### FUNDING

This publication is an output of the projects financed by the Spanish Ministry of Economy, Industry and Competitiveness (AGL2014-54995-P and AGL2016-77633-P), the Government of Principado de Asturias (GRUPIN14-055), FEDER funding through COMPETE (Project: UID/AMB/50017/2013), and by National Funds through the Portuguese Foundation for Science and Technology (FCT) within the Project PTDC/AGR-FOR/2768/2014. ME was supported by a fellowship from the Severo Ochoa Program (BP11117; Government of Principado de Asturias, Spain). JP was supported by a fellowship from the FPU (AP2010-5857; Ministry of Education, Spain). FCT (Fundação para a Ciência e a Tecnologia, Portugal) supported the fellowships of GP (SFRH/BPD/101669/2014). The Spanish Ministry of Economy and Competitiveness supported MM and LV through the Ramón y Cajal program (RYC-2014- 14981 and RYC-2015-17871, respectively). FCT/MEC (Portugal) co-funding by the FEDER, within the PT2020 Partnership Agreement and Compete 2020 provide financial support to Centre for Environmental and Marine Studies (CESAM – UID/AMB/50017).

# SUPPLEMENTARY MATERIAL

The Supplementary Material for this article can be found online at: https://www.frontiersin.org/articles/10.3389/fpls.2018.00485/ full#supplementary-material

Alder, A., Jamil, M., Marzorati, M., Bruno, M., Vermathen, M., Bigler, P., et al. (2012). The path from beta-carotene to carlactone, a strigolactone-like plant hormone. Science 335, 1348–1351. doi: 10.1126/science.1218094

Alfonso, M., Yruela, I., Almárcegui, S., Torrado, E., Pérez, A. M., and Picorel, R. (2001). Unusual tolerance to high temperatures in a new herbicide-resistant

D1 mutant from Glycine max ( L.) Merr. cell cultures deficient in fatty acid desaturation. Planta 212, 573–582. doi: 10.1007/s004250000421

plants via a cytosolic tyrosine:phenylpyruvate aminotransferase. Nat. Commun. 4:2833. doi: 10.1038/ncomms3833

**Conflict of Interest Statement:** The authors declare that the research was conducted in the absence of any commercial or financial relationships that could be construed as a potential conflict of interest.

Copyright © 2018 Escandón, Meijón, Valledor, Pascual, Pinto and Cañal. This is an open-access article distributed under the terms of the Creative Commons Attribution License (CC BY). The use, distribution or reproduction in other forums is permitted, provided the original author(s) and the copyright owner are credited and that the original publication in this journal is cited, in accordance with accepted academic practice. No use, distribution or reproduction is permitted which does not comply with these terms.

# Cyclotide Evolution: Insights from the Analyses of Their Precursor Sequences, Structures and Distribution in Violets (*Viola*)

Sungkyu Park 1 , Ki-Oug Yoo2 , Thomas Marcussen3 , Anders Backlund1 , Erik Jacobsson1 , K. Johan Rosengren4 , Inseok Doo5 and Ulf Göransson1 \*

1 Division of Pharmacognosy, Department of Medicinal Chemistry, Uppsala University, Uppsala, Sweden, 2 Department of Biological Sciences, Kangwon National University, Chuncheon, South Korea, 3 Department of Biosciences, Centre for Ecological and Evolutionary Synthesis, University of Oslo, Oslo, Norway, 4 School of Biomedical Sciences, The University of Queensland, Brisbane, QLD, Australia, 5 Biotech Research Team, Biotech Research Center of Dong-A Pharm Co Ltd., Seoul, South Korea

#### *Edited by:*

Luis Valledor, Universidad de Oviedo, Spain

#### *Reviewed by:*

Jesús Pascual Vázquez, University of Turku, Finland Diego Mauricio Riaño-Pachón, Institute of Chemistry, University of São Paulo, Brazil Monica Escandon, University of Aveiro, Portugal

> *\*Correspondence:* Ulf Göransson<EMAIL_ADDRESS>
#### *Specialty section:*

This article was submitted to Plant Systems and Synthetic Biology, a section of the journal Frontiers in Plant Science

> *Received:* 07 June 2017 *Accepted:* 17 November 2017 *Published:* 18 December 2017

#### *Citation:*

Park S, Yoo K-O, Marcussen T, Backlund A, Jacobsson E, Rosengren KJ, Doo I and Göransson U (2017) Cyclotide Evolution: Insights from the Analyses of Their Precursor Sequences, Structures and Distribution in Violets (Viola). Front.  Their cyclic cystine knot (CCK) motif makes them exceptionally resistant to thermal, chemical, and enzymatic degradation. By disrupting cell membranes, the cyclotides function as host defense peptides by exhibiting insecticidal, anthelmintic, antifouling, and molluscicidal activities. In this work, we provide the first insight into the evolution of this family of plant proteins by studying the Violaceae, in particular species of the genus Viola. We discovered 157 novel precursor sequences by the transcriptomic analysis of six Viola species: V. albida var. takahashii, V. mandshurica, V. orientalis, V. verecunda, V. acuminata, and V. canadensis. By combining these precursor sequences with the phylogenetic classification of Viola, we infer the distribution of cyclotides across 63% of the species in the genus (i.e., ∼380 species). Using full precursor sequences from transcriptomes, we show an evolutionary link to the structural diversity of the cyclotides, and further classify the cyclotides by sequence signatures from the non-cyclotide domain. Also, transcriptomes were compared to cyclotide expression on a peptide level determined using liquid chromatography-mass spectrometry. Furthermore, the novel cyclotides discovered were associated with the emergence of new biological functions.

Keywords: cyclotide evolution, viola phylogeny, sequence signature, cyclotide precursor, neofunctionality, novel cyclotide, precursor domain

# INTRODUCTION

Cyclotides are proteins of ∼30 amino acid residues that are characterized by the cyclic cystine knot (CCK) motif (Craik et al., 1999; Burman et al., 2014). The CCK motif consists of six conserved cysteines that form three disulfide bonds and a head to tail cyclic backbone (**Figure 1A**). The cyclotides have been classified into two main subfamilies, the Möbius and the bracelets, based on a single structural trait: the presence or absence of a conceptual 180◦ twist in the cyclic backbone caused by a conserved cis-Pro residue in loop 5 (Craik et al., 1999; **Figure 1B**). Aside from the CCK, two loops (defined as sequences between adjacent cysteines) have high sequence

similarity between subfamilies (loop 1 and 4), whereas loops 2 and 3 are conserved only within individual subfamilies.

VocA (δ). In the precursor of Cter M, the cyclotide domain replaces the albumin-1b chain.

The discovery of other cyclotides has created a need for a more versatile classification system (Ireland et al., 2006; Nguyen et al., 2013; Ravipati et al., 2015). These varieties include socalled hybrid cyclotides that exhibit sequence characteristics of both the Möbius and bracelet subfamilies (Daly et al., 2006), as well as minor subfamily known as the trypsin inhibitors originating from gourd plants (Hernandez et al., 2000). They contain the CCK motif but do not otherwise exhibit any sequence similarity with the other subfamilies. In addition, linear cyclotide derivatives that exhibit sequence similarity with conventional cyclotides but lack their cyclic backbone have been reported (Ireland et al., 2006; Nguyen et al., 2013). The high sequence diversity of the cyclotides appears to be due to natural selection in angiosperms, the flowering plants, but little is known about the evolutionary mechanisms underpinning the corresponding selection processes or the evolutionary background of cyclotide diversity.

Cyclotides and the CCK motif have only been found in angiosperms, but proteins having one of their two defining structural motifs—cyclic peptides/proteins without the cystineknot (Trabi and Craik, 2002; Arnison et al., 2013) or linear proteins with a cystine-knot (Zhu et al., 2003)—are found in a wide range of organisms across all kingdoms of life. In angiosperms, the occurrences of cyclotides differ between "basal" angiosperms, monocots and eudicots (**Figure 1C**): "linear cyclotides," i.e., peptides that exhibit sequence similarity with cyclotides but lack their head-to-tail cyclic structure, are prevalent in both monocots and eudicots, but true cyclotides have been found only in eudicots (Mulvenna et al., 2006; Zhang et al., 2015a). However, neither linear nor true cyclotides have yet been found in the "basal" angiosperms. It has therefore been proposed that linear cyclotides are ancestral (more primitive) to the true cyclic cyclotides (Mulvenna et al., 2006; Gruber et al., 2008).

To date, cyclotides have been discovered in the eudicot families of Rubiaceae (coffee) (Gran, 1973a), Violaceae (violet family) (Schöpke et al., 1993; Claeson et al., 1998), Fabaceae (legume family) (Poth et al., 2011a), Solanaceae (potato) (Poth et al., 2012), and Cucurbitaceae (cucurbit) (Hernandez et al., 2000), as well as in the monocot family Poaceae (grass family) (Nguyen et al., 2013). Cyclotides appear to have functions in host defense because they exhibit insecticidal (Jennings et al., 2001), anthelmintic (Colgrave et al., 2008a), antifouling (Göransson et al., 2004), and molluscicidal (Plan et al., 2008) activities.  This plethora of activities and the stability of the CCK motif make them of interest for drug development (Northfield et al., 2014).

Many of these activities appear to be due to cyclotides' ability to interact with and disrupt biological membranes (Colgrave et al., 2008b; Simonsen et al., 2008; Burman et al., 2011; Henriques et al., 2011). The membrane disruption is mediated by physicochemical interactions between cyclotides and the lipid membrane, and is governed by the distribution of lipophilic and electrostatic properties over the molecular surfaces of the cyclotides. We recently developed a quantitative structureactivity relationship (QSAR) model for these interactions (Park et al., 2014). However, the relationships between cyclotide sequence diversity, evolutionary selection, and the functions of the cyclotides in planta remain unknown.

Cyclotides are expressed as precursor proteins, which undergo post-translational processing including enzymatic cleavage and subsequent cyclization (Jennings et al., 2001; Harris et al., 2015). The multi-domain architecture of these precursor proteins varies slightly between different types of cyclotides and plant families, but in sequential order from the N- to the Cterminus, they generally feature the following domains: an endoplasmic reticulum (ER) targeting signal, an N-terminal propeptide (NTPP), an N-terminal repeat (NTR), the cyclotide domain (CD), and finally a C-terminal tail (CTR) (**Figure 1D**). In some cases, the modular domains NTR, CD, and CTR are repeated more than once. The cyclotides have been suggested to co-evolve with asparaginyl endopeptidase (AEP) because of its suggested role in cyclization (Mylne et al., 2012). Moreover, the divergent evolution of cyclotides from ancestral albumin domains was suggested based on the architecture of cyclotide precursors found in the Fabaceae plant family (Nguyen et al., 2011; Poth et al., 2011b). However, the relationship between the precursor proteins' architecture and sequences and the evolutionary selection of cyclotides is still unknown.

To date, cyclotide and precursor sequences have been most extensively explored in the family Violaceae Batsch. (Malpighiales), and especially in the genus Viola L. (Burman et al., 2015). The Violaceae are a medium-sized family including ∼1,100 species worldwide. The phylogeny of the Violaceae has recently been inferred from chloroplast and nuclear markers (Tokuoka, 2008; Wahlert et al., 2014), and its systematics has been revised accordingly; currently ∼30 genera are accepted in nomenclature (Wahlert et al., 2014). Viola is the largest genus in the family, with 580–620 species, representing over 50% of all known species (Ballard et al., 1998; Yockteng et al., 2003; Marcussen et al., 2012; Wahlert et al., 2014). Viola is distributed all around the world in temperate regions and at high elevation habitats in the tropics. The genus is old (∼30 million years) and comprises at least 16 extant lineages, referred to as sections, with a complex, reticulate phylogenetic history owing to allopolyploidy (Marcussen et al., 2012, 2015). The species included in this study belong to four northtemperate sections, the diploid sect. Chamaemelanium Ging. (V. canadensis L., V. orientalis W.Becker) and the three allotetraploid sections Melanium Ging. (V. tricolor L.), Plagiostigma Godr. (V. albida Palibin. var. takahashii (Nakai) Kitag., V. mandshurica W.Becker, V. verecunda A.Gray) and Viola (V. acuminata Ledeb.).

In the current study, we explore cyclotide evolution using an integrated approach, exploiting transcriptomics and peptidomics to analyze the sequences of cyclotide precursors and the expression of cyclotides in Viola in light of the phylogeny of the genus. In particular, full precursor sequences are used to obtain insights into the evolutionary history of the cyclotides, and we connect the evolution of new mature cyclotides with the emergence of new functions.

# MATERIALS AND METHODS

#### Collection of Violets

Violets were collected at the sites indicated in **Table 1**, and those sites are their natural habitats. When collected, the plant individuals were in adult stage; however, the exact ages of those plant individuals were not determined. The plant vouchers were deposited at the Kangwon National University herbarium; Viola albida var. takahashii (KWNU93021), V. mandshurica (KWNU93022), V. orientalis (KWNU93023), V. verecunda (KWNU93024), and V. acuminata (KWNU93025). All plant material was collected on September 23, 2014.

#### Sample Collection, RNA Isolation, and RNA Sequencing

For the five Viola species—V. albida var. takahashii, V. mandshurica, V. orientalis, V. verecunda, and V. acuminata—, total RNA was sequenced by Next Generation Sequencing (NGS), outsourced to Macrogen Inc. (Seoul, South Korea). For each Viola species, the RNA sample was prepared from one plant individual; also, the plant tissues were pooled from all of the plant's major organs, i.e., the roots, stems, flowers, and leaves. The collected tissues were immediately frozen in liquid nitrogen, and directly extracted by RNeasy Plant Mini Kit (Qiagen, Hilden, Germany) according to the manufacturers' protocols. Quality and quantity of RNA were measured using an Agilent 2100 Bioanalyzer (Agilent Technologies, Santa Clara, CA, USA) with an RNA Integrity Number (RIN) and rRNA ratio. The measured RIN values (rRNA ratios) are: 6.9 (1.2)

for V. albida var. takahashii, 7.4 (1.6) for V. mandshurica, 7.2 (0.2) for V. orientalis, 6.4 (0.9) for V. verecunda and 5.9 (0.7) for V. acuminata. For mRNA library preparation, a TruSeq RNA preparation kit was used according to the manufacturer's instructions (Illumina, San Diego, U.S.A.). Briefly, the poly-A containing mRNAs were purified using poly-T oligo-attached magnetic beads. The purified mRNAs were fragmented into short sequences by use of divalent cations. Using these short sequences as templates, the first-stranded cDNA was synthesized using random hexamers. A second-strand cDNA was then synthesized using DNA polymerase I and RNase H. The synthesized cDNA went through an end repair process, the addition of a single "A" base, and then ligation of the adapters. The PCR was performed to enrich the selected DNA sequences, and then those selected sequences were sequenced using Illumina HiSeq 2000 Sequencing System (Illumina, USA) that generates paired-end reads with 2 × 100 base pairs (bp) read length.

#### Sequencing Data Analysis and Assembly

FASTQC (version 0.11.3) was used to determine the quality of RNA sequencing data (www.bioinformatics.babraham.ac.uk/ projects/fastqc/). Reads were cleaned using Trimmomatic (v. 0.32; Bolger et al., 2014), and the sequences with Phred score ≥33 and a minimum length of 36 bp were retained for assembly. De novo assembly of these processed reads was performed using the Trinity RNA-Seq assembly (release 17.07.2014; https:// sourceforge.net/projects/trinityrnaseq/) with the default setup, which allows the identification of cyclotide-coding transcripts larger than 200 bp. Also, the reads were assembled separately for the individual species. We summarized the detailed information on the de novo transcriptome assembly for those individual species in Supplementary Table 1. The transcript abundance, Fragments Per Kilobase of transcript per Million mapped reads (FPKM; Trapnell et al., 2010), was calculated using Trinity based on RSEM algorithm (Li and Dewey, 2011).

#### Identification of Precursor Sequences from Transcriptome

The assembled transcriptomes were searched for similar sequences to the cyclotides from Cybase (cutoff date: 3.02.2014; http://www.cybase.org.au/) by the standalone NCBI-blast+ service (2.2.28) (tblastn, E-value cutoff: 50) in the Ugene software package (Okonechnikov et al., 2012). In parallel, a motif search was done using PROSITE for cyclotide precursor sequences [motifs PS60009 and PS60008 for Möbius and bracelet cyclotides, respectively (Sigrist et al., 2002)]. To assist sequence identification, the Fuzzpro service of EMBOSS (v. 5.0.0) (Rice et al., 2000) was used. The combined PROSITE and blast results were filtered to remove duplicates. The result file was further processed by manual inspection after clustal omega (v.1.2.0) alignment. In the manual inspection, we assumed that the sequence is a cyclotide precursor if the sequence contains the six conserved cysteines aligned with previously known cyclotides, and if the conserved Nterminal domain showed sequence similarity to known cyclotide precursors.

#### Collection of Precursor Sequences

In total, 312 (= 157 + 155) precursor sequences were utilised in this study: (i) 157 (=138 + 19) precursor sequences were identified and collected from transcriptome analyses conducted in the current study. Among them, 138 sequences were collected from the five Viola species (i.e., Viola albida var. takahashii, V. mandshurica, V. orientalis, V. verecunda, and V. acuminata) by the RNA sequencing in this study. The remaining 19 sequences were collected the transcriptomic data of V. canadensis obtained from the 1kp-project (www.onekp.com). (ii) Another 155 (=126 + 29) precursor sequences were collected from recent published studies. Among them, the 126 sequences were collected from other seven Viola species (i.e., V. baoshanensis, V. odorata, V. uliginosa, V. adunca, V. tricolor, V. biflora, and V. pinetorum), and the 29 sequences from other Violaceae genera, i.e., Gloeospermum, Melicytus, and Pigea. The name of Viola species where the precursor sequences collected from and the corresponding references are listed in the **Table 2**, and the sequence alignment of those 312 precursors are recorded in the Supplementary Data 1.

#### Nomenclature of Cyclotides and Precursors

Each cyclotide sequence was assigned a tripartite name. The first part is derived from the Latin binomial of the plant in which the corresponding precursor was found, the third part specifies the molecular species of the cyclotide, and the second part specifies the rank of the cyclotide among all the cyclotides having that particular molecular species that were identified in that particular

, *V. verecunda* A.Gray [23]a

, V.

The numbers in parentheses after the species name show the number of precursor sequences derived from that species. Species names marked in bold were used for transcriptomic assay of the current study. Sections are indicated within Viola. Within sect. Plagiostigma, V. mandshurica, and V. albida var. takahashii belong in subsect. Patellares (Boiss.) Rouy & Foucaud, V. verecunda in subsect. Bilobatae (W.Becker) W.Becker and V. baoshanensis in subsect. Diffusae (W.Becker) Chang. Within sect. Viola, V. odorata belongs to subsect. Viola and V. acuminata, V. adunca, and V. uliginosa in subsect. Rostratae (W.Becker) W.Becker. No subsections are applicable to sect. Chamaemelanium. aTranscriptomes from the current study, V. mandshurica, V. albida var. takahashii, V. verecunda, V. acuminata and V. orientalis; bV. baoshanensis (Zhang et al., 2009, 2015b); cV. odorata (Dutton et al., 2004; Ireland et al., 2006); dV. uliginosa (Slazak et al., 2015); eV. pinetorum and V. adunca (Kaas and Craik, 2010); fV. tricolor (Mulvenna et al., 2005; Hellinger et al., 2015); gV. biflora (Herrmann et al., 2008); hTranscriptome data for V. canadensis were obtained from the 1kp project (www.onekp.com), and the cyclotide precursor sequences were determined in this work; iG. blakeanum and G. pauciflorum (Burman et al., 2010); jM. ramiflorus (Trabi et al., 2009); kP. floribunda (as H. floribundus) (Simonsen et al., 2005).

plant species (Supplementary Table 2). The Latin binomials of the violet species considered in this work are Viola albida var. takahashii (valta), Viola mandshurica (viman), Viola orientalis (vorie), Viola verecunda (viver), and Viola acuminata (vacum). Thus, the cyclotide named vacum2-HS4 is the second cyclotide of the HS4 molecular species derived from V. acuminata. The molecular species are named after their NTR signature sequences, i.e., the identities of the two residues at the consensus positions 11 and 12. The precursor sequences were also assigned tripartite names in a similar way to the cyclotides: the first part of the name is derived from the Latin binomial of the species in which the sequence was discovered, the second specifies the numerical rank of the sequence (which is independent of that assigned to the cyclotides), and the third specifies the molecular species of the precursor. Any cyclotide or precursor that had been named in an earlier publication was assigned the same name in this work. However, for previously unnamed precursors, we added the prefix "prc" to the tripartite name to help readers distinguish between precursor and cyclotide sequences. Thus, prc-viul A is the name of the precursor sequence of Viul A from Viola uliginosa Bess.

#### Sequence Alignment for Phylogenetic Analysis

The cyclotide precursor sequences from transcriptome were aligned before the phylogenetic analysis i.e., construction of Bayesian phylogenetic tree, maximum parsimony and splits networks. A total of 92 precursor sequences (two sequences from each of the 46 molecular species) were prepared as DNA sequence for the sequence alignment, and the sequences include only three domains of the precursor, i.e., NTPP, NTR, and cyclotide domains. To guide the alignment of DNA sequences of the precursors, they were translated into protein sequences and aligned independently for each molecular species using Clustal in BioEdit v.7.2.5 (Hall, 1999), and those alignments were in turn combined and realigned (i.e., keeping the indel positions from the alignment of each of the molecular species). The resulting protein-guided alignment of nucleotide sequences was then subject to manual adjustments within reading frames. The aligned DNA sequences are shown in the Supplementary Data 2.

#### The Construction of Bayesian Phylogenetic Tree

Nucleotide substitution model was selected based on the AICc criterion using JModelTest v.2.1.10 (Guindon and Gascuel, 2003; Darriba et al., 2012). A Bayesian phylogenetic analysis was conducted in BEAST v.1.7.4 (Drummond and Rambaut, 2007; Drummond et al., 2012). The analysis file was set up in BEAUti (part of the BEAST package) with the following priors: GTR+G as substitution model with empirical frequencies and four gamma categories, a lognormal relaxed clock prior with rate set to 1.0, and a Yule tree prior (birth-only process). Two MCMC chains were run for 50 million generations each, using a BEAGLE library, and parameters logged every 10,000 generations. We checked the two chains for proper mixing and convergence (i.e., ESS >500) in Tracer v.1.6.0 (http://tree.bio. ed.ac.uk/software/tracer/), removed a visually determined burnin of 1 million generations from each and merged the two chains in LogCombiner v.1.7.4 (part of the BEAST package), and summarized the data in a maximum credibility tree with mean node heights using TreeAnnotator v.1.7.4 (part of the BEAST package). The resulting tree was visualized and edited in FigTree v.1.4.3 (http://tree.bio.ed.ac.uk/software/figtree/). The input file for the BEAST and the resulting phylogenetic tree are found in the Supplementary Data 3, 4, respectively.

#### The Construction of Splits Networks

In order to visualize the data, neighbor splits-networks of uncorrected P distance were produced, both for the nucleotide alignment and for the translated alignment, using SplitsTree (Huson and Bryant, 2006).

### Extraction of Plant Material

Between 250 and 500 mg of dried plant material was homogenized and incubated overnight in 6 ml of 60% acetonitrile in water containing 0.05% triflouoroacetic acid (TFA) to extract cyclotides. Extracts were lyophilized and then dissolved in 2.5 ml of solvent A (Milli-Q H2O with 0.05% TFA). Redissolved extracts were then subjected to gel filtration using PD-10 columns (GE Healthecare) according to the manufacturer's instructions to remove small molecules. The high molecular weight fractions were collected, lyophilized and then dissolved in solvent A to a concentration proportional to the original amount of extracted material (2 µl/1 mg) for LC-MS analysis.

# LC-MS

The samples were analyzed using ultra performance liquid chromatography coupled to quadrupole time-of-flight mass spectrometry (nanoAcquity UPLC/QTof Micro; Waters, Milford, MA). Samples were eluted using a gradient of acetonitrile (1 to 90% over 50 min) containing 0.1% formic acid. A nanoLC column (Waters BEH, 75µm (i.d) × 150 nm) operated at 0.3 µl/min flow rate was used. The capillary temperature was set at 220◦C and the spray voltage at 4 kV. The mass-to-charge (m/z) range was set from 1,000 to 2,000. LC-MS chromatogram and MS spectra were analyzed with the help of MassLyxn V4.1 (Waters, Milford, MA).

#### Reduction and Alkylation

Dried peptide extracts were reduced in a buffer containing 0.05 M Tris-HCl, pH 8.3, 4.2 M guanidine-HCl, and 8 mM DTT. The extraction solutions were incubated at 37◦C for 2 h in the dark, followed after O2 removal with Nitrogen gas. The reduced peptides were then further alkylated in a buffer containing 0.2 M Tris-HCl, pH 8.3, and 200 mM iodoacetamide for 1 h.

#### Identification of Cyclotides from the Plant Extracts

On the original LC-MS chromatogram, all chromatographic peaks were manually investigated if they include a cyclotide-like mass spectrum. We assumed that peaks stemmed from cyclotidelike substances if their masses fell within the range of 2,700–3,300 Da, as deconvoluted from their doubly- and triply charged ions. Then, the presence of three disulfide bonds was used to support their identification as cyclotides: if peaks in the extract showed an increase in mass by 348.18 Da, we considered that the peak stemmed from a true cyclotide (**Figure 2**).

# Matching Transcriptomic Data to Cyclotides on the Protein-Level

Presence of possible cyclotides identified in precursor sequences from transcriptome analyses was determined as follows: putative cyclotide sequences were listed and their monoisotopic mass were calculated. The listed cyclotide sequences were comprised of N-terminal residues ranging from [−3, 2] on the consensus position of their precursor sequences. The C-terminal residues were found differently between cyclic cyclotides and linear cyclotides, i.e., the highly conserved N/D residue located at loop 6 for cyclic cyclotides, and the residue located next to stop codon in their mRNA sequences for linear cyclotides. These calculated monoisotopic mass were then compared to the observed monoisotopic mass from the identified cyclotides in plant extracts using LC-MS. We regarded that the transcriptomic cyclotides were expressed in the protein-level if the difference of the monoisotopic mass <0.40 Da.

### Determination of Cyclotides' Abundance Levels in the Proteome

For each cyclotide-like substance, the chromatographic peaks were identified together with their own retention time from the original LC-MS chromatogram. These chromatographic peaks were further investigated in relation with their mass spectral peaks to estimate the signal intensity (SI) of the cyclotide-like substance. The SI is estimated as a summed signal intensity of triple-charged mass spectral peaks. According to the SI, we assigned the cyclotide abundance into three levels, i.e., the abundance level is: low if SI <250, high if SI > 1,000, and medium if 250 < SI < 1,000.

#### The Calculation of Molecular Descriptors

The physicochemical properties of selected cyclotides (Park et al., 2014), i.e., the total lipophilicity and the exposure ratio, were calculated using scientific vector language (SVL), implemented in MOE 2012 (Chemical Computing Group Inc., Montreal, Canada).

#### Accession Numbers

The accession numbers of Sequence Read Archive (SRA) and Transcriptome Shortgun Assembly (TSA) for transcriptome of the five Viola species are: V. acuminata (TSA: GFWD00000000, SRA: SRR5320546), V. orientalis (TSA: GFXR00000000, SRA: SRR5322130), V. albida var. takahashii (TSA: GFWC00000000, SRA: SRR5320531), V. mandshurica (TSA: GFWG00000000, SRA: SRR5320533), and V. verecunda (TSA: GFWF00000000, SRA: SRR5322180).

# RESULT AND DISCUSSION

The family Violaceae and the genus Viola are an excellent system for studying the evolution of cyclotides. All investigated species express large numbers of these proteins (Burman et al., 2015) as well as their precursor sequences. Furthermore, the specieslevel phylogeny of Viola (Marcussen et al., 2012, 2015) and chloroplast phylogeny of Violaceae (Wahlert et al., 2014) are wellunderstood as a result of recent studies. In the current study, 312 (= 157 + 155) sequences were analyzed in total to investigate the evolution of cyclotides. Among them, 157 precursor sequences were discovered from the transcriptomes of the six Viola species, i.e., V. albida var. takahashii, V. mandshurica, V. orientalis, V. verecunda, V. acuminata, and V. canadensis (listed in the Supplementary Table 3). Those discovered 157 sequences were pooled with the other 155 precursor sequences recently found from the other Viola species (i.e., V. baoshanensis, V. odorata, V. uliginosa, V. adunca, V. tricolor, V. biflora, and V. pinetorum) and from other Violaceae genera (i.e., Gloeospermum, Melicytus, and Pigea).

isotopic peaks. The triple-charged monoisotopic mass of the cyO2 (marked as [cyO2]+3 in the box) is observed as 1047.1Da, and the double-charged monoisotopic mass ([cyO2]+2 ) is observed as 1570.2Da. In each mass spectrum, the distances between isotopic peaks are observed with ∼0.3Da for [cyO2]+3 , and 0.5Da for [cyO2]+2 . The monoisotopic mass of the cyO2 is 3138.4Da (calculated) and 3138.3Da (observed), and the mass difference is 0.1Da. (B) The monoisotopic mass of the cyO2 after the alkylation reaction is 3486.6Da (calculated) and 3486.9Da (observed), and the mass difference is 0.3Da. The calculated mass is based on the triple-charged monoisotopic mass of cyO2 (1163.3Da), and its related mass spectrum is shown in the box marked [cyO2-Alk]+3 .

#### The Sequence Signature of the Prodomain Can Be Used for the Classification of Cyclotides

The structural classification of cyclotides into Möbius cyclotides, bracelets, and hybrids thereof based on the structure and sequence of the mature cyclotide is here replaced by a classification based on the sequences of the cyclotide precursors, including both the prodomain (i.e., the NTPP and NTR domains) and the cyclotide domain. We focused on these three domains (i.e., the NTPP, NTR, and cyclotide domains) because they are present in all of the precursors that have been completely sequenced. The ER domain was not included in the analysis because ER domain sequences often vary strongly with the quality of the sequencing data. Other domains (i.e., the CTR domain and repeats of the NTR and cyclotide domains) are not present in all precursors, and were therefore also excluded.

Precursors were classified by their sequence signatures, i.e., the patterns of insertions and deletions (indels) and conserved sequences in the prodomains (NTPP and NTR). The Nterminal cleavage site of the cyclotide domain was defined as position 0, and a large indel region was detected upstream in the NTPP at positions [−56, −38] (**Figure 3**). Within this indel, the insertion region features sequence variations with minor gaps, and the deletion region has definite sequence gaps [−56, −54], [−50, −38]. Interestingly, the insertions coincide with cyclotide domains of archetypical bracelet cyclotides (e.g., cycloviolacin O2), while deletions coincide with cyclotide domains of archetypical Möbius cyclotides (e.g., kalata B1). Based on these observations, we suggest that these indels in the NTPP domain can be used as a criterion for classifying precursors into the Möbius and bracelet lineages. The term lineage is used instead of subfamily in order to avoid confusion with the structural classification of the cyclotides.

FIGURE 3 | Classification of cyclotide precursors based on the sequence signature. Consensus sequences of the molecular species identified in this work. Based on the sequence signatures in the prodomain, 249 cyclotide precursors were classified into the Möbius and bracelet lineages, 13 molecular series and 46 molecular species (14 Möbius and 32 bracelet species). Species are named for their signature sequences, and numbers within brackets show the numbers of complete and total (complete and partial) precursor sequences for each species. The left box shows the NTPP indel region [−56, −38], which contains defining sequence gaps at [−56, −54], and [−50, −38] in the Möbius lineage. Within a lineage, specific prodomain sequence-traits are associated with different structural subfamilies. In the Möbius lineage, the sequence (Y∧/F∧/H)−9 -(S∧/A/Y)−8 at position [−9, −8] is associated with cyclic cyclotide precursors (where "∧" indicates high occurrence of the labeled residue), while Y−9 -Y−8 is associated with linear cyclotide precursors. Precursors of structural hybrids exhibit have an insertion at [−32, −31], whereas the archetypical Möbius precursors have a deletion at this position. In the bracelet lineage, the sequences (H∧/N∧/S/T/G/K/P)−9 -(L/N/S/F/A)−8 and (Q/E/P/K)−9 -(D/N)−8 are associated with cyclic and linear precursors, respectively. In addition, linear precursors contain a deletion that is flanked by two defining insertions, (P∧/A/L)−49-(N∧/A)−48 and (D/E)−39, in the NTPP at [−49, −48] and [−39]. Conserved amino acids are indicated by their one-letter codes, and variable residues are represented by symbols indicating the physicochemical properties of their side chains: hydrophilic (∼), hydrophobic (\$), flexible (G), rigid (P), disulfide-forming (C) and indels (–). Hydrophilic residues are further divided into positively charged (+), negatively charged (=) and uncharged hydrophilic (\*); hydrophobics (\$) are divided into aromatic (#) and alkyl (<) groups. Numbers are used as consensus symbols for multiple physicochemical categories.

Within each lineage, precursors can be further divided by analogy to their structural classifications. That is to say, it is possible to identify sequence traits associated with linear and cyclic cyclotides within each lineage, and sequence traits associated with archetypical and hybrid subfamilies within the cyclic cyclotides of the Möbius lineage. Signature sequences of residues [−56, −38] were used to subdivide the bracelets, and residues [−32, −31] for the Möbius lineages, together with residues [−9, −8] of the NTR. Within the bracelet lineage, the sequence signatures associated with linear and cyclic cyclotide precursors differ in two ways. First, linear cyclotide precursors contain a signature sequence in the NTPP domain defined by two sequence insertions at [−49, −48] and [−39], together with a sequence deletion between those two sequence insertions. The signature sequences at these two locations are (P∧/A/L)−49 - (N∧/A)−48 and (D/E)−39, where the symbol ∧ denotes a residue that occurs most frequently at the indicated position. Conversely, precursors of cyclic peptides have variable sequences with minor gaps in the [−47, −40] region. Second, whereas positions [−9, −8] in the NTR of precursors of linear cyclotides are conserved, having the signature sequence (Q/E/P/K)−9 -(D/N)−8 , the precursors of cyclic cyclotides show higher diversity, having sequences of (H∧/N∧/S/T/G/K/P)−9 -(L/N/S/F/A)−8 . Different signature sequences were also identified between the three structural subfamilies of the Möbius lineage, i.e., linear, hybrid and archetypical Möbius cyclotides. First, while precursors of both linear and hybrid cyclotides contain a sequence insertion in the NTPP domain at [−32, −31], precursors of archetypical Möbius cyclotides contain a deletion at this position. Second, while the sequence at [−9, −8] in the NTR domain of linear cyclotides is highly conserved (Y−9 -Y−8 ), it is more variable in the precursors of cyclic cyclotides (Y∧/F∧/H)−9 -(S∧/A/Y)−8 .

On the basis of these findings, i.e., the combined sequence signatures of the NTPP and NTR domains and the conservation of residues with similar physicochemical properties in the cyclotide domain, two new classification orders were defined and used to classify the members of each lineage. These new orders were termed the molecular species and molecular series. Precursors with the same sequence signatures in both the NTPP and NTR domains are assigned to the same molecular species, while precursors that only share signature sequences in the NTR domain are assigned to the same molecular series.

Of 283 precursor sequences, 80 were classified into the Möbius lineage and 181 into the bracelet lineage (**Figure 3**; see Supplementary Figure 1). Of these 261 sequences, 249 were classified into 13 molecular series and 46 molecular species. Thus, 78 sequences representing the Möbius lineage were classified into five molecular series and 14 molecular species, while 171 sequences representing the bracelet lineage were classified into 8 molecular series and 32 molecular species. The remaining precursor sequences (34 of the initially examined 283, or 12.0% of the total) could not be grouped with any other sequences, and not be classified with molecular species (Supplementary Figure 2). Most of these sequences (31/283, or 10.9% of the total) were either partial sequences or lone unique sequences. Only three of all the precursor sequences (3/283) did not exhibit features enabling their classification into a given lineage on the basis of their NTPP domain sequences.

Informal hierarchical ranks were then assigned to this system for classifying precursors, with lineage being the highest ranking classification, followed by molecular series and then molecular species. We further investigated the evolutionary relevance of this classification system by performing the phylogenetic analysis using the full DNA precursor sequences (**Figure 4**). We assumed that if the sequence signatures are highly conserved in the course of the evolution, the phylogenetic relationship of the precursor sequences has the hierarchical ranks as classified based on the sequence signatures. In this evaluation, we adopted a Bayesian phylogenetic approach using BEAST, because substitution model-based methods are in most cases better than Maximum Parsimony or distance-based methods at handling homoplasy, which is prominent among the cyclotide precursor sequences (Splits network shown in Supplementary Figure 3), and because BEAST sets a prior on the branch lengths that is particularly suitable also for analysis of short DNA sequences, such as cyclotide precursors.

The two lower hierarchical ranks (i.e., molecular species and series) are consistently recovered in the phylogeny, but not the upper hierarchical rank (i.e., molecular lineage), as indicated by the low posterior support values at the base of the phylogeny (**Figure 4**). Most of precursor sequences belong to the same molecular species are grouped as one monophyletic clade (80%, 74/92) or to the different clades but more closely than the precursor sequences derived from the different molecular series. Among 13 molecular series, the molecular species belong to the four molecular series (i.e., HF, YA, FA, and HS) are not monophyletic. Only one molecular species (DI1) is exceptionally grouped into Möbius lineage in the phylogenetic tree, even though the DI1 is classified as the bracelet lineage based on the sequence signature. In the phylogenetic analysis, we selected a total of 92 precursor sequences by the random selection of two sequences from each of the 46 molecular species. We assumed that such precursor selection would be enough to show a phylogenetic support for the classification system based on the sequence signature approach, because we randomly selected two sequences from each molecular species, and these sequences were mostly paired or grouped into monophyletic clades in accordance with the classification at the level of molecular species and series.

#### The Distribution of Cyclotide Precursors Reflects the Phylogeny of the Genus *Viola*

Precursor sequences were compared to the established phylogenetic relationships between the four infrageneric sections of Viola (i.e., sects. Melanium, Plagiostigma, Viola, and Chamaemelanium), and between the genus Viola and the other Violaceae genera. The presence of a molecular species in different Viola taxa, such as species or sections, indicates that the molecular species evolved prior to the most recent common ancestor of these taxa. However, such inferences in Viola are complicated by the network-like phylogeny of the genus, owing to repeated ancient events of allopolyploidy (Marcussen et al., 2015). Hence, out of four sections in the current analysis, three (i.e., sects. Melanium, Plagiostigma, Viola) are allotetraploids originated by independent hybridizations between the same two parental lineages MELVIO and CHAM around 15 million years ago. The last section, sect. Chamaemelanium, is diploid and descends from the CHAM lineage. Because the parental MELVIO lineage is now extinct, it is impossible to infer with certainty the distribution of the molecular species in the common ancestor of these parental lineages, i.e., the CHAM and MELVIO lineages, especially without diploid outgroups from sister sections (e.g., Rubellium and Andinium). Gene flow by introgression between species potentially occurs between closely related species only, i.e., within the same subsection, and can be ruled out for this dataset (e.g., Marcussen et al., 2015).

The current analysis revealed several points about the distribution of cyclotide precursors within the genus Viola. Firstly, the classification into infrageneric sections is reflected in the distribution of the cyclotide precursors. Some molecular species occurred sporadically or commonly across the four infrageneric sections: 2% of molecular species were found across all sections, and 5% of the studied precursor sequences belong to those molecular species (**Figure 5**; see Supplementary Figure 4). Also, 17% of molecular species (34% of precursor sequences) were found across three sections, and 45% of molecular species (41% of precursor sequences) were found in at least two sections. This indicates that the molecular species likely originated both from the common ancestor of Viola sections and from the hybridization between the parental lineages. Also, there are

some molecular species that could have occurred by the genetic changes from Viola speciation. It should be emphasized that the degree of concurrence is likely higher at the genomic level than at the transcriptomic level estimated from the current study, because all genomic presence of the molecular species can not be captured by the transcriptome. Some molecular species might not

signature.

sequence signature. DI1 is grouped into Möbius lineage in the phylogenetic tree (pp ≥ 0.70), even though it is classified as bracelet lineage based on the sequence

cyclotide molecular species and molecular series found in four different Viola sections. (A) The values shown in the outer rings shows results obtained by considering numbers of molecular species, and while those in the inner rings show results obtained based on numbers of precursor sequences. (B) The values shown in the outer rings shows results obtained by considering numbers of molecular series, and while those in the inner rings show results obtained based on numbers of precursor sequences. Those results are based on data for complete sequences. The result based on both complete and partial sequences is shown in Supplementary Figure 4. (C) Phylogenetic relationships between sections and genera of Violaceae used in the analyses conducted within this work. Viola Sections are abbreviated as follows: Plagiostigma as PLA, Viola as VIO, Melanium as MEL, Chamaemelanium as CHA, Rubellium as RUB, and Andinium as AND. Dotted lines indicate the complex ancestries of the three allotetraploid sections Melanium, Plagiostigma, and Viola, all derived from hybridization between the CHAM and MELVIO lineages 15–20 Ma ago. Genera and sections studied using transcriptomic methods are indicated with asterisks (\*). The total number of species within Viola is estimated to be 580–620, most of which (61–65%) belong to these four sections. The first ancestor of the genus (α) is dated to 31 Mya, and the common ancestor of the four studied sections (β) is dated to 24 Mya (Marcussen et al., 2015). The phylogeny of Viola is based on the work of Marcussen et al. (2015) and that of Violaceae is based on the work of Wahlert et al. (2014).

be expressed, and the RNA sequence of some molecular species might be degraded during the transcriptomic assay.

Secondly, signature sequences are conserved between Viola and other Violaceae genera, i.e., Gloeospermum, Melicytus, and Pigea (Supplementary Figure 5). In particular, the sequence signature of the NTR [−9, −8] region is largely conserved in the Möbius and bracelet lineages: Y−9 -A−8 is found in the Möbius lineages of Viola and Melicytus, G−9 -A−8 in the bracelet lineages of Viola and Gloeospermum, and H−9 -S−8 in the bracelet lineages of Viola and Pigea. In addition, the sequence signature of the NTPP [−56, −38] is conserved in the bracelet lineage. That region shows high sequence similarity between some precursors, i.e., molecular species NS2 and HS3 from the genus Viola and precursors from Melicytus and Pigea, respectively. Interestingly, the protein sequence of Gpc3, a precursor belonging to the bracelet lineage, is identical in both Viola and Gloeospermum. These observations indicate that these sequence signatures have been conserved at least since the divergence of their most recent common ancestor some 50 million years ago (Marcussen, Wahlert et al., in prep.).

Thirdly, the cyclotide sequence diversity of Viola is expected to be dependent on the differentiation into sections. It is estimated that 67–88% of the Viola sections are allopolyploids (Marcussen et al., 2015) that combine genomes from different diploid lineages. Also, increase in ploidy has been shown experimentally to increase also sequence diversity of cyclotides by mutation, e.g., in Oldenlandia affinis (Seydel et al., 2007). However, the cyclotide sequence diversity could not only be very huge by the sequence variations within molecular species, but also be limited by sharing the same molecular species across different sections. The sequence diversity can be large by allowing conservative substitution within the same molecular species. Under the current classification (Marcussen et al., 2015), Viola comprises at least 16 sections and some 600 species, of which 10 sections (∼454 species, 76%) possess at least one CHAM genome, either alone or, as a result of allopolyploidization, in combination with other CHAM genomes or MELVIO genomes. These 10 sections include the four (∼380 species, 63%) inspected by the current study; and the current depth analysis reveals that most of cyclotide precursor sequences across different sections are likely grouped into the molecular species. The precursors identified in the current transcriptome and those reported previously all exhibit striking sequence similarity (Zhang et al., 2009, 2015b; Hellinger et al., 2015; Slazak et al., 2015).

#### Expression Profiles of Cyclotides and Their Structural Diversity

The expression profiles of cyclotides were analyzed at transcriptomic and peptidomic levels both. At the transcriptomic level, expression levels of precursor RNA were evaluated by their FPKM values. The expression of those sequences cyclotide sequences were then assayed at the peptidomic level by LC-MS.

At the transcriptomic level, a structurally diverse set of cyclotides were found, including archetypical, linear, hybrid and the novel types of cyclotides described above. In the bracelet lineage, the majority of cyclotide precursors were archetypical bracelets (73%, 125/171). In the Möbius lineage, the proportion of archetypical Möbius (33%, 26/78) was only second largest after the structural hybrids (35%, 28/78). However, judged by FPKM values, the expression levels of archetypical cyclotides were highest within each lineage. Most archetypical precursors have FPKM values higher than 150, i.e., 64% (=9/14) and 79% (=43/54) of the precursors in Möbius lineage and in bracelet lineage, respectively have values larger than the cutoff at 150 (Supplementary Table 4). On the other hand, precursors of novel cyclotides have low FKPM values, i.e., 53% (=7/13) and 100% (=13/13) of the precursors of novel cyclotides in the Möbius lineage and in the bracelet lineage, respectively. Moreover, in the Möbius lineage, 50% (=6/12) of the precursors encoding structural hybrids exceed the FPKM-cutoff value.

At the peptidomic level, numbers ranging from 19 to 44 cyclotides were detected in each of the Viola species. This number is similar to the number of precursor sequences found at transcriptomic level (23 to 37). Only few of the cyclotides (4– 26%) were detected at both levels (**Figure 6**; see Supplementary Table 5). Most of those cyclotides were archetypical Möbius (kB1 and kS) and archetypical bracelets (cyO2, cyO8, cyO13, mram8, and viba12). Only two of the novel cyclotides were found at both levels, and both were found eluting with very early retention time and at low protein levels of abundance. Cyclotides found at the peptidomic-level only were identified by matching their calculated MW, assuming that they could belong to any type, e.g., be hybrid, archetypical Möbius and bracelet. Discrepancy between expression at transcriptomic- and peptidomic levels has been reported previously (Burman et al., 2010; Hellinger et al., 2015). In those studies cyclotides expressed as peptides were all hybrids, or archetypical Möbius or bracelets. The failure of detection the novel cyclotides using LC-MS could be either from the large hydrophilicity of those novel cyclotides (Supplementary Figure 6) or from their low expression at the peptidomic level. Interestingly, kS (varv A) and kB1 were not found in the transcriptome from two Viola species (V. orientalis and V. albida var. takahashii), whereas these particular cyclotides were found at the protein level in all five Viola species. Such expression differences could be related to the molecular stability and regulation mechanism. However, it is clear that cyclotides accumulate after their expression in plant tissues.

### Major Sequence Changes Are Linked to Neofunctionality

In some cases, the differences between the precursor sequences of closely related molecular species were small (**Figure 7A**) while in other they were quite large (**Figures 7B–D**). In this context, a large difference is defined as an indel of more than five residues. With the exception of these indel regions, the precursor sequences are homologous, and the homologous sequence regions contain unique sequence traits found only in those molecular species. This implies that despite their large sequence differences, those molecular species are closely related. Such an outcome is analogous to that predicted by the theory of punctuated equilibrium (Eldredge and Gould, 1997), which describes evolutionary trends and speciation at the organism level. Another possible explanation for the observations is that mutations gradually accumulate in cyclotide precursors, leading to the loss or inactivation of intermediate sequences.

Some of the cyclotide domains exhibited unusual sequences, which fell into three main groups (**Figure 7E**). The first such group contained sequences whose cyclotide domains have an uneven number of cysteines (five or seven). Having uneven number of cysteines makes the formation of a cystine knot impossible, and can create the possibility of forming disulfide bridges not associated with the cystine knot, as well as the potential for dimerization. The second group consists of sequences that lack a residue required for AEP-mediated cyclisation, i.e., an N or D residue in loop 6. Abnormal precursor sequences of this sort are presumably the result of mutations in the genes encoding cyclic cyclotide precursors rather than being inherited from an ancestral linear cyclotide because they exhibit

FIGURE 6 | Cyclotide expression on peptidomic level. LC-MS chromatogram of the five Viola species used for transcriptome assays in current study (left panel). The abundance of cyclotides is linked to their signal intensity (Right panel). For each Viola species, the LC-MS chromatogram (Base Peak Ion) and mass distribution of the cyclotides are shown with retention times and abundance levels (High/Medium/Low). The symbols + and ◦ denote cyclotides that were not found in the transcriptome and were found in the transcriptome, respectively. The cyclotides (kB1, kS, cyO2, cyO8, and cyO12) were confirmed by comparing the retention time and isotopic mass derived from the LC-MS of V. tricolor and V. odorata as reference chromatogram.

strong sequence similarity with the former class of precursors in both loop 6 and the C-terminal tail. Conversely, typical linear cyclotide precursors have short loop 6 sequences, lack a C-terminal tail, and are quite abundant. The third group of abnormal precursor sequences included sequences with one or more unusually short or long loops. Only five such sequences were identified (among a total of 283), but the three types of abnormality often co-occurred in the same precursor. For example, valt1-FS3U contains one additional cysteine in loop 2, and lacks an N/D residue in loop 6. All of these precursors are only weakly expressed.

The cyclic backbone seems to have emerged by molecular speciation. Comparing sequences from linear and cyclic molecular species, i.e., the cyclic YY1 and linear YY2 in Möbius lineage, and the cyclic PS1 and HS4 to the linear PN1 in the bracelet lineage (**Figures 7F,G**), suggests that linearity is likely to be the primitive (ancestral) trait, rather than linear cyclotides having evolved from cyclic ancestors via the mutational loss of the N/D residue in loop 6 that is needed for cyclisation. Although the cyclotide precursor sequences are too variable at the nucleotide level to resolve the relationships between linear and bracelet (Huson and Bryant, 2006), linearity being ancestral is supported by three facts. First, linear cyclotides that have undergone such mutational losses normally have elongated loop 6 sequences with mutations in the C-terminal tail. Second, the prodomain sequences of these molecular species are well-defined (all of them share the same sequence signature). Third, the molecular species corresponding to these putative primitive linear cyclotides are found in all four studied sections of Viola and also in other genera (Mra13).

Some precursors sequences could also be grouped as rich in prolines (**Figures 7H,I**). The common structure-based criterion used to identify mature Möbius cyclotides is the presence of a cis Pro in loop 5. Classification based on this approach is not generally consistent with the evolutionary classification into lineages: many precursor sequences assigned to the Möbius lineage lack this Pro residue in loop 5 (P32) of the cyclotide domain. Examples include the molecular species YA1, HF1, and YS1-3, which exhibit appreciable sequence similarity with archetypical Möbius cyclotides in loops 2 and 3. These sequence intermediates, which have characteristics of both bracelet and Möbius cyclotides, have been classified as the hybrid subfamily, and it has been suggested that they are genetic chimeras of bracelets and Möbius cyclotides (Daly et al., 2006). Our analysis suggests that this is not the case, at least in the Violaceae. Instead, these structural hybrids and archetypical Möbius cyclotides seem to originate from the molecular species YY1 of the Möbius lineage, which in turn is closely related to the linear ancestor YY2.

A key question is whether the emergence of the cyclic backbone changed the properties of the membrane-interacting surfaces of the linear cyclotide sequences. To answer this question, we analyzed the structural traits (i.e., the exposed ratio and physicochemical properties) of the residues that differ between linear and cyclic cyclotides within Möbius and bracelet lineages. The lineages were considered separately, so the cyclic Möbius cyclotide kalata B1 was compared to the linear cyclotide violacin A from the Möbius lineage. Similarly, the cyclic bracelet cyclotide cyO2 was compared to the linear verec1-PN1 from the bracelet lineage. In both cases, the linear cyclotide appeared to be ancestral to its cyclic counterpart. Linear and cyclic cyclotides from the same lineage exhibit substantial sequence similarity, with extensive conservation of residues' physicochemical properties at key positions. This is illustrated by the example of the residues in loop 6 (**Figure 8B**). In the linear cyclotides, loop 6 is conformationally flexible, and the charged groups of the termini are largely exposed on the molecular surface (**Figures 8A,B**). The exposed ratios and physicochemical properties of the corresponding residues in the cyclic cyclotides are very similar to those in their linear counterparts, implying

FIGURE 9 | A potential link between precursor architecture and neofunctionality. (A) Suggested modes of duplication in precursor genes. In Viola, the precursor architecture differs between lineages. Precursors of bracelet lineages appear to contain only one modular domain (i.e., a concatenated NTR and a cyclotide domain), whereas precursors of the Möbius lineage may contain one to three such domains. The repeats in the Möbius lineage probably originated from internal gene duplication events. All precursors sharing the same architecture, in both the Möbius and bracelet lineages, probably originated from external gene duplication. (B) Sequence alignment between new types of cyclotides and previously identified archetypes. In the Möbius lineage, sequences vary between repeated cyclotide domains within the same precursor and cyclotide domains of different precursors. There are high levels of sequence conservation within repeated cyclotide domains, as exemplified by kB1 and kS within Vok1. However, the sequences of cyclotide domains in precursors containing multiple cyclotide domains in general are frequently very different, as exemplified by kalata S and viul F within prc-Viul F[C]-FS4U, and viba 19 and viba 21 within VbCP14. Of the precursor sequences shown in the figure, the molecular species YS7, which belongs to the Mobius lineage, has a similar sequence to archetypical bracelet cyclotides. For example, valta3-YS4 lacks a proline residue in loop 5 and has an elongated loop 3. Conversely, the molecular species FA1 contains lipophilic residues in loops 2 and 5, while loop 3 consists mainly of negatively charged residues. (C) Comparison of cyclotides' surfaces and electrostatic potentials. New types of cyclotides— viman-FA1, viman1-RS2, and valta3-YS4, and viul-F—are shown with their electrostatic potential surfaces. The distributions of their electrostatic potentials are generally dissimilar to those seen in typical bracelet and Mobius cyclotides. Residues with high exposed ratios are indicated by solid and dashed lines on the molecular surfaces. Residues highlighted with solid lines have dissimilar physicochemical properties to those found in the corresponding positions of archetypical cyclotides. Residues whose physicochemical properties match those of the archetypical cyclotides are indicated by dotted lines. Residues are numbered in accordance with their sequence positions in (B). Negatively charged surface regions are colored in red, positively charged regions in blue, and hydrophobic regions in green. The numbering of residues follow the consensus sequence of the cyclotide domain in Figure 3.

that the structural framework favoring membrane interaction emerged before the cyclic backbone.

Some of the cyclotide domains exhibit very high levels of sequence diversity between different molecular species. Some molecular species (e.g., YS1-4, RS2 and FA1) contain cyclotide domain sequences that are highly dissimilar to archetypical cyclotide sequences (**Figure 9**). We found that these molecular species occur across Viola sections, which means that they must have evolved before the differentiation of these sections some 15 million years ago. Although the biological function of all cyclotide forms is currently unclear, their maintenance within the genome suggests that they fulfill vital and divergent functions in their host plants.

It is likely that some cyclotides may have biological functions based on mechanisms other than membrane binding and disruption (Burman et al., 2011). For example, the cyclotides-FA1 might have different orientations to archetypical cyclotides, even when one compares them to others from the same lineage. They have lipophilic residues in loops 2 and 5 with large exposed areas on the molecular surface. The lipophilicity of this loop differs from that of archetypical cyclotides: the archetypal Möbius cyclotides do not have lipophilic residues in loop 2, and the archetypal bracelets do not have lipophilic residues in loop 5. In addition, loop 3 of the FA1 cyclotides is rich in negatively charged residues, making them unique among the molecular species found in the genus Viola and suggesting that their functions may be unrelated to membrane binding.

Contradicting previous hypotheses, it appears that the cyclotides of the Möbius lineage exhibit a somewhat greater diversity of sequence traits than those of the bracelet lineage. Such large sequence diversity of the Möbius lineage may be correlated with the combined occurrence of internal and external gene duplications in the cyclotide domain. Precursor architecture with such duplication of the cyclotide domain is illustrated in the **Figure 9A**. A large proportion of the precursors (72%; 56/78) belonging to the Möbius lineage contain multiple cyclotide domains that appear to have originated from internal duplication events. Conversely, no bracelet precursors having repeated

#### REFERENCES

cyclotide domains have yet been identified in the Violaceae. There are large sequence differences in certain precursors' cyclotide domains at positions close to the N-terminal prodomains, as can be seen by comparing the sequences of the FA1 and YS4 molecular species to the archetypal Möbius cyclotides varv A and kalata S (**Figures 9B,C**). Pronounced sequence differences also exist between repeated cyclotide domains in some individual precursors, as in the case of the precursor of Viul F. Together with the rest of our findings, this implies that combined internal and external duplications can synergistically produce large changes in cyclotides' sequences and structures, giving rise to new biological functions: neofunctionality.

# AUTHOR CONTRIBUTIONS

SP, TM, KJR, ID, AB, and UG designed experiments. SP carried out interpretation of precursor sequences, modeling and mass spectrometry. UG supervised the experiments. K-OY collected plant material for the transcriptome sequencing. EJ identified cyclotide genes from transcriptomic data. TM and SP performed the phylogenetic analyses. All authors contributed to the writing and revision of the manuscript.

### FUNDING

This research was supported by Swedish Research Council (#2012-5063). KJR is an ARC Future Fellow (FT130100890).

#### ACKNOWLEDGMENTS

We are thankful to Dr. Kyeong-Sik Cheon at Kangwon National University for valuable discussions.

#### SUPPLEMENTARY MATERIAL

The Supplementary Material for this article can be found online at: https://www.frontiersin.org/articles/10.3389/fpls.2017. 02058/full#supplementary-material

outside the box. J. Nat. Prod. 77, 724–736. doi: 10.1021/np40 1055j

with enhanced activity. Chembiochem 9, 1939–1945. doi: 10.1002/cbic.2008 00174

Gran, L. (1973b). Oxytocic principles of Oldenlandia affinis. Lloydia 36, 174–178.

of a novel precursor in Petunia (Solanaceae). J. Biol. Chem. 287, 27033–27046. doi: 10.1074/jbc.M112.370841

**Conflict of Interest Statement:** The authors declare that the research was conducted in the absence of any commercial or financial relationships that could be construed as a potential conflict of interest.

Copyright © 2017 Park, Yoo, Marcussen, Backlund, Jacobsson, Rosengren, Doo and Göransson. This is an open-access article distributed under the terms of the Creative Commons Attribution License (CC BY). The use, distribution or reproduction in other forums is permitted, provided the original author(s) or licensor are credited and that the original publication in this journal is cited, in accordance with accepted academic practice. No use, distribution or reproduction is permitted which does not comply with these terms.

# The Challenge to Translate OMICS Data to Whole Plant Physiology: The Context Matters

Marcelo N. do Amaral and Gustavo M. Souza\*

*Department of Botany, Institute of Biology, Federal University of Pelotas, Pelotas, Brazil*

Keywords: downward causation, emergent properties, hierarchical systems, plant signaling, systems biology

#### INTRODUCTION: SOME CHALLENGES

The exponential development of high-throughput technologies in the last decades, supporting and improving the OMICS science, has allowed uncovering successfully the complexity of the organizational network patterns in the cell's metabolism to the plant phenome, founding the science of system biology (Mochida and Shinozaki, 2011). Further, the huge data sets and growing computational power have stimulated scientists to glimpse about how plants respond to the environmental changes, and how such knowledge could engender new technologies, for instance, to increase crop yields (Edwards and Batley, 2004; Tardieu et al., 2017). Through these technologies, researchers are describing deeply the different hierarchical levels of plant organization, improving the possibility to predict the behavior of whole plant (phenome). Based on extensive analyses of gene expression (genome and transcriptome) and/or metabolic networks (metabolome), it has been possible to monitor and control cellular responses to genetic perturbations or environmental changes (Fukushima et al., 2009).

#### Edited by:

*Paulo Mazzafera, Universidade Estadual de Campinas, Brazil*

#### Reviewed by:

*Anthony Trewavas, University of Edinburgh, United Kingdom*

\*Correspondence: *Gustavo M. Souza<EMAIL_ADDRESS>
#### Specialty section:

*This article was submitted to Plant Systems and Synthetic Biology, a section of the journal Frontiers in Plant Science*

> Received: *13 September 2017* Accepted: *04 December 2017* Published: *13 December 2017*

#### Citation:

*do Amaral MN and Souza GM (2017) The Challenge to Translate OMICS Data to Whole Plant Physiology: The Context Matters. Front. Plant Sci. 8:2146. doi: 10.3389/fpls.2017.02146*

However, different constrains can make both the predictability and the controllability difficult from the bottom-up cause-effect approach that underpins the deterministic view of science based on an upward chain of causality (**Figure 1**) (Noble, 2008; Sheth and Thaker, 2014). The first "bottleneck" is how to integrate the massive datasets from molecular high-throughput technologies with the growing high-throughput information on the crop scale, i.e., plant phenomics (Fukushima et al., 2009; Tardieu et al., 2017), which is a typical problem of finding a proper (if it does exists indeed) scaling law (Souza et al., 2016).

The second constrain comes from the common assumption in biochemical models that the system sampled should be in a metabolic steady-state at a given moment, characterized by constant metabolite levels, and that the different metabolic pathways operate in isolation (Toubiana et al., 2013), which is an obvious oversimplification. For instance, at intracellular level, compartmentalization into organelles enables differences in metabolite concentrations acting as a barrier to passive diffusion between organelles and cytoplasm, creating a non-homogenous cellular metabolic space (Sweetlove and Fernie, 2013; de Souza et al., 2017). Moreover, each metabolic pathway, somehow, is integrated in a dynamical metabolic network (Toubiana et al., 2013), which is challenging for static networks mathematical models that often bypass the network modulation over time. For instance, the stomatal movement depends on a range of environmental and endogenous plant stimuli that affect the internal networks at multiple levels of cellular spatiotemporal organization, generating species-specific responses to combined external stimuli (Merilo et al., 2014). From the modeling of a single guard cell at steady-state, researches seek for elucidate how these interactions determine the phenotype of plants. However, due of this hierarchy of scales, the interpretation of a large set of data from OMICS tools becomes quite difficult, then it is necessary to develop new methods to allow investigations of dynamic aspects of large scale models (Medeiros et al., 2015).

The first "bottleneck" refers to the problem of emergent properties at the higher level organizations of the system that are not fully determined by the properties of the lower levels (Souza et al., 2016), for example, changes in the transcriptome or in the proteome do not always result in respective alterations in the metabolome (biochemical phenotype) that exhibits its own dynamics (Ryan and Robards, 2006). Additionally, there is the influence of downward causation processes (Noble, 2008) (**Figure 1**), when higher levels of organization affect the functioning of the lower levels. For instance, interlocked transcriptional/translational feedback loops are involved in the generation of circadian rhythm in plants, and the functional clock (higher level) controls a wide range of cellular processes such as gene expression (lower level; Fukushima et al., 2009).

The second constraint is related to the different sources of "uncertainty" operating in different levels of plant organization, blurring the predictability from lower levels. For instance, at cellular level, the sources of uncertainty emerge from the spontaneous thermodynamical noise of molecular activity constraining the flux of distributions in metabolic networks (Hoppe et al., 2007), interactions between genes that enable alternative routes for the same phenotype (Kohl et al., 2010), and epigenetic effects changing genes expression (Crisp et al., 2016). Moreover, at the level of whole plant integration, there are many types of long-distance signaling processes (chemical and electrical) over toping each other and engendering a highly complex informational network that feedback on the regulation of cells metabolism (Choi et al., 2016).

Further, it's worth to consider that, especially under stressful conditions, the different sources of external "noise" often affect the way that plants respond to environmental changes (Bertolli and Souza, 2013; Prasch and Sonnewald, 2015). Environmental fluctuations potentiate the accumulation of conserved cellular signals such as reactive oxygen species (ROS) and the modulation of intracellular Ca2+ (Chi et al., 2015; de Souza et al., 2017). Different types of ROS and oxidized molecules produced in different subcellular compartments, together with a spatial and temporal modulation of Ca2+ elicit different transcriptional responses and, in several cases, the expression of nuclear genes can be altered without altering the total concentration of the signaling molecule in the cell as a whole (Tuteja and Mahajan, 2007; Leister, 2012).

#### INTERCELLULAR AND WHOLE PLANT SIGNALING POTENTIATE THE CHALLENGES TO INTEGRATE OMICS INFORMATION ACROSS SCALES

In addition to intracellular complexity, the interaction between neighboring cells plays an important role in the responses to environmental conditions, and the plant metabolism organization as a whole. In plants, plasmodesmata connect the cytoplasts of adjacent cells across the cell wall, allowing intercellular transport and communication to adjacent cells within a tissue or organ, allowing exchange of small molecules, such as ions, sugars, and phytohormones, as well as larger molecules, including proteins, RNA, and viruses (Brunkard et al., 2013). The exchanges of different types of molecules among cells generate, within the same tissue or organ, different gradients of molecules and metabolites, increasing the complexity of physiological processes. In addition, the same type of signal often induces different calcium-dependent responses between two cells of the same type (Gilroy and Trewavas, 2001). Intercellular communication through plasmodesmata plays a crucial role in specifying the fate of cells, as well as in different responses of the same tissue to environmental conditions (Pyott and Molnar, 2015). An example of this can be seen in one of the mechanisms of root development regulation through the short-range cellto-cell movement of miR165/6 (Carlsbecker et al., 2010). The expression of mobile miR165/6 in the endoderm results in a morphogenic gradient, which extends into the xylem layers toward the root center. This generates an opposite PHABULOSA (PHB) expression gradient (regulated by miR165/6), which therefore has a higher concentration in internal xylem tissue. Thus, xylem tissue within the stele is defined, among other factors, by the expression of PHB, which is restricted to xylem and procambium by miR165/6, specifically expressed in the endoderm.

This intercellular communication through non-autonomous mobile signals adds a further challenge to OMICS approaches, because organs and plant tissues present a great heterogeneity in expression patterns and metabolite profiles, and this information can be lost upon tissue homogenization for downstream analyses.

Besides local communication mechanisms, plants developed long distance signaling processes that enable communication and systemic responses. This type of communication responds to a wide range of environmental stimuli in which the perceived signals are transmitted to the distal organs, inducing systemic responses. Several messengers have been proposed to mediate this systemic communication in plants such as ROS, electrical signals and Ca2+, appearing to be integrated, demonstrating a fast, complex, and finely tuned communication system (Gilroy et al., 2014).

The systemic responses increase system complexity (plant as a whole), and thus increase the uncertainties of the bottomup predictability models, since physiological changes in specific tissues may have non-local causes. For instance, a local application of high light results in the activation of a ROS wave, allowing an increase in stress tolerance accompanied by the accumulation of photorespiratory amino acids, including Glycine and Serine, in non-stimulated tissues (Suzuki et al., 2013). Examples also demonstrate that Ca2+ propagation increase in aerial parts of the plant induced by local root treatment with NaCl, showing propagation kinetics differences of leaf-to-leaf (Xiong et al., 2014). The miRNAs also act on systemic responses, such as miR399 that function as a signaling molecule between the aerial tissues and roots to regulate the uptake of inorganic phosphates (Pi) (Chiou et al., 2006).

Actually, the cell is the result of the properties that emerge from the complex interactions and spatial structures among the thousands of molecules and enzymes of which it is composed. In addition, the environmental context, whether from outside or inside the plant, contributes to shape the way that information is processed by each cell (Gilroy and Trewavas, 2001), and these properties expand at different scales within the plant (Souza et al., 2016). According to Vítolo et al. (2012), the observation of different scales of plant organization, under the same circumstances, can show remarkable differences in the responses to the same stimuli, allowing different interpretations if considering each scale as isolated (Stressed or not stressed? It's the question. . . ). For instance, when plants of soybean were subjected to drought, it was observed, in one hand, significant decreases in gas exchanges (reflecting reduction in plant growth) but, on the other hand, non-significant alterations in chlorophyll florescence neither in enzymatic antioxidant activity (Bertolli et al., 2014). Therefore, different scales of organization can show different homeostatic capacities when disturbed, supporting the hypothesis that there is not a privileged level of causation in biological systems (Noble, 2012).

#### CONCLUDING REMARKS

The Cartesian method proposed by René Descartes in the seventeenth century set that the first step in order to understand some natural phenomenon is to analyze it, i.e., to decompose the phenomenon in its constitutive parts and to understand them

#### REFERENCES