Chemotype classification and biomarker screening of male Eucommia ulmoides Oliv. flower core collections using UPLC-QTOF/MS-based non-targeted metabolomics

Introduction

Eucommia ulmoides Oliv. is a multifunctional tree species endemic to China. Due to geological processes during the quaternary glaciation, this monotypic genus remained primarily in Mid-eastern China (Call & Dilcher, 1997). After its introduction to other regions of China, E. ulmoides was planted in 350,000 hm² of land across 28 provincial administrative districts (Du, Hu & Yu, 2013). The bark of E. ulmoides has been used as a traditional herbal medicine for almost 2,000 years (Feng & Liang, 1996). Currently, the tree is fully utilized as a part of the rubber modification, medicine, health food, pollution-free breeding and landscaping industries (Du, Hu & Yu, 2013, 2015).

Eucommia ulmoides is strictly dioecious and its male capitulum inflorescence consists of 5–11 clustered achlamydeous flowers. Several bractlets are attached to its tiny pedicel, and each male flower is composed of 8–12 stamens. In the Yellow River basin and other parts of northern China, the male flower blooming period is between mid-March and early April and lasts approximately 20 days. This is 8–15 days earlier than the bolooming period of the female flower (Liu, 2010; Liang, 2007). The natural compounds that have been extracted from male flowers, including geniposide, geniposidic acid, aucubin, chlorogenic acid, and flavonoids, have important anti-inflammatory, bacteriostatic, antihypertensive, antihyperlipidemic, hepatoprotective, choleretic, and cardiovascular protective qualities (Dong et al., 2005; Zhang et al., 2007). In 2014, the National Health Commission of China added male E. ulmoides flowers to their novel food list (National Health Commission, 2014).

The research and innovation team of the National Germplasm Resources Bank of Major Famous Tree Species in the North, located in Yuanyang, Henan Province, has conserved nearly 2,000 E. ulmoides germplasm resources, including more than 750 male collections, since 2009. Previous studies have evaluated the diversity of the morphological traits, amino acid content, and main active components of 193 male collections (Du et al., 2016, 2017a, 2017b), and 33 male core collections were screened by analyzing the phenotypic and SSR molecular marker data (Li et al., 2018a, 2018b). However, these studies did not explore how metabolic composition is used to classify male flower collections or determine the existence of reasonable chemical indicators when distinguishing traits, which are important to the quality formation and commercialization processes.

The metabolic composition and chemotype of different germplasm resources should be accurately evaluated to effectively expand genetic base, breeding innovation, and direct cultivation. Non-targeted metabolomics is an unbiased approach that comprehensively excavates the metabolic fingerprints of biological samples, despite having limited linear range and poor repeatability (Fan et al., 2017). Because of the particularly high performances in both mass spectrometry (MS) and tandem MS (MS/MS) modes, the combination of ultra-performance liquid chromatography (UPLC) separation and quadrupole time-of-flight mass spectrometry (QTOF-MS) detection has shown absolute superiority in qualitative and quantitative compound analysis (Taamalli et al., 2014; Wolfender et al., 2008; Da Silva et al., 2017). Coupled with multivariate statistical methods, UPLC-QTOF-MS is an optimized compound characterization approach that has been used to analyze ginseng, Ficus deltoidea, Garcinia oblongifolia, and many traditional Chinese medicinal herbs (Zhu et al., 2019; Afzan et al., 2019; Ning et al., 2013; Yu et al., 2018).

In this study, we profiled the metabolite constituents of 22 fresh male E. ulmoides core collections flowers grown under similar environmental conditions and cultivation practices. We processed data using multivariate statistical methods, including principal component analysis (PCA), hierarchical cluster analysis (HCA), and orthogonal partial least squares discriminant analysis (OPLS-DA), to determine the possible phytochemical phenotypes and metabolomic markers across the collections. Our study provides a comprehensive view of E. ulmoides male flower metabolomes and a basis for the future utilization of germplasm evaluation, nutritional properties, and biological analysis.

Materials and Methods

Results

Metabolite annotation

The total ion chromatograms (TICs) acquired in ESI+ and ESI− modes are shown in Figs. S2 and S3. We obtained 11,306 features (including 9,050 high-quality features) and 8,616 features (including 7,487 high-quality features) from the extracted ion chromatograms acquired in ESI+ and ESI− modes, respectively. Since a single isomer feature may match several metabolites from the first-order mass spectra, we calculated the mass feature number corresponding to the candidate metabolites. From the first-order mass spectra, we annotated 6,092 metabolites observed in ESI+ mode, and 3,554 metabolites observed in ESI− mode. From the secondary mass spectra, we annotated 451 and 325 metabolites in ESI+ and ESI− modes, respectively (Table 1). After KEGG pathway enrichment analysis, we found that the metabolites acquired in ESI+ and ESI− modes were involved in 69 and 78 metabolic pathways, respectively. The number of metabolites involved in the top 20 KEGG pathways is displayed in Fig. 1.

Table 1:

Statistics of the annotated metabolites.

Mode	Number of all features	Number of high quality features	Number of MS1 candidate metabolites	Number of MS1 annotated metabolites	Number of MS2 annotated metabolites
ESI+	11,306	9,050	9,686	6,092	451
ESI−	8,616	7,487	6,030	3,554	325

DOI: 10.7717/peerj.9786/table-1

Chemotype classification

Since the metabolites annotated from the MS2 spectra were more accurate and unique than those annotated from the MS1 spectra, we used the MS2 annotation data for chemotype classification of the E. ulmoides samples. We followed two matrices: one using 132 samples × 451 metabolites for the ESI+ mode data, and the other using 132 samples × 325 metabolites for the ESI− mode data. The PCA score scatter plot and HCA dendrogram of the ESI+ mode data matrix can be found in Fig. 2. We used HCA to divide the 22 collections into four groups according to their chemotypes. The seven collections AG-1, BJ-5, MC, SNJ-2, SQ-1, SQ-2, and XS comprised chemotype I; the three collections BJ-3, CL, and HZ comprised chemotype II; the six collections BJ-1, BJ-2, BZ, LC, LY, and SNJ-1 comprised chemotype III; and the six collections AG-2, BJ-4, JA, ZY-1, ZY-2, and ZZ comprised chemotype IV. The total variance of the data using the PCA model was 81.55%, with 26.06% accounted for by principal component 1 (PC1) and 15.82% by principal component 2 (PC2). The two-dimensional PCA score plot provided a visualization of the separation and similarity across the chemical groups. Chemotype I was mainly placed on the left side along PC1, chemotype IV was placed in the center, and chemotype III was placed on the right side. Meanwhile, along PC2, chemotype II was located in the upper half and the other three chemotypes were located in the lower half. Most biological duplications from each collection clustered together, and the boundaries of the classified chemotypes were evident in the PCA score plot, suggesting that the analytical procedures were consistent and repeatable. We did not observe separation based on the geographical origin of the germplasm collections from the analytical results, only that the collections from Shangqiu (SQ-1 and SQ-2) and Zunyi (ZY-1 and ZY-2) clustered in one chemotype.

We applied the PCA model to the ESI− mode data matrix which accounted for 90.40% of the total variance, of which PC1 accounted for 24.36% and PC2 accounted for 16.34% (Fig. 3). Similar to the analysis of the data acquired in ESI+ mode, HCA divided the 22 collections into four chemotypes. Chemotype I was comprised of the five collections AG-1, AG-2, JA, SNJ-2, and ZZ; chemotype II was comprised of the five collections BJ-2, BJ-4, LC, LY, and SNJ-1; chemotype III was comprised of the five collections BJ-5, MC, SQ-1, SQ-2, and XS; and chemotype IV was comprised of the seven collections BJ-1, BJ-3, B2, CL, HZ, ZY-1, and ZY-2. The classification results were slightly inconsistent with those from the ESI+ mode data, especially for chemotype I and chemotype IV. The ESI+ mode identified a higher number of metabolites than ESI− mode,and most metabolites including amino acids, short-chain carboxylic acids and flavonoids were detected in both modes. The chemotype classification boundaries at multiple dimensions in the PCA plots were unclear from the ESI− mode data (Figs. S4 and S5), so we ultimately selected and adopted the classification results originating from the ESI+ mode to establish the chemotype for the collections.

To determine whether a certain classified chemotype of the male Eucommia core collections was unique enough to be distinguished from other types, we constructed four supervised OPLS-DA models to provide further insights into the separation (Fig. 4). To assess each OPLS-DA model, we calculated R²Y (cum) and Q²Y (cum) coupled with two indicators referring to the intercept of R² and Q² from the permutation test (Table 2). We used the generally accepted threshold values (>0.5) of R²Y (cum) and Q²Y (cum) to judge whether the OPLS-DA model was stable and effective (Triba et al., 2015), and clearly found that the constructed models achieved these standards, confirming that each classified chemotype had a distinct metabolite composition. However, the permutation test showed that the OPLS-DA model between chemotype IV and the other chemotypes had slightly higher R² and Q² intercept values, suggesting the possibility of over-fitting in that model.

Table 2:

Main assessment index of each OPLS-DA model.

Model	Orthogonal component	Predictive component	R2Y (cum)	Q2Y (cum)	Intercept of R2	Intercept of Q2
Chemotype I against the other types	1	2	0.904	0.878	0.231	−0.361
Chemotype II against the other types	1	4	0.908	0.818	0.451	−0.619
Chemotype III against the other types	1	3	0.914	0.861	0.326	−0.424
Chemotype IV against the other types	1	5	0.920	0.808	0.513	−0.673

DOI: 10.7717/peerj.9786/table-2

Screening chemical markers

We determined the most relevant chemical markers for discriminating each chemotype based on the VIP value, p-value, and FC value in the OPLS-DA model. In total, we screened 29 chemical markers, 21 of which corresponded to chemotype I, six corresponded to chemotype II, 19 corresponded to chemotype III, and five corresponded to chemotype IV (Fig. 5; Table S3). The chemical markers were flavonoids, amino acids and their derivatives, cinnamic acids and their derivatives, lipids, alkaloids, organooxygen compounds, and organoheterocyclic compounds. The overlapping correlation of the markers is illustrated in Fig. 6. For chemotype I, the exclusive markers were cinnamic acid, L-phenylalanine, glycyl-L-leucine, and jasmonic acid. For chemotype III, the exclusive markers were S-adenosyl-L-methionine and astragalin. Chemotype II and chemotype IV each contained only one exclusive chemical marker, carnosic acid and isorhamnetin, respectively. The chemical marker variations and the peak intensities in the four chemotypes are shown as box-and-whisker plots (Fig. S6). To re-examine the validity of the markers in denoting different chemotypes, we built a data matrix consisting of the 29 markers and 22 collections. We used the same clustering method to draw a heat map (Fig. 7) that showed the samples were well-classified and in agreement with the chemotypes established above. However, the HZ collections were misclassified.

Discussion

Eucommia ulmoides, a monotypic genus of Eucommiaceae, has long been considered to have relatively limited intraspecific variations compared to other species belonging to multi-type genera. In the 1980s, there was high demand for E. ulmoides bark, significantly increasing its market value. Germplasm resources, especially aged E. ulmoides trees, were nearly devastated, leading many researchers to believe that the species’ wild resources had been exhausted in China (Deng, Li & Huang, 2006; Wang et al., 2006). The majority of existing trees were planted from homogeneous provenances and propagated by grafting, which was previously believed to limit genetic diversity (Wang et al., 2006). The sustainable cultivation and utilization of the present resources and effective breeding is one of the major concerns of E. ulmoides researchers. The first steps should be the implementation and evaluation of the species variations at phenotypic, genotypic, and chemotypic levels. A previous study conducted high-resolution metabolic phenotyping using metabolomic approaches to assess a broad range of possible metabolites and enable informed selection in plants (Fernie & Schauer, 2009). Our work was the first attempt at disclosing the metabolome of male E. ulmoides flowers and the chemical variations of the core collections from different geographical origins under similar cultivation measures. These results could provide a foundation for further research on the medical use of these male flowers and their metabolic network and regulation.

Due to LC-MS instruments’ versatility and sensitivity, they can detect hundreds to thousands of features from a single sample, and, collectively, tens of millions of features for large-scale projects in non-targeted metabolomics studies (Rocca-Serra et al., 2016). A significant portion of the detected metabolomes were highly redundant due to the presence of isotopes, adducts, and in-source fragments (Kachman et al., 2019). Robust annotation and identification of metabolites in non-targeted studies remains a major challenge in LC-MS-based non-targeted metabolomics, which might lead to false discoveries and incorrect interpretations (Salek et al., 2013; Scheubert et al., 2017). Other issues, including matching multiple identities with a single feature and multiple disagreeing annotations, often arise when matching the detected features to candidate chemical structures based on their mass and predicting fragmentation patterns based on chemical libraries (Alden et al., 2017). However, it is not realistic for most individual laboratories to establish a comprehensive spectral library that encompasses all possible chemical standards (Broeckling et al., 2016). Although most modern validation strategies for non-targeted approaches are not well defined (Naz et al., 2014), the metabolites annotated in our study cover most well-known chemical constituent categories, provide the metabolomic profiles for male E. ulmoides flowers, and present a potential holistic approach to improve germplasm resource evaluation and understand its physiological mechanisms.

The unsupervised clustering method HCA coupled with PCA reduced the dimensionality of the multivariate data and created new principal components using linear combinations to determine similarities and differences across the samples. These methods have been applied in previous plant metabolomic analyses (Gad et al., 2013; Cevallos-Cevallos et al., 2009). Furthermore, the combined application of PCA and OPLS-DA to spectral datasets yielded valuable insights into general spectral trends and group predictive spectral features (Worley & Powers, 2016). The core collections are small sets of accessions chosen to represent the genetic spectrum of the whole collections and cover high variation levels, representing up to 70% of the genetic diversity of a species. We used multivariate methods to objectively classify the four chemotypes of the male E. ulmoides core collections, referencing the main chemical variation forms in the germplasm resource. It should be noted that the biological replicates from the six collections SNJ-1, BZ, BJ-3, BJ-1, HZ, and LY were not clustered within one chemotype. In addition to errors resulting from using the employed clustering method, we assumed that some random sampling errors, injection problems, and other instrument noises or drifts were the main factors affecting the clear classification of the six collections’ biological replicates. However, according to our Spearman’s correlation analysis based on the peak intensity data, the correlation coefficient was around or above 90%, and we observed a higher correlation in replicates of the same collection than in different collections (Fig. S7). This indicated that there were significant similarities in the collection duplications and they could therefore be identified as belonging to the same chemotype. Due to their different physicochemical properties (e.g., polarity), the two data sets detected using the ESI+ and ESI− modes were also different to some extent, leading to slightly inconsistent HCA chemotype classification. Because of the great dimension and complexity of a non-targeted metabolomics data matrix, it is hard to obtain a strictly precise classification using existing bioinformatics tools. In this study, we verified the effectiveness of the classification results using the OPLA-DA models. However, both classified chemotype sets were inconsistent with geographical origin patterns. We inferred that the chemical variation may have originated from the inherent genetic variation of the germplasm and intrapopulation (not interpopulation) variation. Similar results at the genetic level were found in previous studies (Wang et al., 2006; Wu, 2014).

Other studies have explored the main categories of the nutritive components of male E. ulmoides flowers, including iridoids, phenylpropanoids, flavonoids, amino acid derivatives, and phenolic glycosides (Yan et al., 2018; Ding et al., 2014). In our current study, we screened out 29 chemical biomarkers belonging to seven metabolite classes, the majority of which were flavonoids and amino acid derivatives. Comparison and variation analyses on the 17 amino acids in E. ulmoides male flowers of 193 germplasm resources showed that the content variation degrees, evaluated by coefficient of variation (CV) values, were relatively higher in proline (Pro), cysteine (Cys), lysine (Lys), and methionine (Met) (Du et al., 2016). Another study evaluated the content variations of the eight major active components in the same resources, and found that geniposide, geniposidic acid, and astragaline had the highest variation degrees (CV values above 80%); followed by isoquercitrin, chlorogenic acid, and aucubin (CV values around 50%); and total flavonoids and total amino acids had the lowest variation degrees (CV values below 25%) (Du et al., 2017b). We assumed that the chemical components with the highest degrees of content variation in the germplasm resources would be potential biomarkers for chemotype differences, and biomarkers such as chlorogenic acid, methionine, and S-adenosyl-L-methionine complied with this assumption. However, the most likely biomarkers geniposide, geniposidic acid, and astragaline did not comply with our results. One possible reason for these inconsistencies may be that the potential biomarkers’ features were not fully annotated because the coverage of standard spectra in the reference database was incomplete. Additionally, these were two studies conducted using different experimental conditions, compound extraction procedures, and measuring instruments. Future studies should validate the identification and structure of the screened biomarkers using multiple optimized in silico tools, and confirm their results using classical targeted analytical approaches such as the NMR technique to get more reliable and reproducible results.

Conclusion

We determined the metabolite profiles of male E. ulmoides flowers from different geographical origins using a developed UPLC-QTOF/MS fingerprinting method. We found abundant compositional variations and performed phenotyping at the metabolic level across several male E. ulmoides core collections, providing a basis for understanding the formation mechanism of economically important traits. A total of 451 and 325 metabolites were identified in ESI+ and ESI− modes, respectively, from 22 male flower core collections. We coupled the multivariate statistical method PCA with OPLS-DA models to provide a new technique for classifying E. ulmoides germplasm resources that is efficient in revealing the major metabolites that contribute to chemotype discrimination. Four chemotypes could be determined from the MS2-annotated metabolites that were detected in ESI+ mode. We proposed that compositional variations across the chemotypes were due to genetic composition variations instead of geographical distributions. We screened a total of 29 robust chemical markers that enabled the differentiation of four established chemotypes, most of which were flavonoids and amino acid derivatives. The newly proposed chemical markers, along with existing monograph and molecular markers, could be utilized for more precise germplasm resource identification. It should be noted that our results are preliminary and require further verification using a larger set of samples.

Supplemental Information