Biochemical properties of cofactor and coenzyme metabolism in porcine oviductal epithelial cells – a microarray study

The oviduct represents a key organ in the ultimate mutation of oocytes, transport of gametes, sperm capacitation, fertilization, and early embryo development. The mucosal layer of the oviduct can be exposed to various pathogens and endotoxins entering from the peritoneal cavity, follicular fluid, and uterus. Moreover, there is evidence showing that the presence of gametes in the oviduct or follicular granulosa cells, with associated fluid, can modify gene expression of the oviductal epithelial cells. Close to the time of ovulation, the oviduct epithelium is sensitive to local programming stimuli from both male and female gametes, follicular cells, and fluid [1,2]. Well characterized in vitro models of the oviduct epithelium are of interest for several general reasons. Firstly, they contribute to the understanding of molecular mechanisms associated with oocyte fertilization and early embryo development; additionally, they provide a tool for reproductive toxicity testing.

Ethical concerns and the lack of healthy tissue sources present limitations for the in vitro model of the human oviductal epithelium [3]. Regarding the anatomical structure and physiological processes, pig reveals some close similarities to human. Porcine derived tissues are readily available from local meat processing facilities, and the procedure for the isolation of epithelial cells is relatively simple and inexpensive. Due to these benefits, pig became one of the favored animal models for the study of the oviductal epithelial cells, it has already been used in studies of reproduction biology [4, 5, 6]. The porcine oviductal epithelial cells culture system is easily reproducible, and it represents a long-term culture in the strict sense. Porcine oviductal epithelial cells culture has been reported to be maintained for up to 6 weeks or 10 to 15 passages of cell lines [3].

Affymetrix microarray assays have been previously employed to investigate the expression profiles of genes involved in several important processes in oviductal epithelial cells, such as biological adhesion [7], amino acids metabolism and degradation [8], cell cycle process, proliferation [9], epithelium morphogenesis and oviduct development [10]. The aim of the present study is to further characterize the porcine oviductal epithelial cells model during long-term primary in vitro culture.

The oviducts were collected from 45 crossbred gilts, obtained from a commercial herd. The animals were approximately nine months old and displayed at least two regular oestrous cycles. After reaching the anoestrus phase of the cycle, the animals were slaughtered. Collected oviducts were transported to the laboratory and kept in an isolated container at 38°C.

The method used for OEC harvesting and culture has been described in our previous studies [7,8,11]. Briefly, collected oviducts were washed twice in Dulbecco’s phosphate buffered saline (PBS) and sectioned longitudinally. Then, OECs were gently scratched with surgical blades and subsequently digested for 1 h at 37 °C in 1mg/mL solution of collagenase I (Sigma Aldrich, St. Louis, MO, USA) and Dulbecco’s modified Eagle’s medium (DMEM; Sigma Aldrich). The cell suspension was filtered through 40 μm pore size strainer and centrifuged for 10 min at 200× g. After rinsing with PBS, OECs were incubated for another 10 min at 37°C with 0.5% Trypsin/EDTA (Sigma Aldrich), filtered and centrifuged as described above, resuspended in supplemented DMEM (10% FCS, 100 U/mL penicillin, 100 μg/mL streptomycin, and 1 μg/mL amphotericin B) and seeded onto culture dishes. The cells were cultured for up to 30 days (37 °C, 5% CO₂) and the medium of the cultures was changed every three days. Once OECs reached confluence of around 75-80%, they were digested with 0.025% Trypsin/EDTA (Cascade Biologics, Portland, USA), centrifuged and passaged to another culture dish at a seeding density of 2 × 10⁴ cells/cm².

OECs were suspended in TRI Reagent (Sigma Aldrich) after 24 h and 7, 15, and 30 days of culture. The total mRNA isolation was performed using RNeasy MinElute cleanup Kit (Qiagen, Hilden, Germany). The quantity of isolated mRNA was assessed by measuring the absorbance at a wavelength of 260 nm, whereas the purity of RNA samples was determined as 260 nm/280 nm ratio of absorbance using a spectrophotometer (Bioanalyzer 2100, Agilent Technologies, Inc., Santa Clara, CA, USA). Each mRNA sample was diluted to 100 ng/μL concentration.

Total RNA (100 ng) from each pooled sample was subjected to two rounds of sense cDNA amplification (Ambion® WT Expression Kit). The obtained cDNA was used for biotin labeling and fragmentation using Affymetrix GeneChip® WT Terminal Labeling and Hybridization (Affymetrix, Santa Clara, CA, USA). Biotin-labeled fragments of cDNA (5.5 μg) were hybridized to the Affymetrix® Porcine Gene 1.1 ST Array Strip (48°C/20 h). Microarrays were then washed and stained according to the technical protocol, using the Affymetrix GeneAtlas Fluidics Station. The array strips were scanned employing the Imaging Station of the GeneAtlas System. Preliminary analysis of the scanned chips was performed using the Affymetrix GeneAtlas^TM. Operating Software. The quality of gene expression data was confirmed according to the quality control criteria provided by the software. The obtained CEL files were imported into downstream data analysis software.

All of the presented analyses and graphs were generated using Bioconductor and R programming languages. Each CEL file was merged with a description file. To correct background, normalize, and summarize the results, we used the Robust Multiarray Averaging (RMA) algorithm. To determine the statistical significance of the analyzed genes, moderated t-statistics from the empirical Bayes method were performed. The obtained p-value was corrected for multiple comparisons using Benjamini and Hochberg’s false discovery rate. The selection of significantly altered genes was based on a p-value smaller than 0.05 and expression higher than two-fold.

Differentially expressed genes were subjected to selection through examination of genes involved in cell migration regulation. The differentially expressed gene list (separated for up- and down-regulated genes) was uploaded on the DAVID software (Database for Annotation, Visualization and Integrated Discovery) [12], where genes belonging to the terms of all three Gene Ontology (GO) domains were extracted. Expression data of these genes were also subjected to the hierarchical clusterization procedure, with their expression values presented as a heatmap.

Subsequently, we analyzed the relation between the genes belonging to chosen GO terms with the GOplot package [13]. The GoPlot package calculated the z-score: the number of up- regulated genes minus the number of down- regulated genes divided by the square root of the count. This information allowed to estimate the change course of each gene-ontology term.

Interactions between differentially expressed genes/proteins belonging to the studied gene ontology group were investigated using the STRING10 software (Search Tool for the Retrieval of Interacting Genes) [14]. The list of the names of the genes was used as a query for interaction prediction. The search criteria were based on co-occurrences of genes/proteins in scientific texts (text mining), co-expression, and experimentally observed interactions. The results of such analyses generated a gene/protein interaction network where the intensity of the edges reflected the strength of the interaction score.

Finally, the functional interaction between genes that belong to the chosen GO BP terms were investigated with the REACTOME FIViz application to the Cytoscape 3.7.2 software. The ReactomeFIViz app is designed to find pathways and network patterns related to cancer and other types of diseases. This application accesses the pathways stored in the Reactome database, allowing to perform pathway enrichment analysis for a set of genes, to visualize hit pathways using manually laid-out pathway diagrams directly in Cytoscape, and to investigate functional relationships among genes in hit pathways. The application can also access the Reactome Functional Interaction (FI) network, a highly reliable and manually curated pathway-based protein functional interaction network covering over 60% of human proteins.

The research related to animal use has been complied with all the relevant national regulations and instructional policies for the care and use of animals. Bioethical Committee approval no. 246/1992.

Whole transcriptome profiling with Affymetrix microarrays allowed analysis of the gene expression changes between 7, 15, and 30 days of porcine oviductal epithelial cell culture. Using Affymetrix® Porcine Gene 1.1 ST Array Strip, the expression of 12257 transcripts was examined and genes were considered as differentially expressed when presented a fold change higher than abs (2) and corrected p-values lower than 0.05. This set of genes consisted of 2533 different transcripts.

DAVID (Database for Annotation, Visualization and Integrated Discovery) software was used for the extraction of gene ontology biological process terms (GO BP) that contained differentially expressed transcripts. Up- and down-regulated gene sets were subjected to DAVID searches separately, performed only on selected gene sets with adjusted p-value lower than 0.05. The DAVID software analysis showed that the differently expressed genes belonged to 657 Gene ontology terms 25 genes that belong to “coenzyme metabolic process”, “cofactor biosynthetic process” and “cofactor metabolic process” GO BP terms were chosen for further analysis. These sets of genes were subjected to hierarchical clusterization procedure and presented as heat-maps (Fig. 1). Gene symbols, fold changes in expression, Entrez gene IDs and corrected p values of the genes were shown in table 1.

Figure 1

The enrichment of each GO BP term was calculated as z-score and shown on the circle diagram (Fig. 2).

Figure 2

Chosen GO BP terms contained 166 differently expressed genes. Therefore, the calculated mean value of the mean fold change ratio of each gene between 7, 15 and 30 days of culture. Based on that criteria, 10 most downregulated and 10 most up-regulated genes were chosen for further analysis.

In Gene Ontology database genes that formed one particular GO group can also belong to other GO term categories. Accordingly, we explored the gene intersections between selected GO BP terms. The relation between those GO BP terms was presented as a circle plot (Fig. 3) as well as a heatmap (Fig. 4).

Figure 3

Figure 4

The STRING-generated interaction network was created among differentially expressed genes belonging to each selected GO BP terms. Using such prediction method, we obtained a molecular interaction network formed between the proteins derived from the genes of interest (Fig. 5). Finally, we investigated the functional interactions of the chosen genes with REACTOME FIViz app to Cytoscape 3.7.1 software. The results were shown in Fig. 6.

Figure 5

Figure 6

Whole transcriptome profiling approach was employed to analyse the genes belonging to “coenzyme metabolic process”, “cofactor biosynthetic process” and “cofactor metabolic process” GO BP terms. In general, the chosen ontology groups are associated with chemical reactions and cell signalling pathways resulting in the formation of cofactors, as well as with reactions and pathways that involve coenzymes and cofactors.

The most up-regulated gene was ALDH1L2 which encodes a mitochondrial isoform of a similar cytosolic enzyme, ALDH1L1 possessing 10-formyltetrahydrofolate dehydrogenase/hydrolase activities [15]. While the proposed function of the cytosolic enzyme is the overall regulation of the flux of one-carbon groups through folate pool, its mitochondrial counterpart ALDH1L2 may be also involved in producing formate, a pathway shuttling one-carbon groups from mitochondria to cytosol. Additionally, this enzyme might be involved in the regulation of mitochondrial protein biosynthesis [16]. The reaction catalyzed by ALDH1L irreversibly removes carbon groups from the folate pool completing the folate-dependent biosynthetic pathways, which might fulfill several metabolic functions including clearance of excess of folate-bound one-carbon groups [15]. It should be mentioned that, in contrast to cytosolic ALDH1L1, mitochondrial ALDH1L2 is present at detectable levels in several cancer cell lines, denoting different mechanisms of regulation of these enzymes [15].

The next up-regulated gene was P2RX7 which generates a plasma membrane receptor for extracellular ATP. P2RX7 is expressed at a high level by immune and tumor cells; however recent studies have shown its presence in an increasing number of different cells lineages, including stem, neural, bone, endothelial, muscle, and skin cells [17]. P2RX7 was revealed in the oviductal epthelium, and is a cation-selective plasma membrane channel that undergoes, when over-activated, a channel-to-pore transition responsible for a reversible permeabilization of the plasma membrane to low molecular weight aqueous solutes [18]. P2RX7 induces various downstream cell cascades, including inflammatory molecule release, cell proliferation, apoptosis, metabolic events, and phagocytosis [17]. Previously, P2RX7 has been mostly associated with cytotoxicity and proinflammatory cytokines release. However, recent studies revealed that P2RX7 has a potent growth-promoting activity and has been implicated in tumor progression [18].

The upregulation of genes related to coenzyme A (CoA) metabolism, i.e. PANK1 and ACSS2. PANK1 encodes a member of the pantothenate kinase family, which catalyzes the cytosolic phosphorylation of pantothenate (vitamin B5), N-pantothenoylcysteine and pantetheine, and it controls the overall rate of CoA biosynthesis was revealed by the microarray analysis [19]. CoA plays a central role in the metabolism of carboxylic acids, including short- and longchain fatty acids [20,21]. Fatty acids are incorporated into membranes and signaling molecules, and they participate in energy metabolism. These essential functions require activation of fatty acids by acyl-CoA synthetases, e.g. ACSS2, forming an activating thioester linkage between the fatty acid and CoA [22]. The same as PANK1, ACSS2 gene was up-regulated at D7 and D30 during long-term OECs culture.

Additionally, up-regulation of stearoyl-CoA desaturase gene (SCD), which product is an endoplasmic reticulum enzyme that catalyzes the biosynthesis of monounsaturated fatty acids from saturated fatty acids [23] was observed in the array analysis. SCD plays an important role in linking lipid metabolism with adaptive stress signaling and in disorders such as metabolic syndrome, cardiovascular disease and cancer [24].

The next highly up-regulated gene was AASS, which encodes alpha-aminoadipic semialdehyde synthase, a bifunctional protein that catalyzes the first 2 steps in the lysine degradation pathway [25]. The role of AASS in lysine metabolism has been confirmed by earlier studies, where the inactivating mutations in the AASS gene were shown to be a cause of hyperlysinemia [26].

The similar pattern of expression was observed in the case of PDK3 gene. Pyruvate dehydrogenase kinases (PDKs) characterize the pyruvate dehydrogenase complex (PDC) phosphatase inhibitors, which enhances cell glycolysis and facilitates tumor cell proliferation [27,28]. Pyruvate dehydrogenase kinase 3 (PDK3) is a mitochondrial protein considered as a potential pharmacological target for different types of cancer [29].

The previously mentioned genes encode proteins participating in metabolism activation and cell proliferation. During long-term OECs culture, the up-regulation of these genes may be due to the processes related to cell expansion following the initial seeding on culture vessels.

There are several genes belonging to the chosen GO BP terms which were significantly down-regulated after 15 days of in vitro culture. The most down-regulated gene was ACSL4, which plays an essential regulatory role in the ferroptotic cell death process. Acsl4 is responsible for shaping the cellular lipidome and its expression has been reported to correlate with sensitivity to this form of necrotic cell death [30]. Knockdown of ACSL4 has been shown to inhibit erastin-induced ferroptosis in HepG2 and HL60 cell lines, whereas overexpression of ACSL4 by gene transfection restored sensitivity of LNCaP and K562 cells to erastin [31]. When compared to D1 and D7, the down-regulation of ACSL4 in OECs after 15 days in culture may indicate lower sensitivity to ferroptotic necrosis.

The next gene that was down-regulated during the culture period is HAAO. Known as 3-hydroxyanthranilate-3,4-dioxygenase, the product of this gene is an enzyme that catalyzes the biosynthetic pathway from tryptophan to quinolonate. HAAO hypermethylation, which leads to the repression of gene transcription, appears to be a frequent event in endometrial carcinomas [32,33]. Given this, HAAO down-regulation may be related to OECs proliferation in culture.

In conclusion, the current results suggest differential expression of genes belonging to “coenzyme metabolic process”, “cofactor biosynthetic process” and “cofactor metabolic process” GO BP terms in porcine OECs during 30 days of culturing. Considering the biological roles of the most up- and down-regulated genes, it can be concluded that these changes may indicate the increased metabolic and proliferation activity of studied cells in primary in vitro culture. This assumption is consistent with our earlier studies showing that, after long-term culture, the OECs still proliferate and maintain their tube forming properties [10,34].