Global Gene Networks in 3D4/31 Porcine Alveolar Macrophages Treated with Antigenic Epitopes of Actinobacillus pleuropneumoniae ApxIA, IIA, and IVA

Actinobacillus pleuropneumoniae (App) is the causative agent of porcine pleuropneumonia. Although App produces several virulence factors, Apx toxins, the primary App virulence factors, have been the focus of numerous studies. However, the host response against the Apx toxins has not been elucidated at the transcriptomic level. Therefore, in this study, we examined the response of an immortalized porcine alveolar macrophage cell line (IPAM 3D4/31) to four antigenic epitopes of the App exotoxins, ApxIA, IIA and IVA. The antigenic epitopes of the Apx toxins (ApxIA Ct, ApxIIA Nt, ApxIVA C1 and ApxIV C2) were determined by an in-silico antigenicity prediction analysis. Gene expression in IPAMs was analyzed by RNA-Seq after treatment with the four proteins for 24 h. A total of 15,269 DEGs were observed to be associated with cellular and metabolic processes in the GO category Biological Process and nuclear receptors and apoptosis signaling in IPA analyses. These DEGs were also related to M2 macrophage polarization and apoptosis in IPAMs. These host transcriptional analyses present novel global gene networks of the host response to treatment with Apx toxins.

www.nature.com/scientificreports www.nature.com/scientificreports/ DEGs associated with several signaling pathways including "Intracellular and Second Messenger Signaling" (sirtuin signaling pathway, insulin receptor signaling, and the unfolded protein response), and "Cellular Stress and Injury" (role of BRCA1 in DNA damage response, NRF2-mediated oxidative stress response, and ATM signaling). Significant canonical pathways in ApxIVA C1-treated cells were enriched in phospholipase C signaling, protein kinase A signaling, and the sirtuin signaling pathway in "Intracellular and Second Messenger Signaling" and for iNOS signaling, PI3K signaling in B lymphocytes, and p38 MAPK signaling in "Cellular Immune Response". The canonical pathways in ApxIVA C2-treated IPAMs were primarily associated with actin nucleation by the ARP-WASP complex, integrin signaling and glucocorticoid receptor signaling in "Intracellular and Second Messenger Signaling" and granulocyte adhesion and diapedesis, granulocyte adhesion and diapedesis, and IL-10 signaling in "Cellular Immune Response".
Analysis of signaling pathways in the Apx toxins-treated IpAMs. According to the signaling pathway classification, the DEGs of IPAMs were commonly enriched in the "Apoptosis", "Nuclear Receptor Signaling", "Cellular Immune Response" and "Humoral Immune Response" categories. Based on the z-scores, apoptosis signaling, PTEN signaling involved in "Apoptosis" and LXR/RXR activation, PPAR signaling related to "Nuclear Receptor Signaling" were commonly predicted to be activated, while B cell receptor signaling and HMGB1 signaling involved in the "Humoral Immune Response" and NF-κB signaling and IL-6 signaling related to the "Cellular Immune Response" were commonly predicted to be inactivated (Fig. 3).
Among the signaling pathways of IPAMs treated with ApxIA Ct, ApxIIA Nt and ApxIVA C2, macropinocytosis signaling, MIF regulation of innate immunity and CD40 signaling related to the "Cellular Immune Response" were commonly expressed. However, micropinocytosis signaling and MIF regulation of innate immunity were commonly predicted to be inhibited, although CD40 signaling was predicted to be activated only in the ApxIVA C2-treated cells according to the z-score. The DEGs also significantly mapped to the signaling pathways related to "Apoptosis", comprising toll-like receptors signaling, induction of apoptosis by H1V1 and TWEAK signaling. Induction of apoptosis by H1V1 was commonly predicted to be activated in the ApxIA Ct-, ApxIIA Nt-and ApxIVA C2-treated cells.

Signaling pathways of Apx toxins-treated IPAMs associated with M2 macrophage polarization.
An analysis of significant signaling pathways associated with macrophage polarization was carried out ( Table 1). The results indicated that glucocorticoid receptor signaling, GM-CSF signaling, IL-10 signaling, IL-6 signaling, LXR/RXR activation, NF-κB activation by viruses, NF-κB signaling, PI3K signaling in B lymphocytes, PI3K/AKT signaling and PPAR signaling were commonly expressed and showed consistent activation states in cells treated with ApxIA Ct, ApxIIA Nt, ApxIVA C1 and ApxIVA C2. Among these pathways, LXR/RXR activation and PPAR signaling associated with M2 polarization were commonly activated. However, GM-CSF signaling, IL-6 signaling, NF-κB activation by viruses and NF-κB signaling associated with M1 polarization were suppressed or exhibited a z-score of 0, in the IPAMs treated with ApxIA Ct, ApxIIA Nt, ApxIVA C1 and ApxIVA C2. In particular, NF-κB, IL-1, and ERK were down regulated in the signaling pathways of IPAMs treated with the four Apx recombinant proteins. These molecules can induce the activation of PPARγ, which results in M2 differentiation. PPARγ was commonly predicted to be activated according to the observed z-score for PPAR signaling in ApxIA Ct-, ApxIIA Nt-and ApxIVA C1-treated cells. However, PPARγ was predicted to be inhibited in the ApxIVA C2 treatment group, because NCOR was predicted to be activated based on the z-score (Fig. 4). The real-time PCR results Activation of apoptosis signaling in IPAMs treated with ApxIA Ct, ApxIIA Nt, and ApxIVA C2.
An analysis of signaling pathways involved in apoptosis was performed (Fig. 5). In total, 18, 4, 5, and 5 significant pathways [(−log(p-value) 1 3 ≥ . ] associated with apoptosis were observed in IPAMs treated with ApxIA Ct, ApxIIA Nt, ApxIVA C1, and ApxIVA C2, respectively. Among the signaling pathways of ApxIA Ct-treated cells, apoptosis signaling, PTEN signaling, death receptor signaling, induction of apoptosis by H1V1, and ceramide signaling were predicted to be activated. Of the signaling pathways expressed in ApxIIA Nt-treated cells, apoptosis signaling, ceramide signaling, and PTEN signaling were also predicted to be activated. However, ApxIVA C1 did not predict any predicted pathways associated with "Apoptosis". Among the five signaling pathways in ApxIVA C2-treated cells related to "Apoptosis", toll-like receptor signaling and apoptosis signaling were predicted to be activated. In particular, apoptosis signaling was commonly predicted to be activated and was more significant [−log(p-value) of 6.38, 2.24 and 2.05] than other significant pathways in the IPAMs treated with ApxIA Ct, ApxIIA Nt, and ApxIVA C2, respectively.
Aif, Apaf1, Bad, Bak, Bcl2, Bcl-xl, Caspase 2, Caspase 3, Caspase 6, Caspase 7, Gas2, LaminA and p90rst were the significantly expressed apoptosis signaling genes in ApxIA Ct-treated cells. The genes Aif, Bad, Bak, Bcl2, Caspase 6 and Caspase 7 were significantly expressed in the apoptosis signaling pathway of ApxIIA Nt-treated cells. Bcl2, Caspase 6 and Caspase 7 were significantly expressed in the apoptosis signaling pathway of ApxIVA C2-treated cells. No DEGs in the apoptosis signaling pathway of ApxIVA C1-treated cells were significantly expressed. Among the genes expressed in the apoptosis signaling pathway of ApxIA Ct-, ApxIIA Nt-and ApxIVA C2-treated cells, Caspase 6 was commonly up-regulated and Caspase 7 was commonly down-regulated. The genes Bak, Cycs, Apaf1, Caspase 9, Caspase 6 and LaminA were predicted to be activated in IPAMs treated with ApxIA Ct and ApxIIA Nt. The genes Caspase 9, Caspase 3, Caspase 6 and LaminA were predicted to be activated in ApxIVA C2-treated cells (Fig. 6). Apaf1, Bak, Caspase6, and Lamin A, which are known as key factors of apoptosis, were identified as significant genes in this pathway. In the real-time PCR results, genes (bak, bax, cytochrome c, apaf1, caspase 3, caspase 6 and laminA) related to apoptosis were significantly up regulated after 24 h of infection, except in ApxIVA C1-stimulated IPAMs ( Supplementary Fig. 4).

ApxIA Ct
ApxIIA Nt ApxIVA C1 ApxIVAC2 −log(p-value) Ratio z-score −log(p-value) Ratio z-score −log(p-value) Ratio z-score −log(p-value) Ratio z-score www.nature.com/scientificreports www.nature.com/scientificreports/ Validation of RNA-seq by qRt-pCR. RNA samples were assayed by quantitative real-time RT-PCR (qRT-PCR) to validate the RNA-Seq results. The IPAMs treated with Apx toxins showed significant changes in gene expression compared to those in DPBS-treated cells. The up-regulated gene hmgcs1, the down regulated genes, nts, mgp and abcg1 and the randomly expressed vgf gene were selected to assay by qRT-PCR using the same RNA samples as those used for the RNA-Seq analysis. The correlation coefficient between the two analyses was 0.906. Although the RNA-Seq data for some genes such as abcg1 and vgf showed greater fluctuations (increase or decrease) compared to those of the qRT-PCR data, the differential expression of all selected genes was validated, as the qRT-PCR results showed the same trend, with respect to up-regulation or down-regulation ( Supplementary  Fig. 2).

Discussion
Transcriptional studies of App-infected hosts have been conducted under in vivo and in vitro conditions to identify App infection mechanisms in the host 31,34,36,37 . However, a transcriptomic analysis to identify the effects of Apx toxins in the host had not been performed. Therefore, the goals of this study were to identify the invasive mechanisms of Apx toxins in an immortalized porcine alveolar macrophage cell line using partial epitopes of three types of Apx toxins. Transcriptomic analyses were performed using IPAMs treated with partial epitopes of ApxIA, ApxIIA, and ApxIVA. The results suggest possible roles for ApxIA Ct, ApxIIA Nt, and ApxIVA C2 in impairing the host defense system through the induction of apoptosis in IPAMs.
In this study, gene expression was analyzed in IPAMs treated with ApxIA Ct, ApxIIA Nt, ApxIVA C1, or ApxIVA C2. Among the identified signaling pathways, the DEGs of the ApxIA Ct-treated IPAMs were associated with "Apoptosis" and "Intracellular and Secondary Messenger Signaling", those of ApxIIA Nt-treated cells were associated with "Intracellular and Second Messenger Signaling" and "Cellular Stress and Injury", and those of ApxIVA C1-and ApxIVA C2-treated cells were enriched in "Intracellular and Second Messenger Signaling" and "Cellular Immune Response". Signaling pathways of IPAMs treated with ApxIA Ct, ApxIIA Nt, ApxIVA C1, and ApxIVA C2 were commonly enriched in "Apoptosis", "Nuclear Receptor Signaling", "Humoral Immune Response" and "Cellular Immune Response". Apoptosis signaling and PTEN signaling associated with "Apoptosis", LPS/ www.nature.com/scientificreports www.nature.com/scientificreports/ IL-1 mediated inhibition of RXR function, LXR/RXR activation and PPAR signaling associated with "Nuclear Receptor Signaling" were predicted to be activated, while B cell receptor signaling and HMGB1 signaling associated with "Humoral Immune Response", GM-CSF signaling, IL-6 signaling, NF-κB signaling and PI3K/AKT signaling associated with "Cellular Immune Response" were predicted to be inhibited.
The activation states of macrophages are classified into classically activated (M1) and alternatively activated (M2), and these phenotypes are affected by different cytokines and microbial products 38 . M1 macrophages produce pro-inflammatory cytokines, such as TNF-α, IL-6, and IL-12. In contrast, M2 macrophages release anti-inflammatory cytokines, such as transforming growth factor (TGF)-β, IL-10, and IL-1 receptor antagonist [22][23][24] . In our study, polarization was predicted by the inhibition of pathways associated with M1 macrophage differentiation and the activation of pathways associated with M2 macrophage differentiation. In addition, nuclear receptor signaling was significantly enhanced in cells treated with the three Apx antigenic epitopes. The pathways associated with M1 macrophage polarization, including GM-CSF signaling, IL-6 signaling NF-κB activation by viruses and NF-κB signaling, were commonly predicted to be inhibited or had a z-score of 0 in the four antigenic epitope-treated IPAMs. However, the pathways associated with M2 macrophage polarization including LXR/RXR activation and PPAR signaling, were predicted to be activated. Of these signaling pathways, PPARγ was commonly predicted to be activated in the PPAR signaling pathway as a result of treatment with three Apx antigenic epitopes but not ApxIVA C2. In the PPAR signaling of ApxIVA C2-treated cells, PPARγ was predicted to be inhibited, because NCOR was predicted to be activated. However, inactivation of NF-κB and IL-1 can induce activation of PPARγ and these pathways were commonly down regulated with respect to PPAR signaling in IPAMs treated with the four antigenic epitopes. STATs, NF-κB and PPARs, known as leading transcription factors, drive a particular phenotype of macrophage polarization. Additionally, PPARγ, TGF-β, STAT3, IL5 and STAT6, which are related to alternatively activated macrophages were significantly up-regulated after 24 h of stimulation with the four antigenic epitopes according to the quantitative real-time PCR results.
Nuclear receptors such as peroxisome proliferators activated receptor (PPAR) can induce the biological effects caused by toxins. PPARs function as important regulators of cell differentiation, proliferation, and apoptosis in macrophages. Importantly, of the various PPAR subtypes, PPARγ inhibits the expression of inflammatory cytokines and induces anti-inflammatory effects in the lung [39][40][41][42] . The importance of PPARγ in regulating the M2 phenotype of macrophages has been demonstrated by the fact that M2-type responses were compromised in the absence of PPARγ 26,27 . PPARγ expression is also important for the expression of characteristic genes in M2 macrophages 27,28 . These PPAR pathways induce anti-inflammatory responses to inhibit immune inflammatory genes (encoding IL-2, IL-6, IL-8, TNF-α, and metalloproteases) through the inhibition of NF-κB expression [43][44][45][46][47][48] . In particular, PPARγ can repress the activation of inflammatory response-associated genes by a SUMOylation-dependent pathway in macrophages, which induces the transrepression of NF-κB target genes 49 . The prediction of down-regulated NF-κB and up-regulated PPARγ in the PPAR pathway of IPAMs treated with the four Apx epitopes indicated that these pathways may be associated with M2 macrophage differentiation, and the inhibition of cellular immune responses in the analysis of signaling pathways also supports these results. www.nature.com/scientificreports www.nature.com/scientificreports/ Among the signaling pathways observed to be associated through the Ingenuity Pathway Analysis, the apoptosis signaling pathway was significantly enriched in IPAMs treated with the ApxIA Ct, ApxIIA Nt, and ApxIVA C2 epitopes. However, all pathways associated with apoptosis were not enriched in ApxIVA C1-treated cells. The genes caspase 9, caspase 6, and laminA were commonly activated in the apoptosis signaling pathway in IPAMs treated with ApxIA Ct, ApxIIA Nt, and ApxIVA C2. In addition, the genes bak, cycs, and apaf1 were activated in ApxIA Ct-and ApxIIA Nt-treated cells. We identified apaf1, bak, caspase6, and laminA as genes in the apoptosis pathway exhibiting significant differential expression in response to these Apx epitopes. Furthermore, bak, bax, cytochrome c, apaf1, caspase3, caspase6 and laminA were significantly up-regulated after 24 h of stimulation with ApxIA Ct, ApxIIA Nt and ApxIVA C2.
Several studies have shown that toxins can kill target cells by inducing caspase activity [50][51][52][53][54] , in particular, toxins stimulate intracellular signals that induce apoptosis by intrinsic signaling pathways. Bak, a proapoptotic protein, is known as a crucial protein in the induction of intrinsic apoptosis 55 . Activation of bak results in mitochondrial membrane permeabilization 56,57 , which causes cytochrome c, a proapoptotic mitochondrial protein, to be translocated to the cytosol 58 . Cytochrome c induces the assembly of Apaf1 as well as procaspase-9, forming an www.nature.com/scientificreports www.nature.com/scientificreports/ "apoptosome", which causes activation of the initiator caspase 9 59-61 . Activated caspase 9 causes DNA fragmentation and apoptosis through activation of caspase 3. Caspase 6 may be activated by the activation of caspase 3 and its activation is required for the formation of apoptotic bodies during apoptosis [62][63][64][65][66] . However, several studies showed that caspase 6 can activate procaspase 3, a downstream effector that causes activation of the cell death program 67,68 , and can act upstream of caspase 3 67 . In our study, caspase 6 was commonly activated in the apoptosis signaling of IPAMs treated with ApxIA Ct, ApxIIA Nt, and ApxIVA C2, while caspase 3 was only activated in ApxIIA Nt-, and ApxIVA C2-treated cells. However, Lamin A is activated by caspase 6, which eventually leads to apoptosis. Degradation of lamins, which are protease substrates, is mediated by the proteolysis of caspase 6. Lamins are a major component of the nuclear lamina and have nuclear functions during apoptosis 69,70 . In particular, lamin A/C plays an important role in nuclear apoptotic events such as in acting shrinkage, disassembly of the nuclear membrane and the formation of apoptotic bodies 62,63,71,72 . This analysis of signaling pathways showed that apoptosis through the intrinsic signaling pathway was commonly induced in IPAMs treated with ApxIA Ct, ApxIIA Nt, and ApxIVA C2. The genes bak, bax, cytochrome C, apaf1, caspase3, caspase6 and laminA were also commonly activated in the apoptosis signaling pathway of ApxIA Ct-, ApxIIA Nt-and ApxIVA C2-treated cells.
In conclusion, the results of this study indicated that the antigenic epitopes of three types of Apx toxins showed different immune responses in an immortalized porcine alveolar macrophage cell line. The signaling pathways were enriched in the categories Apoptosis, Nuclear Receptor, Humoral Immune Response, and Cellular Immune Response. Of these pathways, activation of apoptosis, nuclear receptor signaling and inactivation of both humoral and cellular immune responses were commonly predicted. In particular, the pathways associated with M2 macrophage polarization were activated, including PPAR signaling pathways. Among the PPAR signaling pathways, PPARγ was activated in IPAMs treated with the four antigenic epitopes. In addition, activation of LXR/ RXR and, PPAR signaling and inactivation of GM-CSF signaling, IL-6 signaling and NF-κB signaling were commonly predicted among the DEGs of IPAMs treated with the four Apx antigenic epitopes. This result indicates that the four antigenic epitopes of Apx toxins may induce M2 macrophage differentiation in IPAMs. The apoptosis signaling pathway was also commonly predicted to be activated in the IPAMs treated with ApxIA Ct, ApxIIA Nt, and ApxIVA C2. In particular, the Bak-APAF1-Caspase-9-Caspase-6-LaminA pathway was activated in the apoptosis of IPAMs induced by three Apx antigens. This result shows that ApxIA Ct, ApxIIA Nt and ApxIVA C2 may induce apoptosis in porcine alveolar macrophages and impair host defense systems. Future work will explore the role of the signaling pathways and the DEGs related to the pathogenesis of Apx exotoxin treatment to clarify the host response mechanism in porcine alveolar macrophages. This study will be helpful in understanding the host response to Apx toxins.
Transformed E. coli cells were cultured in LB broth with ampicillin at 100 µg/ml to an absorbance of 0.6 nm (OD 0.6nm ), at which time isopropyl β-D-1-thiogalactopyranoside (IPTG, 1 mM) was added and the cells were cultured for an additional 4 h. The cells were harvested and re-suspended in lysis buffer (20 mM Tris-HCl, 500 mM sodium chloride, 8 M urea, and 50 mM imidazole, pH 7.0). Next, nickel-nitrilotriacetic acid (Ni-NTA) chelate affinity chromatography was performed according to the manufacturer's instructions. The bound protein was eluted with elution buffer (20 mM Tris-HCl, 500 mM sodium chloride, 8 M urea, and 500 mM imidazole, pH 7.0). The purified recombinant protein was analyzed by sodium dodecyl sulfate-polyacrylamide gel electrophoresis (SDS-PAGE) and Western blot using an anti-histidine antibody (Jackson ImmunoResearch, USA) and an alkaline phosphatase kit (BioRad, USA) following the manufacturer's protocol.

RNA sequencing. After 24 h of incubation, total RNA was extracted from the cells using an RNeasy Mini
Kit (Qiagen, Germany) according to the manufacturer's instructions. For quality control, RNA purity and integrity were evaluated via denaturing gel electrophoresis, determining the OD 260/280 ratio, and analysis using an Agilent 2100 Bioanalyzer (Agilent Technologies, USA). The total RNA concentration was measured using Quant-IT RiboGreen (Invitrogen, USA). To determine the integrity of the total RNA, samples were evaluated by using TapeStation RNA ScreenTape system (Agilent Technologies, USA). Only high-quality RNA preparations, with an RIN greater than 7.0, were used for RNA library construction.
A library was prepared with 1 µg of total RNA for each sample using an Illumina TruSeq mRNA Sample Prep Kit (Illumina, USA). The first step in this workflow involves purifying the poly-A containing mRNA molecules using poly-T-attached magnetic beads. Following purification, the mRNA is fragmented into small pieces using (2019) 9:5269 | https://doi.org/10.1038/s41598-019-41748-3 www.nature.com/scientificreports www.nature.com/scientificreports/ divalent cations under elevated temperature. The cleaved RNA fragments are copied into first strand cDNA using SuperScript II reverse transcriptase (Invitrogen, USA) and random primers, followed by second strand cDNA synthesis using DNA Polymerase I and RNase H. These cDNA fragments then go through an end repair process, the addition of a single ' A' base, and subsequent ligation of the indexing adapters. The products are then purified and enriched via PCR to create the final cDNA library. The libraries were quantified using qPCR according to the qPCR Quantification Protocol Guide (KAPA Library Quantification Kits for Illumina Sequencing Platforms) and were qualitatively assessed using a TapeStation D1000 ScreenTape system (Agilent Technologies, USA). Indexed libraries were then sequenced using a HiSeq4000 platform (Illumina, USA) by Macrogen Incorporated.

Read mapping and gene expression analyses.
To estimate gene expression levels and to identify alternatively spliced transcripts, the RNA-Seq reads were mapped to the genome of Sus scrofa using TopHat 73 , which is capable of reporting split-read alignments across splice junctions, and evaluated using Cufflinks 74 with the default options. The transcript counts at the isoform level were calculated, and the relative transcript abundances were reported as FPKM (Fragments Per Kilobase of exon per million fragments mapped) using Cufflinks. In addition, novel transcript and novel alternative splicing transcripts were identified per sample. These results were obtained using the Cufflinks reference annotation based transcript (RABT) assembly method, which allows for the discovery of reference transcripts and novel transcripts using the-g option.
Raw data were calculated as FPKM (Fragments Per Kilobase of exon model per million mapped fragments) of each transcript in each sample using Cufflinks. We excluded the transcripts with zeroed FPKM values that were more than the total observed in the samples. Next, we added 1 to the FPKM values of the filtered transcripts to facilitate a log2 transformation. Filtered data were transformed by logarithm and normalized using the quantile normalization method. For each transcript, we conducted an independent t-test between the treatment and the control. Each transcript was determined using an independent t-test.
Screening and functional analysis of DEGs. DEGs were obtained by analyzing the mRNA transcripts of IPAMs treated with Apx toxins. Genes showing significantly altered expressions were subjected to gene set enrichment analysis using the PANTHER classification database (http://www.pantherdb.org). DEGs were categorized by biological processes and molecular functions using the PANTHER classification database and Fisher's exact test, to detect coordinated changes in predefined sets of related genes.
Biological system analysis. Data were analyzed through the use of IPA (QIAGEN Inc., https://www. qiagenbioinformatic s.com/products/ingenuitypathway-analysis) 75 . Differentially expressed genes with adjusted p values of less than 0.05 were uploaded into the IPA program. Each gene was mapped to its corresponding gene object in the Ingenuity Knowledge Base. A biological function analysis was performed using IPA to compare the DEGs associated with disease and disorders, molecular and cellular functions, and physiological system development and function in Apx toxin-treated IPAMs. A right-tailed Fisher's exact test was used to calculate the p-value for each biological function. Canonical pathways from the IPA library of canonical pathways were investigated to identify major biological pathways associated with the Apx toxin treatment of IPAMs. The significance of the association between the data set and the canonical pathway was determined based on two parameters, including (1) the ratio of the number of genes from the dataset that mapped to the pathway to the total number of genes that mapped to the canonical pathway and (2) the p-value calculated using Fisher's exact test to determine the probability that the association between the genes in the data set and the canonical pathway results from chance alone. expression analysis of selected genes by real-time pCR. After treatment with the four recombinant proteins for 2, 12, 24, and 48 h, total RNA was extracted from the cells using an RNeasy Mini Kit (Qiagen, Germany) according to the manufacturer's instructions. The expression of 22 genes was also studied by real-time PCR to verify the findings from the RNA-Seq data analysis or to acquire additional information about major immune responses (qRT-PCR, Supplementary Table 6). The qRT-PCR reactions were performed using 1 µL of cDNA, a Rotor-Gene SYBR Green PCR Kit (Qiagen, Germany) and a Rotor-Gene Q real-time PCR cycler (Qiagen, Germany). The cycling parameters were as follows: 95 °C for 3 min, followed by 45 cycles of 95 °C for 15 sec, 30 sec at 60 °C and 30 sec at 72 °C with fluorescence detected during the extension phase. The expression level was determined via the 2 −ΔΔCt method using glyceraldehyde-3-phosphate dehydrogenase (GAPDH) as a reference gene. The relative expression level was compared to that observed in control cells to determine the fold change in expression for each gene. statistical analysis. Statistical significance of internalization was analyzed by Student's t-test or repeated measures ANOVA using GraphPad Prism version 7.00 for Windows, GraphPad Software, La Jolla California USA, www.graphpad.com. Genes were considered differentially regulated when p was < 0.05. When differences were determined to be significant, fold changes with respect to the control condition were represented as follows: fold change = the mean ratio of gene expression in toxin-treated cells/the mean ratio of gene expression in DPBS-treated cells.

Data Availability
Raw files and normalized datasets are available from the Gene Expression Omnibus (GEO) https://www.ncbi. nlm.nih.gov/geo under the accession number GSE116263.