Effect of Different Processing Methods on the Chemical Constituents of Scrophulariae Radix as Revealed by 2D NMR-Based Metabolomics

Scrophulariae Radix (SR) is one of the oldest and most frequently used Chinese herbs for oriental medicine in China. Before clinical use, the SR should be processed using different methods after harvest, such as steaming, “sweating”, and traditional fire-drying. In order to investigate the difference in chemical constituents using different processing methods, the two-dimensional (2D) 1H-13C heteronuclear single quantum correlation (1H-13C HSQC)-based metabolomics approach was applied to extensively characterize the difference in the chemical components in the extracts of SR processed using different processing methods. In total, 20 compounds were identified as potential chemical markers that changed significantly with different steaming durations. Seven compounds can be used as potential chemical markers to differentiate processing by sweating, hot-air drying, and steaming for 4 h. These findings could elucidate the change of chemical constituents of the processed SR and provide a guide for the processing. In addition, our protocol may represent a general approach to characterizing chemical compounds of traditional Chinese medicine (TCM) and therefore might be considered as a promising approach to exploring the scientific basis of traditional processing of TCM.

In TCM, processing is a common practice, aiming at enhancing efficacy and reducing toxicity in clinical application [15][16][17]. According to the Chinese Pharmacopoeia (2020), SR needs to be processed by "sweating" [14]. Sweating processing is one of the common methods in TCM raw material processing. The detailed procedure for SR is as follows. The raw material is baked with a low fire to semi-dry and then piled up to sweat to make the internal moisture spillover. The procedure is repeated several times until the material is dry. The sweating process makes the SR become soft and change color to increase fragrance and reduce irritation [18]. There are two other important methods for SR processing recorded in ancient classical TCM processing books such as Lei Gong's Moxibustion Theory and the Compendium of Materia Medica, using steaming and hot-air drying methods [19]. Although there are two articles describing the chemical constituents in SR after different processing methods [8,10], they only investigated the change of a limited number of known components in SR. However, for traditional Chinese medicine with complex chemical components, the information obtained from a few known chemical components is not comprehensive enough to clarify the impact of traditional processing on SR. Therefore, it is essential to develop a method to comprehensively characterize compounds in SR and its processed products.
Untargeted metabolomics is a technology allowing an unbiased and hypothesis-free assessment of the metabolome [20]. Nuclear magnetic resonance (NMR) is a powerful platform in untargeted metabolomics due to its ability to analyze essentially all types of known and unknown natural products [21]. Recently, two-dimensional (2D) NMR has attracted attention because it can provide more detailed metabolite fingerprints [22,23]. The 2D NMR experiment commonly used in metabolomics is 2D 1 H-13 C heteronuclear single quantum coherence (HSQC) [24]. The large 13 C chemical shift dispersion in the 2D 1 H-13 C HSQC spectra allowed us to analyze a wide range of metabolites without prior separation [25], and the information provided by specific atomic correlations can assist in the identification and assignment of known metabolites as well as illuminate the NMR information of unknown metabolites [26]. Based on the rather inclusive and publicly available databases (HMDB, BMRB), multiple well-resolved 1 H-13 C cross-peaks in 2D 1 H-13 C HSQC spectra have been frequently used for metabolite identification. A few works have studied the characteristic metabolic changes of the crust from dry-aged beef and metabolites of six different laboratory Escherichia coli strains by using this method [27,28]. However, few studies have used 1 H-13 C HSQC spectra in untargeted metabolomics analysis of TCM.
In our current study, 2D 1 H-13 C HSQC-based metabolomics was used to extensively characterize the chemical components in SR processed by freeze-drying, steaming for different durations (S), hot-air drying (HD), and "sweating" (SW). Firstly, the extraction conditions of SR metabolomics samples were optimized to fully extract metabolites and save data acquisition time. Secondly, we developed two in-house TopSpin AU programs to automate the batch processing of 1H NMR and an in-house NMRPipe script to automate the processing of all 2D 1H-13C HSQC spectra in-batch. Finally, the differential features selected were analyzed using NPid and preliminarily assigned by comparing with the customized NMR database of SR obtained through literature data mining and public spectral databases (HMDB, BMRB). This study provides a useful and quick strategy to comprehensively analyze chemical compounds of SR with the 2D 1 H-13 C HSQC features and therefore might be considered as a promising approach for exploring the scientific basis in traditional processing of SR.

Samples Collection and Processing
A total of 99 samples were obtained by processing fresh roots of SR with four different methods (FD, S, SW, HD). The detailed information on sample processing is shown in Table S1.

Selecting the Optimum Extraction Solvent
Dried samples of SR (250 mg) were exactly weighed and extracted in parallel with five different solvents, namely 20%, 40%, 60%, 80%, 100% v/v methanol/water solutions. The 1D 1 H NMR spectra of the five extracts were obtained. Segmented absolute integrals of three regions were calculated except for the signals from residual water and protonated solvents, including the aromatic region (10.00-5.60 ppm), the polar (carbohydrates) region (5.60-3.33 ppm), and the aliphatic region (3.30-0.05 ppm) ( Table S2). The absolute segmental integrals of the extracts with five different extraction solvents were compared to evaluate the extraction efficiency of the SR metabolites. The 60% methanol/water solutions could extract more aromatic compounds and carbohydrates, whereas 80% methanol/water solutions could extract more aliphatic compounds. However, 60% methanol/water solutions could extract more global metabolites. Finally, 60% methanol/water solutions were selected as the optimum extraction solvents.

Optimization of the Solid/Liquid Ratio of the Extraction
The effect of solid/liquid ratio on the extraction yield was investigated with a SR samples vs. different volumes of 60% methanol/water solutions at 20, 60, 100, 140, 180, 200 mL/g. With the increasing volume of solvent, the total integrals of 1 H NMR spectra increased gradually and reached the plateau when the solvent to solid ratio was 140 mL/g ( Figure S1), indicating that the volume of solvents used was enough to completely extract the SR metabolites. In practice, considering the possible variation of the quantity of metabolites among different SR samples, the optimum solvent to solid ratio was appropriately relaxed to 160 mL/g to ensure complete extraction of SR metabolites.

Optimum Sample Amount Used per NMR Sample
The use of excessive SR per NMR sample may lead to oversaturation of the extracts in the NMR deuterated solvents, which may bring about non-linear response or even insensitivity of the intensity of NMR signals to the real quantity of metabolites in the original SR sample. In the aromatic region, the range of SR sample amount of 100 to 500 mg was located in the linear regime. In the sugar region, the range of SR sample amount of 100 to 300 mg was located in the linear regime. In the aliphatic region, the range of SR sample amount of 100 to 500 mg was located in the linear regime ( Figure S2). The integral of 1D 1 H NMR spectra showed that the optimum amount of SR used per NMR sample was 250 mg, and the concentration of the NMR samples made the acquisition time of 1 H-13 C HSQC spectra reasonably short (1 h 52 min).

Multivariate Statistical Analysis
In order to study the effects of different steaming durations on the metabolites of SR, PCA and OPLS-DA were performed on FD and steamed samples (S01, S02, S04, S08, S12, S24, S48, S72). Figure 1 showed that the difference between steamed and FD samples became more and more significant with the extension of steaming time. That is, the longer the steaming time, the more significant the change of chemical composition of SR. Comparison of FD and S01 in the PCA score plot of 2D 1 H-13 C HSQC spectra showed no significant separation between FD and S01 ( Figure 1A), whereas FD and S02 were more dispersed ( Figure 1B). Comparison of FD and S04-S72 ( Figure 1C-H) showed a gradually separation of the two groups. In the supervised OPLS-DA score plots (Figure 2A-H), clear discrimination was observed between FD and steamed samples (S01, S02, S04, S08, S12, S24, S48, S72). All OPLS-DA models are valid with R 2 (cum) and Q 2 (cum) values exceeding 0.8 (Table S3). Variables with VIP > 1.1 and |p(corr)| > 0.9 were selected as significantly differential features and highlighted in magenta in the volcano plot (Figure 2a-h). Comparison between FD and steamed samples (S01, S02, S04, S08, S12, S24, S48, S72), with total 215 features (10, 10, 20, 28, 29, 61, 123, 150 for samples at each steaming time, respectively, calculate the intersection of eight groups), showed significantly differential features (Table 1). Univariate statistical analysis was performed to further confirm that the selected features were statistically significant with p value < 0.01 using Student's t-test.

Change of Chemical Constituents of SR Steamed for Different Durations
Out of the 215 differential features, 131 features were assigned to 20 compounds ( Table 1). The identified differential metabolites were further confirmed by overlaying and comparing the 2D 1 H-13 C HSQC spectra of the mixture with the standards of the candidate metabolites. The boxplots of the intensities of 20 differential metabolites in SR steamed for different durations are shown in Figure 3. Taking FD samples as a reference, the relative contents of 15 differential metabolites, namely harpagoside, cinnamic acid, harpagide, aucubin, 6-O-methyl-catalpol, amino acids (leucine, threonine, tyrosine, methionine, valine, arginine, isoleucine), and oligosaccharides (sucrose, raffinose, and stachyose), decreased significantly (p value < 0.01) at varying steaming duration, while the relative contents of three differential monosaccharides, namely β-glucose, fructose, and glucuronic acid, increased with the steaming time. The relative contents of glutamine and 4-aminobutyric acid increased at the beginning (from 1 to 24 h and from 1 to 2 h, respectively) and then gradually decreased with steaming time.
With the extension of steaming time, the decrease in oligosaccharides and glycosides and increase in monosaccharides indicated possible hydrolysis of the oligosaccharides and glycosides during the steaming. We speculate that the decrease in the amino acids with the extension of steaming could be related to thermally induced Maillard reaction [29], which might be largely responsible for the browning sesames of the processed SR [30]. The boxplots of the intensities of the other 84 unidentified differential spectral features are shown in Figure S3. Further phytochemical study on the steamed SR might be helpful to determine their chemical identities.
The majority of chemical constituents of SR showed significant change after being steamed for 4 h, while the active components iridoids did not decrease significantly. Therefore, S04 samples were selected to represent samples processed by steaming and compared with SW and HD samples. In order to study the effect of different processing methods on chemical constituents of SR, PCA and PLS-DA were performed on the samples processed by SW, HD, and S04. PCA gave a model with five principal components (PCs), which cumulatively accounted for 74.7% of the total variance. In the PCA score plot, HD and S04 overlapped with each other, suggesting that these were relatively similar in their compositional patterns in 2D 1 H-13 C HSQC spectra. Compared with HD and S04, SW showed significant metabolic differences in the PCA model ( Figure 4A). PLS-DA generated a two-dimensional model with goodness of fit R 2 (cum) of 0.934 and predictive power Q 2 (cum) of 0.901. Further, the permutation test checked the validity and the degree of overfit for the PLS models. As shown in Figure S4, the intercept (R 2 = 0.209, Q 2 = −0.23) in the permutation plot is a measure of the overfit, which suggests that the models are reliable and not overfitted. The LV1 and LV2 components of the model explained 42.2% and 10.7% of the total variance, respectively. The differently processed samples were clearly divided into three groups in the score plot of PLS-DA ( Figure 4B). A total of 81 features were selected as preliminary differential variables based on their VIP (>1.0) and |p(corr) | (>0.8) values as shown in the volcano plot ( Figure 4C,D). Univariate statistical analysis was performed to further confirm that the selected features (Table 2) were statistically significant with p value < 0.01 according to ANOVA.

Identification of Differential Metabolites of SW, HD, and S04
The effects of different processing methods on chemical constituents of SR can be characterized using 81 differential features, of which 49 features were assigned to seven compounds ( Table 2). The boxplots of the intensities of these seven differential metabolites are shown in Figure 5. Compared with SR processed by HD and S04, the relative content of six differential metabolites, harpagoside, cinnamic acid, 6-O-methyl-catalpol, aucubin, methionine, and tyrosine, in SR processed by SW was generally lower, while the relative content of α/β-glucose was higher in SR processed by SW. Figure 5. Boxplots of the relative intensities of seven identified differential compounds in samples processed using different methods, sweating (SW, claret), hot-air drying (HD, peach), and steaming for 4 h (S04, black).
The boxplots of the intensities of the other 32 unidentified differential features in 1 H-13 C HSQC spectrum are shown in Figure S5. Among these 32 unidentified features, 16 of them were characteristic of SR processed by SW, while three of them did not show in SR processed by SW. Further phytochemical study to elucidate the chemical identities of these unidentified features is necessary and of significance to deepening our understanding of the effects of different processing methods on the chemical constituents of SR. Table 2. Detailed information and assignments of the differential features (F) among differently processed SR with VIP > 1.0 and |p(corr)| > 0.8 by PLS-DA and p value < 0.01 by univariate ANOVA.

Chemicals and Reagents
Methanol-d 4

SR Raw Material Samples Processing
Samples went through processing with different methods: vacuum freeze-drying, steaming (S) for different durations, sweating (SW), and hot-air drying (HD). The FD samples (n = 9) were first frozen in an ultra-low temperature refrigerator at −80 • C and then placed in a freeze dryer with a cold-trap temperature of −80 ºC to dry to constant weight. Since the steaming time was extremely important to the processing, the steamed samples were divided into different subgroups according to the steaming time for 1 (S01, n = 9), 2 (S02, n = 9), 4 (S04, n = 9), 8 (S08, n = 9), 12 (S12, n = 9), 24 (S24, n = 9), 48 (S48, n = 9), and 72 h (S72, n = 9), and all steamed samples were subsequently dried to constant weight at 55 • C in the oven. The SW samples (n = 9) were processed as follows: dried for 24 h at 55 • C in the oven, placed into a bag for 48h as the first sweating then removed from the bag, dried for 24 h at 55 • C in the oven, placed into a bag for 48 h as the second sweating, and finally removed from the bag and dried to constant weight at 55 • C in the oven. The HD samples (n = 9) were processed with oven drying at 55 • C. All processed samples were then ground into rough powder and kept in a refrigerator at 4 • C for subsequent extraction and NMR sample preparation.

NMR Sample Preparation and Data Acquisition
NMR samples were prepared as follows: 250 mg of ground SR was weighed and extracted with 40 mL 60% methanol/water solutions (v/v) in a 100 mL Erlenmeyer flask by ultrasonication at room temperature for 1 h. After centrifugation (10 min, 3077 g), the obtained supernatant was collected and evaporated using a rotary vacuum evaporator at 45 • C. Then, 600 µL of deuterated solvents (360 µL CD 3 OD + 240 µL D 2 O + 0.36 mg TSP) was added to redissolve the dried extracts of SR. The mixtures were then ultrasonicated for 3 min at room temperature and vortexed for 1 min. Finally, after centrifugation (10 min, 13,225 g), 500 µL of supernatant was transferred into a 5 mm NMR tube for NMR measurement.
For each sample, 1D 1 H and 2D 1 H-13 C HSQC spectra were recorded at 298 K on a Bruker Avance III spectrometer (operating at 700.2034 MHz for 1 H and 176.0795 MHz for 13 C, respectively) equipped with a 5 mm QCI-F CryoProbe. The 1D 1 H NMR spectra were acquired with a transmitter frequency of 4.798 ppm, relaxation delay of 2 s, 64 scans, 8 dummy scans, and a spectral width of 16 ppm. The 2D 1 H-13 C HSQC spectra were acquired with a 1.5 s relaxation delay; 8 transients; 512 increments; spectral widths of 12.98 and 150.00 ppm in F2 and F1, respectively; and offset 4.798 and 80 ppm for the 1 H and 13 C dimension, respectively.

NMR Spectrum Processing and Data Preprocessing
Two in-house TopSpin AU programs were developed to automate the batch processing of 1D 1 H NMR spectra for all samples. One was for automatic Fourier transformation (EFP), zeroth order phase adjustment (APK0), baseline correction (ABS), and chemical shift calibration (SREF) of all 1D 1 H NMR spectra. Considering the variation of the overall biomass among samples, the other AU program was developed to automate the segmented integration of all 1D 1 H NMR spectra for the regions of 12.0-4.90, 4.78-3.38, and 3.30-0.05 ppm with the solvent peaks of water and residual proton of methanol-d4 excluded, and the integral of each 1 H NMR spectrum was used as the normalization factor for the processing of the corresponding 2D 1 H-13 C HSQC spectra.
An in-house NMRPipe script was used to automate the processing of all 2D 1 H-13 C HSQC spectra in-batch. To obtain unbiased and comparatively meaningful absolute intensities, the scaling parameter (NC) was read and applied to 2D 1 H-13 C HSQC raw data besides normalization by the integral of the corresponding 1 H NMR spectrum.
Peaks were identified and quantified in all 2D 1 H-13 C HSQC spectra using a Javabased program, Newton, which employed an algorithm called fast maximum likelihood reconstruction (FMLR) for 2D 1 H-13 C HSQC spectral deconvolution to extract peaks and to obtain accurate signal quantitation [31]. Newton supports batch analysis of a series of related datasets, and all processed 2D 1 H-13 C HSQC spectra can be imported into Newton together. An appropriate threshold was set, and a prototype peak was selected as model peak for peak fitting. The batched spectral peaks obtained from Newton analysis were then preliminarily matched, forming peak clusters across samples using an in-house developed R program (manuscript in preparation). After comparing and uniquely designating the ambiguous peaks to a specific peak cluster and merging multiple redundant peaks from the same sample in the same peak cluster, the obtained peak clusters were further filtered to generate the valid peak clusters (spectral features). The valid peak clusters were then refilled. Finally, the regularized matrix of peak intensity (M) with the different features and samples as columns and rows, respectively, was obtained for the following chemometric analysis. Excluding the signals from the reference TSP and the residual protonated solvents, the final data matrix of peak intensity has 474 features and 99 measurements (that is, number of samples) for each feature; that is, the size of the matrix is 99 × 474.

Chemometric Analysis
All data were mean centered and scaled to unit variance (UV) as a preprocessing step then analyzed with principal component analysis (PCA) and orthogonal partial least squares discriminant analysis (OPLS-DA) or partial least squares discriminant analysis (PLS-DA) based on the matched and aligned spectral peaks obtained by Newton analysis, using SIMCA-P (version 14.1, Umetrics, Umea, Sweden). Seven-fold internal cross-validation was also performed to verify these models [32][33][34]. The quality of the OPLS-DA and PLS-DA models was revealed by the parameters of R 2 and Q 2 [32,35]. R 2 indicates the goodness of fit, and Q 2 represents the predictive power of the model. The values of R 2 and Q 2 close to 1.0 indicate an excellent model [32,35]. The values of variable importance for projection (VIP) and p(corr) were used to evaluate the contribution or importance and the relevance of individual X-variables to the OPLS-DA or PLS-DA model [27,36]. The variables with a VIP value higher than 1 and high p(corr) correlation value are considered as potential discriminant variables associated with the classification [36,37].

Metabolite Identification and Verification
We previously developed an automatic approach, NPid to rapid identification of known natural products in a mixture of extracts based on a 2D 1 H-13 C HSQC spectrum [38]. For NPid analysis, a customized NMR database of SR was constructed by using the chemical shifts and structural information of natural products in SR retrieved from published phytochemical research. In 2D 1 H-13 C HSQC spectra, each chemical constituent could be characterized by multiple well-resolved 1 H-13 C cross-peaks.
The differential features selected were analyzed using NPid and preliminarily assigned by comparing with the customized NMR database of SR obtained through literature data mining and public spectral databases (HMDB, BMRB). Each metabolite was identified and represented by multiple resolved 1 H-13 C cross-peaks, and each 1 H-13 C HSQC cross-peak contained paired chemical shifts of 1 H and 13 C. When standard compounds were available, the preliminarily identified metabolites were further verified by comparing the 2D 1 H-13 C HSQC spectra of the mixture with the standards of the candidate metabolites.

Conclusions
In this study, a 2D 1 H-13 C HSQC-based untargeted metabolomics approach was successfully applied to extensively characterize the difference in chemical components of SR processed by steaming for different durations and using different processing methods. These findings could elucidate the change of chemical constituents of the processed SR and provide a guide for the processing. Comprehensive characterization of the effect of the processing methods on the chemical constituents of TCM is a premise for follow-up research, such as exploring the therapeutic basis and establishing quality control standards and disclosing the underlying secret of the effectiveness and mechanisms of TCM processing.
Supplementary Materials: The following supporting information can be downloaded at: https: //www.mdpi.com/article/10.3390/molecules27154687/s1, Figure S1: Results of optimum ratio of sample weight to solvent volume; Figure S2: Results of optimum amount of Scrophulariae Radix used for per NMR sample of 60% methanol: water solutions; Figure S3: Boxplot of intensities of unidentified differential features at different steaming time; Figure S4: The permutation plot for the PLS model; Figure S5: Boxplot of 32 unidentified differential variables in Scrophulariae Radix processed by different methods. Table S1: Detailed information on Scrophularia ningpoensis Hemsl sample processing; Table S2: Absolute segmental integrals information of 1 H NMR spectra of SR extracts from five solvents; Table S3: Performance parameters R 2 (cum) and Q 2 (cum) of OPLS-DA models.  Data Availability Statement: Compiled data is reported in the tables above. The raw data files are available from authors upon request.

Conflicts of Interest:
The authors declare that they are not aware of competing financial interests or personal relationships that could have appeared to influence the work reported in this paper. The manuscript described has not been submitted elsewhere for publication, in whole or in part, and all the authors listed have approved the manuscript that is enclosed.