LLC-PK1 (ATCC CL101) was obtained from the American Type Culture Collection (Manassas, VA, USA). Cells were cultured at 37°C with 5% CO₂ in Dulbecco’s modified Eagle’s medium (Hyclone, Logan, UT, USA) with 8% Fetal Bovine Serum, 100 U/mL penicillin, and 100 μg/mL streptomycin (Welgene, Daegu, Korea). Treatments with GABA and PIP (Sigma Aldrich, St. Louis, MO, USA), MAPK inhibitors SB 202190 and SP600125 (Calbiochem, San Diego, CA, USA) was performed when cells reached 80% to 90% confluence. Cells were incubated at 37°C with 5% CO₂ in DMEM without fetal bovine serum for 1 hour before treatments with 0.1 mM GABA and/or 5 μM PIP in final concentration.

Prepared cells were harvested and RNAs were isolated using TRIzol reagent (Ambion, Carlsbad, CA, USA) or RNeasy min kit (Qiagen, Germantown, MD, USA). The purity and integrity of total RNA were checked using NanoDrop One (Thermo Scientific, Waltham, MA, USA).

cRNA was prepared from 1 to 5 μg total RNA using Agilent’s Quick Amp labeling kit (Agilent, Santa Clara, CA, USA). Labeling of probes, in vitro transcription, hybridization, and dyeing of chips were performed by Macrogen Gene Company Limited (Seoul, Korea). Then 1.65 μg of cRNA was hybridized to Agilent Porcine Gene Expression microarray (43,803 probes). Arrays were scanned for data export and processing using Agilent Feature Extraction v11.0.1.1 software on R 2.15.1. Shortly after, signal intensity was extracted. Quality check and filtering were performed by flag. Probe signal intensity was converted to log (base 2) value. Normalization for untreated, GABA, PIP, or both was performed using quantile normalization. Fold changes (FC) compared to untreated sample were calculated. Probes that yielded >1.5 or <−1.5 folds change with a p-value <0.05 in either treatment relative to untreated LLC-PK1 were retained, totaling 2,393 probes.

Heatmapper [13] was used to generate heatmap plot representing genes using average linkage clustering method with Pearson measurement method. To evaluate biological functions of differentially expressed gene, gene ontology analysis was performed using statistical analysis tool PANTHER Overrepresentation Test [14] (Released 20190308) with GeneOntology database released 2019-01-01 and Sus scrofa (all genes) as reference. Fischer’s exact test with Bonferroni correction for multiple testing (p<0.05) was performed with the set of validated genes from microarray data. Functional groups were ordered by their overall gene numbers and overall enrichment score [15].

To validate selected differentially expressed genes, real-time reverse-transcription polymerase chain reaction (RT-PCR) was performed using CFX96 Real-Time System (BIO-RAD, Hercules, CA, USA). Briefly, 1 μg of RNA was treated with DNA-free Kit (Thermo Scientific, USA) following the manufacturer’s protocol. Subsequently, first-strand DNA synthesis was conducted with oligo [13] and AccuPower CycleScript RT PreMix kit (BIONEER, Daejeon, Korea). Amplification was done using IQ SYBR Green Supermix (BIO-RAD, USA). Primers used are listed in Table 1. The results were reported as FC compared to control expression normalized to glyceraldehyde 3-phosphate dehydrogenase by the Pfaffl method [16]. Pairwise comparisons between means of different treatments were performed using a two tailed Student t-test, the null hypothesis being that the means are equal.

Differential expressed genes heatmap has been plotted in Figure 1A. For validation purpose, FCs of genes of interest were retrieved as differential signals from multiple probes might reflect differential isotype expression upon treatment [17]. Probes were located and mRNA isotypes of the genes of interest were analyzed using UCSC genome Browser (http://genome.ucsc.edu/) based on Sus Scrofa genome (Sscrofa11. 1/susScr11) (Figure 1B). High mobility group protein B2 (HMGB2) was represented by five probes. However, validation showed that only four probes corresponded to HMGB2. Three of these probes were located in the 3′ UTR and one probe was located in the third exon of identified HMGB2 isotypes. Of these four probes, A_72_P675158 was upregulated by GABA and A_72_P077771 was synergistically upregulated by GABA and PIP at significant levels while the remaining two showed synergistic upregulation at lower magnitude (FC<1.5). LIM domain-binding protein 1 (LDB1) was represented by two probes. A_72_P540123 located in the seventh exon of AB242619 or the second exon of AK344041 showed upregulation upon GABA treatment. On the other hand, A_72_P185256 located in the 5′ UTR of AB242619 showed no significant change upon GABA and/or PIP treatment. Only one probe for Wilms tumor protein homolog (WT1) was available. A_72_P077696 located in the 3′ UTR of the longest isotype showed upregulation by GABA alone and in combination with PIP. To validate these results, we designed quantitative RT-PCR (RT-qPCR) primers in the 3′UTR of reference sequences of HMGB2, LDB1, and WT1 according to NCBI RNA reference sequences collection (RefSeq). After GABA and/or PIP treatment at the same condition as for microarray analysis, mRNAs were extracted and RT-qPCR was performed to evaluate expression levels of HMGB2, LDB1, and WT1 after GABA and/or PIP treatment compared to untreated cells (Figure 1C). HMGB2 and LDB1 showed increase by GABA treatment as A_72_P675158 and A_72_P540123 probes. WT1 showed expression pattern that concurred with its microarray probe.

From microarray data, we identified 302 genes upregulated by GABA alone and 26 genes upregulated by PIP alone (FC >1.5). There were 202 downregulated genes by GABA alone and 43 downregulated genes by PIP alone. Synergistically upregulated genes were defined by two criteria. First, we filtered for genes with FCs greater than 1.5 upon GABA and PIP co-treatment. Second, among these genes, we selected 227 genes that has higher FC after co-treatment compared to treatment by GABA alone or PIP alone. Twenty genes were defined as synergistically downregulated as FCs after co-treatment were <−1.5 while FCs after GABA or PIP treatment were between 1.0 and 1.1. GO analysis was performed for genes upregulated by GABA, GABA/PIP, and downregulated by GABA alone. We were unable to perform GO analysis for lists of genes affected by PIP treatment and genes downregulated by GABA/PIP in combination due to low number of genes in these lists (Figure 2). GABA-upregulated genes were related to metabolic process (31.5%), biological regulation (26.2%), cellular process (15.6%), localization (12.3%), and cellular component organization or biogenesis (11.3%). GABA- and PIP- upregulated genes were related to metabolic process (30.4%), biological regulation (23.2%), cellular component organization or biogenesis (11.9%), and cellular process (9.3%). Genes downregulated by GABA alone were related to metabolic process (31.5%), biological regulation (24.9%), response to stimulus (13.0%), and cellular process (11.4%).

Fold enrichment analysis can evaluate the number of genes that represent a specific pathway versus the number expected at random case [14]. At a p-value cut-off of 0.05, 15 biological processes were enriched at >1.89- fold for GABA upregulated genes, 76 biological processes were enriched at >2.00-fold for GABA downregulated genes, and 15 genes were enriched at >1.80-fold for synergistically upregulated genes (Figure 3A). The most upregulated biological pathway by GABA alone treatment was erythrocyte differentiation, while the majority of downregulated biological pathways were involved in responses to various stimulus. The most upregulated biological pathway in synergistically upregulated was related to reproductive development. In contrast to fold enrichment higher than one, fold enrichments that were lower than one represented biological processes that were underrepresented and therefore putatively downregulated in the gene list of interest (Figure 3B). At p-value cut-off of 0.05, nine biological processes were enriched at <0.18-fold for GABA upregulated genes and seven biological pathways were enriched at <0.10-fold for GABA/PIP upregulated gene. Both lists showed underrepresentation of biological processes related to stimulus.

Based on microarray data and enrichment analysis, we identified seven genes that were upregulated by GABA and related to biological process controlling positive regulation of erythrocyte differentiation. We also identified 31 genes that were synergistically upregulated and related to immune systems and processes. FCs of these genes of interests were retrieved from microarray data (Table 2). Validation was performed for inhibin subunit beta A (INHBA), signal transducer and activator of transcription 1 (STAT1), hypoxia-inducible factor 1-alpha (HIF-1alpha), and SPF36 ring finger protein like 1 (ZPF36L1) among GABA upregulated genes. Our previous study has shown that EPO and EPO-R upregulation by GABA/PIP in combination can be downregulated by inhibitor of JNK (SP600126) and inhibitor of p38 (SB203580) MAPK [11]. Therefore, we investigated if MAPK inhibition could affect GABA upregulated genes (Figure 4A). INHBA, STAT1, and HIF1A were validated to be upregulated by GABA. Inhibition of JNK activation significantly downregulated expression of these genes. However, inhibition of p38 activation affected expression of these genes at lesser extent. Validation was performed for tripartite motif containing 8 (TRIM8), GATA binding protein 3 (GATA3), BPI fold containing family A member 1 (BPIFA1), and C-C motif chemokine ligand 4 (CCL4) among synergistically upregulated genes involved in biological process related to immune system process. Expression analysis upon addition of inhibitors was also performed for those genes (Figure 4B). Results showed that TRIM8, GATA3, BPIFA1, and CCL4 were synergistically upregulated by GABA and PIP in combination, while the addition of MAPK inhibitors downregulated their expression.