The effects of 2,3,7,8-tetrachlorodibenzo-p-dioxin (TCDD) on the transcriptome of aryl hydrocarbon receptor (AhR) knock-down porcine granulosa cells

Culture of porcine granulosa cells (AVG-16 cells)

AVG-16 cell line obtained from granulosa cells of pigs was purchased from The European Collection of Authenticated Cell Cultures (06062701; Salisbury, UK; Horisberger, 2006). Previously, we demonstrated that AVG-16 cells are an useful model for studying dioxin effects on ovarian functions (Sadowska et al., 2015). AVG-16 cells were cultured and passaged as previously described (Sadowska et al., 2015; Sadowska et al., 2017). Specifically, one day prior to siRNA transfection, cells were seeded in six-well culture plates with density of 0.7 × 10⁶ cells/three mL Dulbecco’s modified Eagle’s medium (Sigma Aldrich, St. Louis, MO, USA). At ∼70% confluency, the AVG-16 cells were washed (PBS) and the medium was exchanged. The cells were transfected with small interfering RNAs (siRNAs). Untransfected cells (C_UT) were used as control.

AhR gene knock-down in porcine granulosa cells

The granulosa cells were transfected as described previously (Orlowska et al., 2019). Specifically, for the transfection we used Viromer^® BLUE (Lipocalyx GmbH, Halle, Germany) and the mixture of three different siRNAs (anti-AhR 1 + anti-AhR 2 + anti-AhR 3; Sigma Aldrich; Table S1). Negative control cells (C_NEG) were transfected with nontargeted siRNA (Invitrogen, Carlsbad, CA, USA). The transfection mixture was added to the cells in a drop-wise manner and then the cells were cultured for 24 h (37 °C, 5% CO₂, 95% air). Untransfected cells (C_UT) were used as controls. To check the efficacy of the AhR knock-down, the AhR gene expression level and protein abundance were determined in C_UT cells, C_NEG cells and cells transfected with the three relevant siRNA sequences (C_S) by real-time PCR and western blotting, respectively.

TCDD treatment of granulosa cells

In the current study, we compared the transcriptomes of the TCDD-treated AhR knock-down porcine granulosa cells and the untreated AhR knock-down cells. To better recognize the TCDD mode of action, the comparison was made at three time points: 3, 12 and 24 h of cell incubations with TCDD (Sadowska et al., 2017; Nynca et al., 2019). Due to the fact that two examined cell types (TCDD-treated and untreated cells) were both AhR knock-down, the effect of AhR knock-down was excluded. Such an approach was designed to capture only TCDD-induced changes in gene expression. After 24 h of cell transfection with the mixture of siRNAs (anti-AhR 1 + anti-AhR 2 + anti-AhR 3) fresh medium was added and the AhR knock-down cells were treated with 0 nM (control cells) or 100 nM of TCDD (TCDD-treated cells; Sigma Aldrich) for 3, 12 or 24 h (n = two biological replicates per one time point), yielding 12 examined samples. In the present study, we were focused on identifying all possible pathways (molecules) involved in the mechanism of TCDD action in porcine granulosa cells and, therefore, the selected dose of TCDD (100 nM) moderately exceeded its environmentally relevant concentrations. The concentration of 100 nM of TCDD was previously reported to affect porcine granulosa cells steroidogenesis (Gregoraszczuk et al., 2000; Gregoraszczuk, 2002; Jablonska et al., 2014), gene expression (Piasecka-Srader et al., 2016; Sadowska et al., 2017; Nynca et al., 2019), lncRNA expression (Ruszkowska et al., 2018) and proteome (Orlowska et al., 2018) without affecting cell viability and morphology. Moreover, 100 nM was one of the two most frequently used TCDD concentrations in in vitro experiments performed on porcine ovarian cells. After incubation, the medium was removed, cells were washed (PBS) and total RNA was isolated.

Total RNA isolation and evaluation of RNA integrity

Total RNA isolation and evaluation of RNA integrity were performed as described previously (Sadowska et al., 2017; Ruszkowska et al., 2018). The RNA samples with integrity number greater than 8.0 were designated for NGS.

Construction and sequencing of Illumina cDNA libraries

Paired-end cDNA libraries were prepared from purified RNA using Illumina TruSeq Stranded mRNA Sample Preparation Kit (Illumina Inc., San Diego, CA, USA) according to the manufacturer’s protocol. In brief, poly(T) magnetic beads were used to purify mRNA from total RNA sample. The first strand of cDNA was synthesized from the cleaved mRNA fragments using random-sequence primers and reverse transcriptase. In the next step, the RNA template was removed and the second strand of cDNA, with dUTP in place of dTTP, was synthesized by DNA polymerase I. Then, a single “A” nucleotide was added to the 3′ end of double-stranded cDNA, and cDNA was subjected to the adapter ligation. To create final libraries, cDNA samples were amplified through 15 cycles of PCR. The libraries’ profiles were validated using 2100 Bioanalyzer (Agilent Technologies, Santa Clara, CA, USA) to check the size distribution of the libraries. The Qubit High Sensitivity kit (Thermofisher Scientific, Waltham, MA, USA) was used to check the concentration of the samples. Finally, libraries were sequenced on HiSeq2000 instrument (Illumina) using 100 bp paired-end sequencing. The library construction and sequencing were performed by Source BioScience (Nottingham, UK).

Bioinformatic analysis of gene expression

Following sequencing, the obtained cDNA fragments were evaluated and trimmed as previously described (Sadowska et al., 2017). The obtained fragments were mapped to the porcine reference genome with GTF annotation files (Sus Scrofa.Scrofa 11.1.91) using TopHat version 2.0.14 (Trapnell et al., 2012). The alignment files (in BAM format) were used to calculate raw counts per gene (HTSeq package version 0.9.1) (Anders, Pyl & Huber, 2015). The raw counts per gene were normalized according to gene length and library size using DESeq2 version 1.12.3 (Love, Huber & Anders, 2014) package in Bioconductor (Lawrence et al., 2013) and customized algorithms of the authors (R version 3.4.3). A Principal Component Analysis (PCA) and a Pearson’s correlation analysis between the two biological replicates of TCDD-treated and untreated AhR knock-down cells were performed for logarithmic normalized counts with R statistical software (Wickham, 2009). Differentially expressed genes (DEGs) were identified by comparing the expression levels of genes in AhR knock-down cells treated with TCDD to their respective expression levels in AhR knock-down untreated control cells at each time point (3, 12 and 24 h). The corresponding P-values were determined by means of R statistical software using three binominal statistical packages: limma 3.28.6 (Ritchie et al., 2015), edgeR 3.14.0 (Robinson, McCarthy & Smyth, 2010) and DESeq2 for each incubation time. The threshold for the significantly different expression was set at P- adjusted <0.05 and log2 fold change (log2FC) ≥ 1.0 or log2FC ≤ −1.0. Only genes that were classified as DEGs by each of the three binominal statistical packages were used for further analysis.

Functional enrichment analysis

The potential functions of all DEGs were estimated by the Gene Ontology (GO) enrichment analysis (g.Profiler software; P- adjusted <0.05) (Reimand et al., 2007). DEGs were classified into three categories of GO database: “biological processes”, “cellular components” and “molecular functions”.

To establish possible interactions between DEGs, the Search Tool for the Retrieval of Interacting Genes (STRING) S. scrofa database was used. The STRING is a large database of known and predicted protein interactions that cover more than 1,100 organisms (Franceschini et al., 2013) and includes direct (physical) and indirect (functional) associations. Because the aim of the current study was to identify genes associated with the response of AhR knock-down cells to TCDD, DEGs belonging to the GO “response to stimulus” term were selected to be analyzed in detail by STRING database. The gene interaction network created by STRING is built of nodes and edges, where nodes are defined as individual genes, and edges as interactions between the genes. The number of edges indicate the number of interactions for a particular gene in the network. A protein-protein interaction P-value (PPI enrichment P value) was calculated, and the cut-off criterion of the combined score was set to >0.7 (high confidence). In addition, the false discovery rate (FDR) test was applied to determine if all DEGs were enriched.

Real-time PCR

To check the efficacy of the AhR gene knock-down in porcine granulosa cells after the siRNA transfection, real-time PCR was performed. The AhR gene expression level was determined in C_UT, C_NEG and C_S cells after 27 (24 h of transfection + 3 h of cell culture), 36 (24 h of transfection + 12 h of cell culture) and 48 (24 h of transfection + 24 h of cell culture) hours of culture. In addition, real-time PCR was used to validate the results of NGS by analysing the expression of three selected DEGs identified in AhR knock-down granulosa cells treated with TCDD for 3 h. The validation was performed on the same RNA samples which were used for NGS (n= two biological replicates per one time point). The reverse transcription reaction and real-time PCR were performed as previously reported (Sadowska et al., 2017; Ruszkowska et al., 2018). Primers and probes (Thermofisher Scientific, Waltham, MA, USA) for particular genes are presented in Table S2. Glyceraldehyde 3-phosphate dehydrogenase (GAPDH) and β-actin were used as reference genes.

Western blotting

To determine the AhR protein abundance in AhR knock-down cells, the cells were cultured in 25 cm² flasks (seeding density: 3 ×10⁶ cells/12 mL culture medium) and transfected with the siRNAs mixture (anti-AhR 1 + anti-AhR 2 + anti-AhR 3; C_S) and siRNA duplex with an irrelevant sequence (C_NEG) using Viromer^® BLUE as described above. Control cells were cultured without treatments (untransfected cells; C_UT). Twenty four hours after transfection of the cells with the siRNAs mixture medium was exchanged and the cells were incubated with fresh medium for 3, 12 or 24 h (n = four biological replicates per one time point) to match the experimental design of the TCDD study. After incubation, the cells were designated for total protein isolation. The AhR abundance was determined in C_UT, C_NEG and C_S cells after 27 (24 h of transfection + 3 h of cell culture), 36 (24 h of transfection + 12 h of cell culture) and 48 (24 h of transfection + 24 h of cell culture) hours of culture. Protein isolation and western blot analysis were performed as we previously described (Orlowska et al., 2018). Protein samples (15 µg) were separated by SDS-PAGE (9% polyacrylamide gel) and transferred to nitrocellulose membranes (GE Healthcare, Little Chalfont, UK) using the Owl™ VEP-2 Mini Tank Electroblotting System (Thermofisher Scientific) for wet electroblotting (90 mA, 70 min). Following the transfer, the membranes were blocked in TBST containing 5% skimmed milk powder (overnight, 4 °C) and incubated (3 h) with rabbit polyclonal AhR antibodies (ENZO Life Science, New York, NY, USA; diluted in TBST, 1:1,000) or goat polyclonal β-actin antibodies (Santa Cruz Biotechnology, Dallas, TX, USA; diluted in TBST, 1:400). Then the membranes were washed with TBST and incubated (1 h, RT) with HRP-conjugated goat anti-rabbit IgG (Santa Cruz Biotechnology, 1:5,000) or donkey anti-goat secondary antibodies (Santa Cruz Biotechnology; diluted in TBST, 1:10,000). Immobilon chemiluminescent HRP substrate (Millipore, 205 Billerica, MA, USA) was employed to visualize immunolabeled bands. The results of the western blotting were quantified by densitometric scanning of immunoblots (n = four replicates for each protein) with Image Studio Lite (version 5.2).

Statistical analysis

Statistical analysis was performed using Statistica Software (Tulsa, OH, USA). To present the real-time PCR data as arbitrary units of relative expression, gene expression level was normalized to the expression of two reference genes: GAPDH and β-actin. This was done by using comparative cycle threshold (C_T) method and the Quantity Based Active Schematic Estimating (Q-BASE) model (Hellemans et al., 2007). For statistical analysis of the AhR protein abundance, the raw data were normalized using log10 transformation. The densitometric analysis of AhR protein abundance was performed in relation to the reference protein (β-actin). The differences in the AhR gene expression level as well as in the AhR protein abundance between untransfected cells, cells transfected with irrelevant siRNA sequence and cells transfected with the mixture of three anti-AhR siRNAs were evaluated using one-way ANOVA followed by Tukey test. The difference in gene expression level between control and TCDD-treated cells was evaluated using Student’s t-test. The level of significance was set at P < 0.05 for all analyses. All data are presented as mean ± SEM.

The effect of the AhR-targeted siRNAs mixture transfection on the AhR gene expression level and protein abundance in porcine granulosa cells

AhR-targeted siRNA transfection significantly reduced (P < 0.05) the expression level of AhR gene in porcine granulosa cells at all examined time points (3 (Fig. 1A), 12 (Fig. 1B) and 24 h (Fig. 1C). The time points were selected to match the experimental design of the planned TCDD study. The AhR expression level was reduced by 80%, 57% and 72% at 27 (24 h of transfection + 3 h of cell culture), 36 (24 h of transfection + 12 h of cell culture) and 48 (24 h of transfection + 24 h of cell culture) hours of culture, respectively. Nontargeted siRNA transfection (the negative control) had no effect on the AhR expression level (P >  0.05; Figs. 1A–1C).

AhR-targeted siRNA transfection significantly reduced (P < 0.05) the AhR protein abundance in porcine granulosa cells cultured for 3 (Fig. 1D), 12 (Fig. 1E) and 24 h (Fig. 1F). No difference was observed in AhR protein abundance between control untransfected cells (C_UT) and cells transfected with irrelevant siRNA sequence (C_NEG) (P >  0.05). The representative (cropped) blots for AhR and β-actin protein are presented in Figs. 1G–1I. The whole blots can be found in (Fig. S1).

The effects of TCDD on the transcriptome of AhR knock-down porcine granulosa cells

The sequencing data from the study have been submitted to the European Nucleotide Archive database (https://www.ebi.ac.uk/ena) under accession number: PRJEB29985. Sequencing of the AhR knock-down porcine granulosa cell transcriptome provided 779.9 million raw reads, ranging from 51.6 million to 79.5 million per sample. After removing reads containing adapters and low-quality reads (reads length <  90 bp; Phred score Q <  20), the remaining 746.6 million high-quality reads (48.5 to 76.9 million per sample) were mapped to the Ensembl porcine genome with GTF annotation file (Sus Scrofa.Scrofa 11.1.91). The number of reads aligned to the reference genome ranged from 44.5 to 69.5 million per sample, and an average of 96% of these reads were mapped to unique locations. The total number of genes expressed in AhR knock-down granulosa cells of all examined samples ranged from 15,100 to 15,913 (Table S3). The Pearson correlation coefficients (P <0.0001) between the two biological replicates of TCDD-treated and untreated AhR knock-down granulosa cells collected at each time point ranged between 0.971 and 0.999 (Table S4). The results of the PCA revealed that the most marked differences in gene expression pattern between control and TCDD-treated cells were found after 3 h of the treatment. In contrast, the least differences were observed after 12 h. Moreover, genes with expression significantly affected by TCDD after 3-h culture were clustered separately from those of 12- and 24-h cultures (Fig. 2). The MA and Volcano plots present the significant changes (P-adjusted <0.05, log2FC ≥ 1.0 or log2FC ≤ −1.0) in gene expression profiles of AhR knock-down granulosa cells treated with TCDD compared to AhR knock-down control cells (Fig. 3).

A total of 360 DEGs were identified in the current study (Table S5). We found 354 (148 up- and 206 down-regulated) DEGs after 3 h of TCDD treatment and only 10 (1 up- and 9 down-regulated) DEGs after 24 h of the treatment. No DEGs were found after 12 h of the TCDD treatment (Fig. 4). Among all DEGs, four genes (odontogenesis associated phosphoprotein [ODAPH]; ankyrin repeat domain 66 [ANKRD66]; synaptosome associated protein 25 [SNAP25]; gamma-aminobutyric acid type A receptor pi subunit [GABRP]) were identified at two time points, i.e., after 3 and 24 h of TCDD treatment. The log2FC value for DEGs ranged from -2.54 (CXC chemokine receptor type 4; CXCR4) to +2.02 (cytochrome P450, family 26, subfamily A, member 1; CYP26A1). TCDD did not affect AhR as well as CYP1A1 gene expression at any of the examined time points. The expression profile of the top 50 DEGs (i.e., DEGs with the highest log2FC value) is presented in Fig. 5.

Functional enrichment analysis of DEGs

To investigate the possible significance of the identified DEGs in the AhR knock-down porcine granulosal response to TCDD, the genes were classified to three main categories (“biological processes”, “cellular components” and “molecular function”) according to GO database as described in Materials & Methods. Three hundred twenty four out of 360 DEGs were assigned to 172 GO terms (P- adjusted <  0.05); within these GO terms 145 were ascribed to “biological processes”, 21 to “cellular components” and 6 to “molecular function” category (Table S6). The most DEGs of “biological processes” category were allocated to “regulation of biological process” (GO:0050789; 206 DEGs), “regulation of cellular process” (GO:0050794; 199 DEGs), “cellular macromolecule metabolic process” (GO:0044260; 164 DEGs) and “response to stimulus” (GO:0050896; 162 DEGs) (Fig. 6). Most of the genes ascribed to the “cellular component” category were annotated to “nucleus” (GO:0005634; 157 DEGs), “non-membrane-bounded organelle” (GO:0043228; 103 DEGs) and “intracellular non-membrane-bounded organelle” (GO:0043232; 103 DEGs). The most enriched GO term in the “molecular function” category was “binding” (GO:0005488; 227 DEGs) (Fig. 6).

The “response to stimulus” was one of the most enriched GO term containing 162 DEGs. Functional classifications of these genes performed with the use of STRING produced a gene interaction network with 119 nodes and 105 edges (PPI enrichment P - value: 1.0 × 10⁻¹⁶; Fig. 7). The network included genes related to: (i) cell cycle regulation (cell division cycle 45 (CDC45), minichromosome maintenance complex component 2 (MCM2), RAD51 associated protein 1 (RAD51AP1), cyclin-dependent kinase 1 (CDK1), polo like kinase 1 (PLK1), aurora kinase B (AURKB) and DNA topoisomerase 2-alpha (TOPOII)), (ii) interferon signaling pathway (interferon induced protein with tetratricopeptide repeats 2 (IFIT2), 2′-5′-oligoadenylate synthetase 2 (OAS2), interferon regulatory factor 1 (IRF1), interferon regulatory factor 7 (IRF7), inflammatory response protein 6 (RSAD2)) and (iii) cytokine-cytokine receptor interaction (C-C motif chemokine ligand 20 (CCL20), C-C chemokine receptor type 7 precursor (CCR7)). Among the most interacting nodes were CDK1 (15 edges), TOPOII (13 edges), CDC45 (10 edges) and IRG6 (6 edges).

Validation of NGS data by real-time PCR

To validate the NGS results, CYP26A1, CXCR4 and NR4A2 (nuclear receptor subfamily 4, group A, member 2) were chosen for real-time PCR. The first two DEGs were selected based on their log2FC value (top up- and down-regulated genes) as well as their potential importance for granulosa cell functions. The NR4A2 was selected because the gene was also identified in our previous study concerning the TCDD effects on intact porcine granulosa cells (Sadowska et al., 2017). The expression patterns of these three DEGs obtained by real-time PCR were in agreement with NGS results (Fig. 8).