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
Ethical concerns and the lack of healthy tissue sources present limitations for the
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
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% CO2) 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 × 104 cells/cm2.
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 GeneAtlasTM. 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 (
Gene symbols, fold change in expression ratio, Entrez gene IDs, corrected p values and mean value of fold change ratio of studied genes
GENE SYMBOL | RATIO D7/D1 | RATIO D15/D1 | RATIO D30/D1 | ADJUSTED P VALUE D7/D1 | ADJUSTED P VALUE D15/D1 | ADJUSTED P VALUE D30/D1 | ENTREZ GENE ID | MEAN RATIO |
---|---|---|---|---|---|---|---|---|
2,078475 | -1,049840 | -1,389865 | 0,000888 | 0,778847 | 0,027679 | 448980 | -0,120410 | |
2,260462 | 1,012229 | -1,484548 | 0,000495 | 0,948252 | 0,010911 | 100514283 | 0,596048 | |
2,156741 | 1,463697 | 1,431760 | 0,001697 | 0,036412 | 0,040484 | 397239 | 1,684066 | |
2,562460 | 1,513475 | 1,395239 | 0,000490 | 0,019946 | 0,042725 | 100511474 | 1,823725 | |
3,100575 | 1,215495 | 1,317360 | 0,000301 | 0,274540 | 0,110445 | 100512013 | 1,877810 | |
2,424519 | 1,680756 | 1,691567 | 0,000276 | 0,002924 | 0,002153 | 397566 | 1,932281 | |
2,751177 | 1,737387 | 1,421065 | 0,000059 | 0,000779 | 0,006133 | 100157582 | 1,969876 | |
2,691975 | 1,528067 | 1,911722 | 0,000358 | 0,017206 | 0,001723 | 397538 | 2,043921 | |
2,049555 | 1,676912 | 2,729973 | 0,002097 | 0,008605 | 0,000218 | 100233169 | 2,152147 | |
2,109918 | 2,722032 | 2,420194 | 0,000678 | 0,000117 | 0,000152 | 100322873 | 2,417382 | |
3,272247 | 1,857686 | 2,239022 | 0,000043 | 0,000649 | 0,000117 | 100154310 | 2,456318 | |
2,234415 | 2,697009 | 2,514316 | 0,000560 | 0,000143 | 0,000144 | --- | 2,481914 | |
2,978245 | 3,699473 | 2,573760 | 0,000132 | 0,000034 | 0,000122 | 100158056 | 3,083826 | |
4,856476 | 1,981458 | 2,658224 | 0,000100 | 0,004418 | 0,000537 | 100286871 | 3,165386 | |
2,831266 | 3,234674 | 3,436204 | 0,000047 | 0,000015 | 0,000007 | 100155605 | 3,167382 | |
3,411749 | 3,730045 | 3,501388 | 0,000588 | 0,000307 | 0,000279 | 397324 | 3,547727 | |
3,123906 | 3,962945 | 3,669010 | 0,000048 | 0,000012 | 0,000009 | 396592 | 3,585287 | |
2,203769 | 4,264528 | 4,968930 | 0,000233 | 0,000008 | 0,000003 | 100514587 | 3,812409 | |
5,037882 | 1,402767 | 5,796847 | 0,000092 | 0,087483 | 0,000025 | 396670 | 4,079166 | |
PDK3 | 2,773770 | 4,066151 | 5,486628 | 0,000182 | 0,000024 | 0,000006 | 100153858 | 4,108850 |
2,332931 | 4,119242 | 6,132447 | 0,007329 | 0,000476 | 0,000092 | 100513962 | 4,194873 | |
6,349718 | 2,503393 | 5,620930 | 0,000007 | 0,000095 | 0,000002 | 100153866 | 4,824680 | |
7,855562 | 3,491674 | 7,158650 | 0,000007 | 0,000033 | 0,000002 | 100154650 | 6,168629 | |
6,165142 | 8,522042 | 5,795398 | 0,000027 | 0,000007 | 0,000010 | 497623 | 6,827528 | |
4,535793 | 8,626347 | 22,102755 | 0,000038 | 0,000004 | 0,000001 | 100151976 | 11,754965 |
The enrichment of each GO BP term was calculated as z-score and shown on the circle diagram (F
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 (
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 (
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
The next up-regulated gene was
The upregulation of genes related to coenzyme A (CoA) metabolism, i.e.
Additionally, up-regulation of stearoyl-CoA desaturase gene (
The next highly up-regulated gene was
The similar pattern of expression was observed in the case of
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
The next gene that was down-regulated during the culture period is
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