Literature DB >> 30646498

Assessing Parameter Suitability for the Strength Evaluation of Intramolecular Resonance Assisted Hydrogen Bonding in o-Carbonyl Hydroquinones.

Maximiliano Martínez-Cifuentes1, Matías Monroy-Cárdenas2, Juan Pablo Millas-Vargas3, Boris E Weiss-López4, Ramiro Araya-Maturana5,6.   

Abstract

Intramolecular hydrogen bond (IMHB) interactions have attracted considerable attention due to their central role in molecular structure, chemical reactivity, and interactions of biologically active molecules. Precise correlations of the strength of IMHB's with experimental parameters are a key goal in order to model compounds for drug discovery. In this work, we carry out an experimental (NMR) and theoretical (DFT) study of the IMHB in a series of structurally similar o-carbonyl hydroquinones. Geometrical parameters, as well as Natural Bond Orbital (NBO) and Quantum Theory of Atoms in Molecules (QTAIM) parameters for IMHB were compared with experimental NMR data. Three DFT functionals were employed to calculated theoretical parameters: B3LYP, M06-2X, and ωB97XD. O…H distance is the most suitable geometrical parameter to distinguish among similar IMHBs. Second order stabilization energies ΔEij(2) from NBO analysis and hydrogen bond energy (EHB) obtained from QTAIM analysis also properly distinguishes the order in strength of the studied IMHB. ΔEij(2) from NBO give values for the IMHB below 30 kcal/mol, while EHB from QTAIM analysis give values above 30 kcal/mol. In all cases, the calculated parameters using ωB97XD give the best correlations with experimental ¹H-NMR chemical shifts for the IMHB, with R² values around 0.89. Although the results show that these parameters correctly reflect the strength of the IMHB, when the weakest one is removed from the analysis, arguing experimental considerations, correlations improve significantly to values around 0.95 for R².

Entities:  

Keywords:  DFT; NBO; QTAIM; hydrogen bond; hydroquinone; polyphenol

Mesh:

Substances:

Year:  2019        PMID: 30646498      PMCID: PMC6359028          DOI: 10.3390/molecules24020280

Source DB:  PubMed          Journal:  Molecules        ISSN: 1420-3049            Impact factor:   4.411


1. Introduction

Hydrogen bonding (HB) is generally described as a stabilizing electrostatic interaction with a partly covalent character. Some of the most direct experimental evidence of the HB formation is the deshielding of the hydrogen atom involved in the interaction, which is easily observed by NMR spectroscopy [1,2]. It is well known that intramolecular hydrogen bonding (IMBH) has a key role in determining the structure and properties of biologically-active molecules, and their study in model compounds is relevant in drug design [3,4,5]. For example, IMBH formation is determinant for passive membrane permeability of cyclic peptides [6] and in lipophilicity [7]. Any intramolecular hydrogen bonding IMHB is governed by the equilibrium between closed (IMHB) and open (no IMHB) forms. Consequently, the implementation of IMHB considerations in drug discovery programs requires experimental methodologies to assess the position of these equilibria, which are dependent on the strength of each IMBH [8]. The 1H-NMR spectroscopy is the most common tool to measure the propensity of compounds to form IMHBs, nevertheless, up to now the throughput of NMR, and other traditional methods, is limited in obtaining structure–activity relationships (SAR) and driving rapid cycles of compound optimization [8]. However, current computational chemistry allows one to assess, theoretically, the strength of IMBH through parameters such as the distance hydrogen-acceptor, the distance donor-acceptor, and the hydrogen bond angle [9]. Among the theoretical estimated parameters for the HB interaction, those from Natural Bond Orbital (NBO) and Quantum Theory of Atoms in Molecules (QTAIM) analysis [9,10,11,12,13] are the most popular. In this field, DFT methodologies have proven to be valuable tools to calculate molecular properties with an appreciably low computational time demand compared with ab initio methodologies. Previously we have studied the properties of some ortho-carbonyl hydroquinones with capability to cross membranes and to reach the inner mitochondrial membrane. These compounds act either as oxidative phosphorylation (OXPHOS) uncouplers, complex I inhibitors, or exert both activities; as a consequence they exhibit antitumor [14,15,16] or antiplatelet effects [17], or both. Interestingly, small structural changes on the hydroquinone scaffold determine their biological activities [15]. All these molecules exhibit strong resonance assisted intramolecular hydrogen bonding, assessed by 1H-NMR spectrometry [18]. These compounds act either as weak acids into the mitochondrial membrane, or by losing a hydrogen atom (or alternatively losing an electron followed by deprotonation) to afford the corresponding semiquinone radical, which can interfere in cell signaling. The formation of IMHB plays a key role [19,20] on both properties. IMHBs have shown appreciable effects on the antioxidant and electrochemical properties of hydroquinones and related phenols [21,22]. On the other hand, the o-carbonyl hydroquinone motif is present on natural products with different biological activities, such as the anticancer drugs doxorubicin and daunorubicin [23], acylnaphtohydroquinones [24], shikonin [25], and in peyssonol A [26]. The precise assessment of the strength of the IMBH interaction has many difficulties. Determination of the H-bond energy is possible only for the intermolecular hydrogen bond: the energy of the intermolecular interaction is usually equal to the energy required to separate the two interacting molecules [1]. However, there is no clear definition of the hydrogen bond energy for the IMHB because the interacting entities cannot be fully separated [27]. Moreover, it is relatively easy to roughly calculate the relative strength of a series of IMBHs when they take place in molecules with large structural differences; however, to assess the relative strength of IMBHs in very similar molecules, establish correlations with experimental parameters, and use these data to make correlations with biological activities, requires an adequate choice of theoretical methodologies. The aim of this work is to study the IMHBs experimentally and theoretically in a series of structurally similar o-carbonyl hydroquinones. Calculated NBO and QTAIM parameters of the IMHBs obtained by different DFT functionals (B3LYP, M06-2X, and ωB97XD) will be correlated with experimental values obtained by NMR.

2. Results

The studied compounds were synthesized according to Scheme 1 [28,29]. Starting from hydroquinone (1) or 2,3-dimethylhydroquinone (2), the corresponding acylhydroquinones (3–5) were obtained by Fries rearrangement [30]. These compounds were oxidized with silver (I) oxide to obtain the corresponding quinones (6–8), which were not isolated. These quinones yield benzofurans (12–17) by reaction with the enamines (9–11). The last step involves an acid rearrangement of these benzofurans, yielding the corresponding bicyclic hydroquinones (18, 20, 22–24). Monomethyl ethers of hydroquinone 19 and 21 were obtained by methylation of hydroquinones 18 and 20, respectively; they were synthesized to assess the effect of blocking the hydroxyl group not involved in IMHB.
Scheme 1

Synthesis of the bicyclic hydroquinones studied in this work.

Three hybrid functional were employed in this work: B3LYP, composed of Becke three parameters hybrid functional B3 and the non-local correlation functional of Lee, Yang, and Parr (LYP) [34,35,36]. The second corresponds to a hybrid meta-generalized gradient-approximation (hybrid meta-GGAs) functional, M06-2x, which have been described as one the of the best Truhlar group functionals for non-covalent interactions [37]. The last functional corresponds to the long-range corrected functional ωB97XD, which is designed to avoid rapid die-off of the non-coulomb part of the exchange functionals, which affects the modelling of long distance processes [38]. In this first part, we study main geometrical parameters from the IMHBs and their correlation with 1H-NMR chemical shifts. The results are presented in Table 1. Calculations were carried out at DFT level using the B3LYP, M06-2X, and ωB97XD functionals with the 6-311++g(d,p) basis set. Among the geometrical parameters, only the O15…H16 distances displayed enough differentiation to be used as comparison parameters for correlation with the chemical shifts of chelated protons. Compound 21 displayed the largest chemical shift and a more de-shielded proton, indicating the strongest IMHB, followed in strength by 20. The later indicates that methyl substituents on C1 and C2 increase the IMHB strength. Additionally, the methylation at O12 and also a methyl group at C9 increase the strength of the IHB in 21 and 22, respectively. The three functionals give the shortest O15…H16 distance for 21, indicating the strongest IHB, which is in agreement with the largest downfield chemical shift. The smallest chemical shift for 18 indicates the weakest IHB in the series. Nevertheless none of the three functional gives the largest O15…H16 distance for 18. These discrepancies can reflect the presence of intermolecular interactions competing with the IMHB, contributing in a non-negligible degree to the 1H-NMR spectra at the probe temperature, especially between pairs of molecules where one phenolic hydroxyl is replaced by a methoxy group. The same considerations can be applied to correlations with the next parameters. Global correlation among the chemical shift and O15…H16 distance for all molecules were carried out (Figure 1). Results indicate that the ωB97XD functional gives the best correlation with an R2 = 0.89, followed by M06-2X with a R2 = 0.85 and B3LYP, and the worst correlation with R2 = 0.80.
Table 1

Experimental chemical shifts (δ) in ppm of protons involved in hydrogen bonding, and main geometric calculated parameters. Distances in angström, (Å), and angles in degree (°).

Moleculeδ ppmDistance (O15H16)Distance (O11-H16)Angle (O15H16-O11)
B3LYPM062xωB97XDB3LYPM062xωB97XDB3LYPM062xωB97XD
18 12.5401.6571.6941.6740.990.980.982148146147
19 12.6901.6571.6951.6750.990.980.981148146147
20 13.2001.6391.6621.6450.9920.9830.984149147148
21 13.4401.6201.6581.6400.9940.9840.985149148149
22 12.9301.6601.6911.6660.990.980.982148146147
23 12.7401.6611.7011.680.990.980.982148146147
24 12.8001.6541.691.6720.990.9810.982148146147
Figure 1

Correlations between H16…O15 distance and 1H-NMR chemical shifts (δ) for the hydrogen bond.

In order to estimate the energy of IMHB, the open-close method was employed, which defines the H-bond energy as the difference between the open and close conformers [39,40]. Calculations were carried out with the functional ωB97XD, which presents the best correlation between experimental chemical shift and O16…H15 distance. Table 2 shows that energy for IMHB in the studied hydroquinones are around 15 kcal/mol. Hydroquinone 20 presents the highest value, 16.37 kcal/mol, while 19 presents the lowest IMHB energy, 14.22 kcal/mol. These energy differences are large enough to discard the contribution of the open form to the experimental NMR shieldings. Figure 2 shows the correlation between Hydrogen bond energy (EHB) and 1H-NMR chemical shifts (δ), which presents a R2 = 0.71, a value lower than the correlation with the O15…H16 distance (R2 = 0.89), calculated with the same functional.
Table 2

Hydrogen bond energy (EHB), calculated as the difference between open and close conformers. Energies in kcal/mol.

MoleculeEHB
18 14.34
19 14.22
20 16.37
21 15.71
22 14.44
23 14.36
24 14.52
Figure 2

Correlation between Hydrogen bond energy (EHB) and 1H-NMR chemical shifts (δ).

The results from NBO analysis are presented in Table 3. Second order stabilization energies ΔEij(2) between donor (φi) and acceptor orbital (φj) allow one to study orbital interactions responsible for the IMHB. For all studied compounds, the HB donor atom is always an oxygen (two lone-pairs) and the acceptor is the antibonding sigma orbital of an O-H groups (σ*OH). As we found in the geometrical analysis, strong IMHB presence in 21 (the largest chemical shift) is in agreement with the largest ΔEij(2) obtained from the calculations with the three functionals. Nevertheless, ΔEij(2) did not predict well the weakest IMHB, corresponding to 18, and the lowest stabilization energy was not predicted by any of the three calculations using the studied functionals. Global correlation among the chemical shift and ΔEij(2) for all HQs were carried out (Figure 3). As in the geometrical analysis, correlation results indicate that the best correlation is obtained with the ωB97XD functional, with R2 = 0.89. However, in this case it was followed by B3LYP with R2 = 0.87 and then M06-2X, the one that gave the worst correlation with R2 = 0.85.
Table 3

Second order stabilization energies (ΔEij(2)) between donor LP O15 and acceptor σ* O11-H16 orbitals involved in the IMHB. Energies in kcal/mol.

MoleculeΔEij(2) LP O15 → σ* O11-H16
B3LYPM06-2XωB97XD
18 27.1423.1225.11
19 27.0222.9524.87
20 28.9926.1327.99
21 31.1526.628.62
22 27.6523.2325.68
23 26.6722.5424.59
24 27.4723.4425.27
Figure 3

Correlations between second order stabilization energies ΔEij(2) LP O15 → σ* O11-H16 and 1H-NMR chemical shifts (δ) for IMHB.

The QTAIM analysis is presented in Table 4. It shows negative values for ∇2ρ and H for all studied compounds, indicating the covalent nature of their IHBs. |V|/G ratio is also a sensitive parameter for evaluating the covalency [41,42]. V can be interpreted as the pressure exerted on the electrons by the HB system, while G can be interpreted as the pressure exerted (as a reaction) by the electrons on the HB system, at the critical point. For negative values of ∇ρ, as occurred in all our cases, values of |V|/G > 2 were indicative of covalent interaction [2]. Hydrogen bond energy (EHB) was obtained from the potential energy density (V) [13]. All values were slightly above 30 kcal/mol, while NBO analysis gave all energies slightly below 30 kcal/mol for the same interaction. These values were approximately twice those calculated by the open-close method, which indicates that NBO and AIM overestimate the interaction energy for the IMHB. EHB values calculated by different methods have shown to be generally inconsistent [43]. As we found in both geometrical and NBO analysis, the strongest IHB is obtained for compound 21 and is well represented by QTAIM analysis for the three functionals, giving the highest EHB for 21. On the other hand, similarly to NBO analysis, here QTAIM EHB also did not adequately reflect the weakest IMHB present in 18. Global correlations between QTAIM EHB and chemical shifts are presented in Figure 4. A better performance is shown by ωB97XD and B3LYP with a R2 = 0.88. As in the above cases, M06-2X gave the worst correlation with R2 = 0.85.
Table 4

Atoms-in-molecule parameters for IHB of 18–24. Electron density at the critical point ρBCP (a.u), its Laplacian ∇2ρ (a.u.), electron kinetic energy density G (a.u.), potential energy density V (a.u.), total electron energy density H (a.u.), and hydrogen bond energy EHB (kcal/mol).

FunctionalMoleculeρBCP2ρGV|V|/GHEHB
B3LYP 18 0.0537−0.03770.0459−0.10132.2070−0.147231.91
19 0.0537−0.03780.0460−0.10152.2065−0.147531.96
20 0.0560−0.03830.0478−0.10512.1987−0.152933.12
21 0.0587−0.03920.0501−0.11002.1956−0.160134.65
22 0.0546−0.03830.0469−0.10332.2025−0.150232.54
23 0.0533−0.03740.0455−0.10032.2044−0.145731.59
24 0.0541−0.03780.0463−0.10202.2030−0.148332.14
M06-2X 18 0.0454−0.04080.0444−0.09902.2297−0.143431.19
19 0.0454−0.04090.0444−0.09902.2297−0.143431.19
20 0.0490−0.04270.0479−0.10662.2255−0.154533.57
21 0.0495−0.04280.0483−0.10742.2236−0.155733.83
22 0.0458−0.04120.0449−0.10012.2294−0.145031.53
23 0.0446−0.04020.0435−0.09702.2299−0.140530.57
24 0.0459−0.04100.0448−0.09992.2299−0.144831.48
ωB97XD 18 0.0510−0.03800.0447−0.09892.2125−0.143531.14
19 0.0508−0.03810.0446−0.09872.2130−0.143431.10
20 0.0546−0.03940.0479−0.10572.2067−0.153633.30
21 0.0553−0.03960.0485−0.10692.2041−0.155333.66
22 0.0518−0.03860.0457−0.10102.2100−0.146631.81
23 0.0502−0.03760.0439−0.09732.2164−0.141230.64
24 0.0512−0.03810.0449−0.09932.2116−0.144231.29
Figure 4

Correlations between hydrogen bond energy EHB (kcal/mol) and experimental chemical shifts δ (ppm) of proton involved in IMHB.

The results shown in Figure 3 corroborates that correlations with the weakest IMHB of this series are not well represented by the computational calculations employed. In order to test what is the impact of the weakest IMHB on the correlations studied previously, compound 18 was eliminated from the data and correlations were obtained again. The results showed a significant improvement, without change in the relative position among different functionals and parameters studied. Figure 5 shows correlations between 1H-NMR chemical shifts and O11-H16…O15 distance for the IHBs. R2 for B3LYP changed from 0.80 to 0.88, from 0.85 to 0.91 for M06-2x, and from 0.89 to 0.95 for ωB97XD.
Figure 5

Correlations between O-H…O distance and 1H-NMR chemical shifts (δ) for the IMHB, excluding the weakest value.

Similarly, correlations between 1H-NMR chemical shifts and NBO energy for the IHBs are presented in Figure 6. R2 changed from 0.88 to 0.95 for B3LYP, from 0.85 to 0.92 for M06-2X, and from 0.89 to 0.96 for ωB97XD.
Figure 6

Correlations between second order stabilization energies ΔEij(2) LP O15 → σ* O11-H16 and experimental chemical shifts (δ) of protons involved in IMHB, excluding the weakest value.

Finally, correlations between 1H-NMR chemical shifts and QTAIM energy for the IHBs, presented in Figure 7, showed a change in R2 from 0.88 to 0.94 for B3LYP, from 0.85 to 0.91 to M06-2X, and from 0.88 to 0.94 for ωB97XD.
Figure 7

Correlations between hydrogen bond energy EHB (kcal/mol) and experimental chemical shifts δ (ppm) of protons involved in IMHB, excluding the weakest value.

3. Materials and Methods

3.1. Synthetic Procedures

Melting points are uncorrected. All NMR spectra were acquired using a Bruker AVANCE DRX 300 spectrometer (Bruker BioSpin GmbH, Rheinstetten, Germany) operating at 300.13 MHz (1H) or 75.47 MHz (13C). Measurements were carried out at a probe temperature of 300 K.

3.1.1. 1-(5-Hydroxy-3,3,6,7-tetramethyl-2-morpholino-2,3-dihydrobenzofuran-4-yl)ethan-1-one (14)

To a solution of 1-(2,5-dihydroxy-3,4-dimethylphenyl)ethan-1-one (400 mg, 2.2 mmol) in 20 mL of dichloromethane, Ag2O (1.28 g) was added and stirred with magnetic stirring for one hour. Then, the filtered solution was added over a solution of 4-(2-methylprop-1-en-1-yl)morpholine (314 mg, 2.2 mmol), in CH2Cl2 (10 mL), yielding 637 mg of the furan (14) (2.0 mmol, 90% yield). 1H-NMR δ(CDCl3): 1.37 (s, 3H, 3-CH3); 1.45 (s, 3H, 3-CH3); 2.14 (s, 3H, Ar-CH3); 2.19 (s, 3H, Ar-CH3); 2.33–2.60 (m, 4H, 2× CH2-O); 2.61 (s, 3H, CH3CO); 3.70 (m, 4H, 2× CH2-N); 4.68 (s, 1H, O-CH-N); 7.09 (s, 1H, 5-OH). 13C-NMR δ(CDCl3): 11.82; 12.51; 19.80; 31.17; 32.81; 45.07; 49.28; 66.75; 108.78; 120.94; 122.53; 123.74; 128.69; 147.01; 150.86; 205.43. M.P.: 140.1–142.2 °C. IR (KBr) ν/cm−1: 842.77; 1684.89; 2836.66; 2909.00; 2949.52; 2981.35; 3357.56. HRMS calculated: 319.1784; found: 319.1787.

3.1.2. 5,8-Dihydroxy-4,4,6,7-tetramethylnaphthalen-1(4H)-one (21)

A solution of furan 14 (451 mg, 1.4 mmol) in ethanol (15 mL) was refluxed for 2 h, then the solution was poured in a mixture ice and water, and the solid precipitate was filtered and dried yielding hydroquinone 21 (256 mg, 78% yield). 1H-NMR δ(CDCl3): 1.61 (s, 6H, 2 × 4-CH3); 2.23 (s, 3H, Ar-CH3); 2.26 (s, 3H, Ar-CH3); 4.55 (s, 1H, 5-OH); 6.23 (d, 1H, J = 10 Hz, 2-H), 6.82 (d, 1H, J = 10 Hz, 3-H), 13.20 (s, 1H, 8-OH). 13C-NMR δ(CDCl3): 11.56; 13.04; 25.17; 38.03; 112.68; 123.43; 124.04; 131.51; 131.64; 143.52; 155.12; 160.73; 191.11. P.F.: 200.5–201.6 °C. IR (KBr) ν/cm−1: 1658.84; 2963; 2998.71; 3406.75. HRMS calculated: 232.1099; found: 232.1113.

3.1.3. 8-Hydroxy-5-methoxy-4,4,6,7-tetramethylnaphthalen-1(4H)-one (22)

To a solution of hydroquinone 21 (100 mg, 0.43 mmol) in acetone (20 mL), powdered K2CO3 (85 mg, 0.6 mmol) and dimethyl sulfate (0.1 mL, 0.6 mmol) were added and stirred with magnetic stirring at reflux for 3 h. To the filtered mixture, a solution of KOH (5% in methanol) was added and kept at room temperature overnight. Then, it was neutralized with HCl 10% and extracted with CH2Cl2 5 × 20 mL. The organic phase was washed with water and dried with anhydrous sodium sulfate, then was filtered and the solvent was eliminated at vacuum, by purification in flash column chromatography with hexane/ethyl acetate (5:1) as eluent, the product 22 was obtained (79 mg, 74% yield) 1H-NMR δ(CDCl3): 1.57 (s, 6H, 2 × 4-CH3); 2.20 (s, 3H, Ar-CH3); 2.31 (s, 3H, Ar-CH3); 3.72 (s, 3H, 8-OCH3); 6.24 (d, 1H, J = 10 Hz, 2-H); 6.78 (d, 1H, J = 10 Hz, 3-H); 13.44 (s, 1H, 5-OH). 13C-NMR δ(CDCl3): 11.39; 14.77; 27.09; 38.15; 61.49; 112.36; 124.13; 124.84; 137.88; 139.43; 149.29; 157.10; 160.76; 190.92. P.F.: 82.1–84.3 °C. IR (KBr) ν/cm−1: 1655.95; 2848.23; 2920.58; 2972.67; 3001.61; 3042.12; 3441.48. HRMS Calculated: 246.1256; found: 246.1265.

3.2. Computational Calculations

The calculations were carried out using the Gaussian 09 [44] program package (Revision a.01; Gaussian, Inc.: Wallingford, CT, USA). No symmetry constraints were imposed on the optimizations, which were performed at the DFT B3LYP/6-311++G (d,p) level. No imaginary vibrational frequencies were found at the optimized geometries, indicating that they are true minima of the potential energy surface. NBO analysis was performed using the NBOPro 6.0 [45] program package (NBO 6.0; University of Wisconsin: Madison, WI, USA). QTAIM analysis was performed using the AIM2000 [46] program package (Büro für innovative software, Bielefeld, Germany).

4. Conclusions

In this work, we found that O…H distance is the most suitable geometrical parameter for distinguishing among IMHBs of similar strengths. Second order stabilization energies ΔEij(2) from NBO analysis and hydrogen bond energy (EHB) obtained from QTAIM analysis were also able to properly distinguish among IMHBs in the studied molecule series. In all cases, ΔEij(2) from NBO give values for the IMHB lower than EHB obtained from QTAIM analysis. ΔEij(2) are slightly below 30 kcal/mol, while EHB from QTAIM gave values somewhat above 30 kcal/mol, which are approximately twice of those calculated by the open-close method. In all cases, parameters calculated using ωB97XD gave the best correlations with experimental 1H-NMR chemical shifts for the IMHBs, with values of R2 around 0.89. Despite this, the results showed that the employed parameters correctly described the strength of the IMHB, when the weakest one was removed from the analysis, correlations improved significantly to values around 0.95 for R2. These results can help to select the most suitable parameters for modeling molecules with IMHBs of similar strengths and achieving an adequate distinction between them through computational calculations.
Hydroquinone R1 R2 R3 R4
18 [31]HHMeH
19a [31]HHMeMe
20 MeHMeH
21 b MeHMeMe
22 [32]HMeMeH
23 [32]HH-CH2-(CH2)-CH2-H
24 [15,33]HHEtH

Note: a Obtained by methylation of 18, b obtained by methylation of 20.

  2 in total

1.  Introduction to "Intramolecular Hydrogen Bonding 2018".

Authors:  Goar Sánchez
Journal:  Molecules       Date:  2019-08-07       Impact factor: 4.411

2.  Glycerol as Alternative Co-Solvent for Water Extraction of Polyphenols from Carménère Pomace: Hot Pressurized Liquid Extraction and Computational Chemistry Calculations.

Authors:  Nils Leander Huamán-Castilla; María Salomé Mariotti-Celis; Maximiliano Martínez-Cifuentes; José Ricardo Pérez-Correa
Journal:  Biomolecules       Date:  2020-03-20
  2 in total

北京卡尤迪生物科技股份有限公司 © 2022-2023.