The reference standards of macrolide antibiotics (erythromycin, roxithromycin, clarithromycin, azithromycin, impurity J, midecamycin, josamycin, and acetylspiramycin) and azithromycin impurities (impurities F, J, I, L, Q, R, and S) were obtained from the National Institutes for Food and Drug Control (NIFDC). The purity of all compounds was >95%.

Wild-type AB strain zebrafish (Danio rerio) and transgenic Tg (fabp10: dsRed) zebrafish in this study were reared in a standard laboratory environment with a 14-h light/10-h dark cycle at 28 ± 1°C (Chen et al., 2019). As previously described (Zhang et al., 2020), 3 days post-fertilization (3 dpf), zebrafish larvae were exposed to compound solutions at different concentrations for 72 h. For visual assessment of liver toxicity in zebrafish, the tested zebrafish larvae were imaged by using a fluorescence microscope to observe the liver phenotype of zebrafish larvae, 20 larvae per group.

As previously described (Han et al., 2019; Zhang et al., 2020), the ADMET parameters were predicted by using the pKCSM online protocols (http://biosig.unimelb.edu.au/pkcsm/prediction). The topological polar surface area (TPSA) and the octanol–water partition coefficient (logP) were calculated by Molinspiration tools (http://www.molinspiration.com/. molinspiration.com/), the two chemical descriptors were primary factors of absorption. In this study, the absorption level depends on factors including water solubility and the apparent permeability coefficient indicated by colon cancer cell line (Caco-2). The distribution level depends on factors including the volume of distribution (VDss) and the blood–brain barrier (logBB). The metabolism of drugs was predicted based on two main isoforms of cytochrome P450 (CYP2D6 and CYP3A4 substrates). The excretion level indicated by total clearance was measured by the proportionality constant Cltot. The toxicity was predicted based on the hepatotoxicity. These pharmacokinetic parameters were calculated and examined for compliance with standard ranges based on human data.

The docking procedure in the present investigation was performed by CDOCKER algorithm on Discovery Studio (DS) 2020 software. For ligand preparation, the 3D structures of macrolide compounds were extracted from the PubChem Compound database (https://www.ncbi.nlm.nih.gov/pccompound/). CHARMM force field was applied to minimize the energy of ligands structure. For protein preparation, the X-ray crystal structures of human FosB/JunD bZIP domain (PDB ID: 5VPB) were downloaded from Protein Data Bank (PDB) database (http://www.rcsb.org/pdb/home/home.do). Before docking, the structure of this domain was prepared by adding up the hydrogen atoms and removing water molecules. Ten molecular docking poses saved for each ligand were ranked according to -CDOCKER energy (Supplementary Tables S3, S4). The lowest energy structure was considered as the most stable conformation, thus, the optimum pose with the highest -CDOCKER interaction energy was selected to detect the interaction between macrolide compounds and FosB/JunD bZIP domain.

Zebrafish larvae at three dpf were exposed to azithromycin, and its impurities in 2 mM and dissected liver of larvae were collected at six dpf for transcriptome sequencing. As previously described (Han et al., 2018a), the process of ribonucleic acid sequencing (RNA-seq) analysis was carried out by CapitalBio Corporation (Beijing, China). To identify the differentially expressed genes (DEGs) between the control and compound treatment group, the threshold of fold change ≥2 (upregulated gene) were established. Functional enrichment analysis was performed by Gene Ontology (GO) terms and Kyoto Encyclopedia of Genes and Genomes (KEGG) pathways under a significance threshold of p-value < 0.05. Cytoscape software was employed to map the drug–pathway network.

As previously described by Zhang et al. (2020), the treated zebrafish larvae liver and human L02 cells were collected for the qRT-PCR analysis. Total RNA were extracted from L02 cells and the dissected liver of larvae by using TRIzol reagent (Sigma–Aldrich, Saint Louis, MO, United States) and synthesized cDNA using M-MLV reverse transcriptase (Promega, M1708). The reaction of qRT-PCR was performed by LightCycler480. The relative mRNA levels were calculated with β-actin as the reference gene using the comparative 2^ΔΔCT quantization method. The specific primers (5′-3′) for amplification in zebrafish are as follows: fosab (F: TTA​CCA​GCC​TTA​ACG​CCG​AC, R: TGG​ACC​ATC​CAC​TGC​AAG​TC), β-actin (F: CCG​TGA​CAT​CAA​GGA​GAA​G, R: ATA​CCG​CAA​GAT​TCC​ATA​CC). The primer sequences (5′-3′) for qRT-PCR in L02 cells are as follows: c-fos (F: TAC​TAC​CAC​TCA​CCC​GCA​GA, R: GGC​CTC​CTG​TCA​TGG​TCT​TC), β-actin (F: CAC​CAT​TGG​CAA​TGA​GCG​GTT​C, R: AGG​TCT​TTG​CGG​ATG​TCC​ACG​T).

As previously described by Zhang et al. (2020), the human L02 cells were incubated with each macrolide compounds at different concentrations for 24 h. For cell viability, Cell Counting Kit-8 (CCK-8; Dojindo, CK04) was used by following the manufacturer’s protocol and the absorbance of each well was measured at 450 nm using a microplate reader. For cell toxicity, the level of lactate dehydrogenase (LDH) release was determined by LDH activity assay kits (Dojindo, CK12) according to the manufacturer’s specification, and the absorbance at 490 nm in each well was measured by a microplate reader.

All experiments were performed in triplicate. All the data were shown as mean ± SD, and statistical analyses were performed by GraphPad Prism software. The statistical comparisons were carried out by one-factor ANOVA, and p-values < 0.05 were considered significant.

Considering that c-fos is involved in some drug’s hepatotoxicity and its basic leucine zipper domain (bZIP) contributes to AP-1 binding to DNA (Zhang et al., 2020; Ellenberger et al., 1992; Landschulz et al., 1988), we infer that the interaction between the FosB/JunD bZIP domain and structures of drugs probably is a key access for understanding hepatotoxicity of macrolide drugs. Here, we used the FosB/JunD bZIP domain as c-fos receptor domain for macrolide molecular docking analysis. We selected erythromycin, roxithromycin, clarithromycin, azithromycin, midecamycin, josamycin, and acetylspiramycin as ligands to docking into the active site of the FosB/JunD bZIP domain. The most optimum docking poses in this domain are shown in Figure 1, and the highest -CDOCKER interaction energy scores are given in Table 1. The results showed that the scores of these seven macrolides were closely. -CDOCKER interaction energy score implies that the binding affinity between the ligand and the protein (Rampogu and Lemuel, 2016). The affinity order is azithromycin > erythromycin > acetylspiramycin > roxithromycin > midecamycin > clarithromycin > josamycin, which is in good agreement with the reports in literatures. The 16-membered ring macrolides were less hepatotoxic than erythromycin (Periti et al., 1993). Roxithromycin and clarithromycin have less hepatotoxicity than erythromycin, and azithromycin has higher hepatotoxicity than erythromycin (Björnsson, 2017). Our previous experimental data indicated that the liver toxicity induced by these seven compounds were similar (Zhang et al., 2020).

Then, we further confirmed the correlation between the docking energy scores and liver toxicity caused by macrolide antibiotics. It is well known that telithromycin is a ketolide antibiotic derived from macrolide structure (Supplementary Figure S1), has been pulled from the market due to severe acute liver injury (Author Anonymous, 2012). Therefore, we investigated the liver toxicity of telithromycin compared to azithromycin. First, employing transgenic zebrafish Tg (fabp10: dsRed), the liver phenotype in vivo showed that telithromycin could cause hepatomegaly at 0.1 mm and a dark coloration and amorphous liver at 0.5 mm; a dark and amorphous liver indicates degeneration or necrosis of liver cells. While azithromycin showed no toxicity to zebrafish liver at 0.1 mm, but hepatomegaly at 0.5 mm (Figure 2A). Thus, telithromycin has a higher liver injury than azithromycin in zebrafish. Second, the CCK-8 assay and LDH release assay were used to detect the cytotoxicity of telithromycin and azithromycin on human normal liver cells (L02 cells) (Figure 2B), the results suggested that telithromycin exerted a higher cytotoxicity than azithromycin in the human liver cells. Finally, the expression of fosab in zebrafish and c-fos in L02 cells were tested by using qRT-PCR analysis (Figure 2C), the results showed that telithromycin significantly elevated the mRNA levels of fosab and c-fos than azithromycin. These data demonstrated that the liver toxicity of telithromycin was higher than that of azithromycin.

Meanwhile, telithromycin as a ligand docked into the active site of FosB/JunD bZIP domain, and the docking score was 66.31 kcal/mol, which was much higher than that of azithromycin (Supplementary Figure S2; Supplementary Table S1). Therefore, the highest -CDOCKER interaction energy score with FosB/JunD bZIP domain may reflect the hepatotoxicity induced by macrolide antibiotics.

To quickly assess the liver toxicity of diverse impurities of macrolides compared with their active pharmaceutical ingredients (APIs) in a time- or cost-saving manner, we formulated a prediction interval for hepatotoxicity according to the docking scores of macrolides binding with the FosB/JunD bZIP domain. The scope of the prediction interval was the mean ± 2 SD of docking scores of the selected seven representative macrolides (erythromycin, roxithromycin, clarithromycin, azithromycin, midecamycin, josamycin, and acetylspiramycin) commonly used clinically, which is 48.1∼53.1. When the score is less than 48.1, it is considered that the chemical has a little or no liver toxicity (less than API); when the score is in 48.1∼53.1, the liver toxicity is approximately equal to API; when the score is greater than 53.1, the liver toxicity is greater than API, and this impurity requires strict quality control.

Therefore, we might be able to put forward a hypothesis for the hepatotoxicity assessment model of macrolides: the virtual screening of impurities with high liver toxicity could be carried out quickly according to this prediction interval, and the effect on c-fos gene expression in L02 cells can be used as a complementary verification, and then the positive can be controlled as a specific impurity.

In order to prove the aforementioned hypothesis and test the practicality of the hepatotoxicity assessment model of macrolides, first, we employed the prediction interval to virtually screen azithromycin impurities with high hepatotoxicity. We chose azithromycin and its impurities with different structures (impurities E, F, J, H, I, K, L, Q, R, and S) (Figure 3) as ligands to dock into the FosB/JunD bZIP domain. The most optimum docking poses of the azithromycin impurities are shown in Figure 4, and the highest -CDOCKER interaction energy scores are presented in Table 2. These scores revealed that the docking score of impurity H was the highest, followed by impurity L, both of which were greater than 53.1; the docking scores of impurities E, I, J, and Q were 48.1∼53.1; and the docking scores of impurities F, K, R, and S were less than 48.1. According to the prediction interval for hepatotoxicity of macrolides, we infer that impurities H and L might display higher hepatotoxicity than azithromycin; the hepatotoxicity of impurities E, I, J, and Q maybe equal to that of azithromycin, while the hepatotoxicity of impurities F, K, R, and S perhaps lower than that of azithromycin.

Next, we detected the expression of c-fos gene in L02 cells to verify the results of virtual screening. The impurities F, J, I, L, Q, R, and S were successfully synthesized except for impurities E, H, and K, and were examined by qRT-PCR analysis. The results showed that impurities I, J, L, Q, and R significantly upregulated the expression of c-fos gene in L02 cells at a concentration dependent manner, and impurity L had the greatest effect, impurity R moderately lower than azithromycin, while F and S did not affect (Figure 5). These were consistent with the virtual screening results. Therefore, the hepatotoxicity prediction of azithromycin impurities may be accurate, implying feasibility of the hepatotoxicity assessment model for macrolides.

To comprehensively analyze the hepatotoxicity induced by azithromycin impurities and further verify the results drawn by the assessment model, we investigated the liver toxicity of azithromycin impurities in vivo and in vitro. First, the ADMET properties of these impurities are predicted in Table 3, the values of topological polar surface area (TPSA) and logP showed that azithromycin impurities presented high lipophilicity, indicating these impurities might have good membrane permeability or oral absorption (Fonteh et al., 2015; Qidwai, 2016). Among them, impurity H has the highest TPSA and log P values. The absorption levels were assessed by water solubility and Caco2 permeability, the Caco2 permeability value of impurity L was the highest, and impurity H was the second. Therefore, impurity H and L had a better absorption than other impurities, which might lead to higher concentration in vivo. The volume of distribution at steady state (VDss) and blood–brain barrier membrane permeability (logBB) were used to assess the distribution level of these impurities, the VDss value of impurity K was the highest, and impurity H was the second. A high VDss value indicates that this compound is more inclined to accumulate in tissues. The logBB of all impurities was < −1, indicating these impurities did not easily cross the blood–brain barrier. The metabolism levels were evaluated by the two main subtypes of cytochrome P450, namely, CYP2D6 and CYP3A4. Impurities F, H, I, L, and Q were likely to be metabolized by CYP3A4, while all these impurities were not metabolized by CYP2D6. About the drug elimination level, the total clearance of impurity H was the lowest. For prediction of toxicity, all selected azithromycin impurities might cause hepatotoxicity.

Second, we detected the cytotoxicity of these impurities in L02 cells in vitro via CCK-8 assay to test the cell viability (Figure 6A). The results showed that all the impurities and azithromycin reduced the cell viability in a dose dependent manner. There were no significant differences in cytotoxicity of these impurities compared with azithromycin, except for impurity F and impurity S displaying less cytotoxicity in L02 cells, respectively. Furthermore, we chose impurity L, a potentially higher hepatotoxicity impurity, and impurity F, a possibly lower hepatotoxicity impurity, to assess the liver toxicity compared with azithromycin using the transgenic zebrafish (Figure 6B). The results demonstrated that impurity F was nearly non-toxic to zebrafish liver at 0.5 and 1 mm; azithromycin initiated hepatomegaly at 0.5 mm and blackened liver tissue at 1 mm; meanwhile, impurity L triggered hepatomegaly and caused the liver lost transparency with a brown coloration at 0.1 mm, and the liver tissue became black at 0.5 mm. These results demonstrated that impurity F and impurity S had a lower liver toxicity than azithromycin, whereas the liver toxicity of impurity L was higher than azithromycin. The results of liver toxicity in vitro and in vivo were consistent with the results evaluated by the assessment model. Importantly, these experimental data also provided proof for the accuracy and practicality of the assessment model for hepatotoxicity induced by macrolides.

Furthermore, RNA-seq and qRT-PCR analyses were performed to examine the expression of c-fos gene in vivo. According to GO enrichment analysis, we found fosab mRNA level were upregulated significantly by azithromycin impurities J, I, L, Q, and R, respectively, while impurity F and S had no significant effect on the expression of fosab in zebrafish (Table 4). The fold change of fosab gene expression promoted by impurity L was the highest. The results of qRT-PCR analysis also demonstrated that these impurities could significantly elevate the expression of fosab gene in zebrafish liver at a concentration dependent manner except for impurities F and S, while impurity L has the greatest effect (Figure 6C). These data were similar to the expression of c-fos gene in L02 cells. On the other hand, our previous study has shown that the liver toxicity induced by macrolides was related to the high expression of c-fos (Zhang et al., 2020), these data might explain why impurity L has a higher liver toxicity, probably due to its greater effect on c-fos gene expression. Therefore, the assessment model is accurate, which can be illustrated by the experimental confirmation via detection of the expression of c-fos gene in L02 cells.

Furthermore, we performed comparative analyses of the structures of azithromycin impurities on the basis of the structural characteristics (Figure 3). The structures of these impurities have mainly changed in the side chain at C3, C5, and C9 position of lactone rings. Among them, the structures of impurities F, L, and H are different only in the nitrogen (N)-linked substituents on the side chain group at C5 position compared with azithromycin. The N-linked substituents in impurity F is an aldehyde group, which is an electron-withdrawing group from N atom, while the N-connected substituents in impurity H and impurity L are electron-donating groups, which lead to different charges of N atoms (Figure 7). Consequently, the hepatotoxicity of azithromycin impurities may be related to the charge of N atoms on the side chain group at C5 position through the analysis of the structure–hepatotoxicity relationship of azithromycin impurities. Simultaneously, the structural modification of C9 position to form lactone or oximido cannot exacerbate the hepatotoxicity according to the structures of impurities Q, R, and S, and the oximido of impurity S dramatically reduces hepatotoxicity. Similarly, the changes at C3 side chain also cannot aggravate the hepatotoxicity according to the structures of impurities K and J.

All the aforementioned results indicate that the assessment model for evaluating hepatotoxicity of macrolides may be accurate and practical, besides, the charge of N atoms on the side chain group at C5 position could also reflect the hepatotoxicity of macrolides.