- Open Access
Identification and characterization of long non-coding RNAs in porcine granulosa cells exposed to 2,3,7,8-tetrachlorodibenzo-p-dioxin
Journal of Animal Science and Biotechnology volume 9, Article number: 72 (2018)
Long non-coding RNAs (lncRNAs) may regulate gene expression in numerous biological processes including cellular response to xenobiotics. The exposure of living organisms to 2,3,7,8-tetrachlorodibenzo-p-dioxin (TCDD), a persistent environmental contaminant, results in reproductive defects in many species including pigs. The aims of the study were to identify and characterize lncRNAs in porcine granulosa cells as well as to examine the effects of TCDD on the lncRNA expression profile in the cells.
One thousand six hundred sixty-six lncRNAs were identified and characterized in porcine granulosa cells. The identified lncRNAs were found to be shorter than mRNAs. In addition, the number of exons was lower in lncRNAs than in mRNAs and their exons were longer. TCDD affected the expression of 22 lncRNAs (differentially expressed lncRNAs [DELs]; log2 fold change ≥ 1, P-adjusted < 0.05) in the examined cells. Potential functions of DELs were indirectly predicted via searching their target cis- and trans-regulated protein-coding genes. The co-expression analysis revealed that DELs may influence the expression of numerous genes, including those involved in cellular response to xenobiotics, dioxin metabolism, endoplasmic reticulum stress and cell proliferation. Aryl hydrocarbon receptor (AhR) and cytochrome P450 1A1 (CYP1A1) were found among the trans-regulated genes.
These findings indicate that the identified lncRNAs may constitute a part of the regulatory mechanism of TCDD action in granulosa cells. To our knowledge, this is the first study describing lncRNAs in porcine granulosa cells as well as TCDD effects on the lncRNA expression profile. These results may trigger new research directions leading to better understanding of molecular processes induced by xenobiotics in the ovary.
Recent advances in RNA sequencing (RNA-Seq) technologies led to the discovery of tens of thousands non-coding RNAs (ncRNAs) across the entire transcriptomes of many species. Previously, ncRNAs were considered as evolutionary junk or transcriptional noise, but recently they were reported to play a crucial role in the regulation of gene expression [1, 2]. The class of ncRNAs includes short ncRNAs (e.g., microRNA, piwiRNA), housekeeping ncRNAs (e.g., ribosomal RNA, transfer RNA) as well as long non-coding RNAs (lncRNAs), which constitute the largest group of ncRNAs [3, 4]. In broad terms, lncRNAs are defined as non-coding transcripts longer than 200 nucleotides (nt). Similar to mRNAs, most known lncRNAs: 1) are transcribed by RNA polymerase II, 2) contain an intron-exon structure and 3) undergo post-transcriptional modifications, e.g., 5′ capping, 3′ polyadenylation or alternative splicing . lncRNAs originate from intergenic regions of the genome (lincRNAs), introns of protein-coding genes (linRNAs) or the opposite strand of DNA (antisense lncRNAs) [5, 6]. It was demonstrated that lncRNA sequences are not highly conserved and their expression profiles differ among species, developmental stages and cell types. Due to their variable nature, lncRNAs are still poorly characterized [7, 8].
It was found that lncRNAs regulate gene expression via numerous mechanisms. lncRNAs may be implicated in e.g., chromatin organization, epigenetic modification, transcriptional regulation, RNA processing and modification, regulation of miRNA activity as well as protein stability or cellular protein localization [3, 8, 9]. lncRNAs seem to be involved in many fundamental biological processes including maintenance and induction of stem cell pluripotency, cell differentiation, cell proliferation and apoptosis [6, 8]. The mutations of lncRNAs and deregulation of lncRNA expression have been associated with a number of human diseases, such as different types of cancer as well as neurobiological, cardiovascular and autoimmune disorders [6, 8]. The results of recent studies also indicate that lncRNAs take an active part in cellular response to various xenobiotics .
2,3,7,8-tetrachlorodibenzo-p-dioxin (TCDD), a chlorinated xenobiotic and a member of the polycyclic aromatic hydrocarbon (PAH) family, is a by-product of human activity including herbicide production, waste incineration and fossil fuel burning . Because of TCDD lipid solubility and high resistance to degradation, its half-life is long and ranges from 7 to 10 years in human bodies and from 25 to 100 years in the environment. As a result, TCDD persists in soil and water sediments as well as in plant and animal organisms [11, 12]. The main intracellular mechanism of TCDD action involves the activation of the aryl hydrocarbon receptor (AhR) pathway. AhR is a ligand-activated transcription factor and a member of the basic helix-loop-helix/PER-ARNT-SIM (bHLH-PAS) family. Upon activation, the ligand-AhR complex translocates to the nucleus, dimerizes with AhR nuclear translocator (ARNT) and induces the expression of different genes including xenobiotic-metabolizing enzymes . Exposure to TCDD may result in a variety of harmful short- and long-term effects, such as wasting syndrome, cancer and neurological dysfunctions. TCDD has also been demonstrated to cause endocrine disruption and reproductive defects in many species including pigs [14,15,16,17,18,19].
Granulosa cells play a critical role in maintaining ovarian function and, in consequence, female fertility. These cells protect and nurture oocytes as well as produce steroid hormones (estradiol [E2], progesterone [P4]) responsible for creating an optimal environment for follicular development, fertilization, implantation and embryo growth . Due to the AhR presence in porcine granulosa cells [21, 22] and the fact that TCDD affects granulosal production of E2 and P4 in pigs [14, 15, 19], it is necessary to indicate molecular targets of TCDD in these cells. The effects of TCDD on gene expression in porcine granulosa cells were described in our recent paper , but information concerning lncRNAs in the porcine ovary is scarce . We assumed that the transcriptome of porcine granulosa cells contains lncRNAs and that some lncRNAs are involved in TCDD action in the cells. Therefore, the aims of the current study were to: 1) identify and characterize porcine granulosa lncRNAs, 2) examine the effects of TCDD on the lncRNA expression profile (i.e., to identify differentially expressed lncRNAs [DELs]) and 3) predict target genes for DELs in order to tentatively evaluate their regulatory role in granulosal response to TCDD.
Culture of AVG-16 cells
AVG-16 cell line originating from granulosa cells of medium porcine follicles was used in the current study (06062701; The European Collection of Authenticated Cell Cultures; UK; ). Previously, we have showed that AVG-16 cells are morphologically and physiologically similar to primary porcine granulosa cells and are a good model for studying xenoestrogen effects on ovarian functions . Before the experiment, the cells were thawed and cultured in six-well plates with seeding density of 1 × 106 cells/3 mL culture medium (Dulbecco’s modified Eagle’s medium [DMEM] with 2 mmol/L L-glutamine, 10% fetal bovine serum [FBS], 0.1 mmol/L non-essential amino acid [NEAA], 2.5 ng/mL fibroblast growth factor-basic human [bFGF] and antibiotic mixture [100 U penicillin, 100 μg streptomycin and 0.25 μg amphotericin B/mL]; Sigma Aldrich, St. Louis, MO, USA) [22, 23]. After reaching 60–65% confluency, the cells were treated with TCDD (100 nmol/L; Sigma Aldrich) for 3, 12 or 24 h (n = 2 biological replicates per one time point). To bring out the potential of TCDD to transduce intracellular signaling in the examined cells, the selected dose of TCDD moderately exceeded its environmentally relevant concentration. Moreover, 100 nmol/L of TCDD was found to be effective in porcine granulosa cells (steroidogenesis: [26,27,28]; gene expression: [23, 29]). At the end of the culture, medium was removed, cells were washed twice with a phosphate-buffered saline (PBS; Sigma Aldrich) and designated for total RNA isolation.
Total RNA isolation, evaluation of RNA integrity, cDNA library construction and sequencing
Total RNA was isolated from cells using peqGold TriFast (Peqlab Biotechnologie GmbH, Erlangen, Germany) according to manufacturer’s instructions. RNA concentration and quality were determined spectrophotometrically (NanoVue Plus, GE Healthcare, Little Chalfont, UK). The integrity of total RNA was assessed by Agilent 2100 Bioanalyzer using RNA 6000 Nano LabChip Kit (Agilent Technologies, Santa Clara, CA, USA). Only samples with RNA integrity number (RIN) values higher than 8.0 were used for RNA-Seq.
Depleted RNA obtained from 400 ng of total RNA was used to construct cDNA libraries (TruSeq RNA Sample Preparation Kit; Illumina, San Diego, CA, USA). Following RNA purification and fragmentation, first and second cDNA strands were synthesized. Next steps included 3’ end adenylation, adapter ligation and library amplification (PCR). Quantification of the cDNA library templates was performed using KAPA Library Quantification Kit (Kapa Biosystem, Wilmington, MA, USA). Library profiles were estimated with the use of DNA High Sensitivity LabChip Kit (Agilent Technologies) and 2100 Bioanalyzer. The cDNA library templates were sequenced using HiSeq 2500 high throughput sequencing instrument (Illumina) with 100 paired-end sequencing. The RNA-Seq analyses of protein-coding genes were described in a separate manuscript . This article is focused on the identification and characterization of long non-coding RNAs in porcine granulosa cells exposed to TCDD.
Sequencing data analysis and transcriptome assembly
The sequencing data from this study have been submitted to the NCBI BioProject database (http://www.ncbi.nlm.nih.gov/bioproject) under accession number: PRJNA429720. The quality of cDNA fragments obtained after sequencing (raw reads) was first evaluated using FastQC program (http://www.bioinformatics.babraham.ac.uk/projects/fastqc/). Next, the raw reads were trimmed (Trimmomatic tool version 0.32) by removing the adapter sequences and reads shorter than 90 nt . After trimming to the same length (90 nt), the fragments were mapped to the porcine reference genome (Sus_scrofa.Sscrofa10.2; Ensembl database) using STAR version 2.4 . The mapped reads were assembled into contigs with Cufflinks version 2.2.1  and StringTie version 1.0.4 [33, 34].
Identification of lncRNAs
A customized multi-step pipeline (Fig. 1a) was employed to identify putative novel lncRNAs in porcine granulosa cells. The detailed identification steps included removal of: 1) protein-coding transcripts and transcripts that overlap, in sense orientation, with at least one base of all porcine protein-coding sequences annotated in the Ensembl database; 2) transcripts with a single exon [35, 36]; 3) transcripts shorter than 200 nt; 4) transcripts which were predicted by Coding Potential Calculator (CPC) (score < 0) , Coding-Non-Coding Index (CNCI) (score < 0) , FlExible Extraction of LncRNAs (FEELnc) (coding potential > 0.558) , Pfam Scan (v1.3) (E-value < 0.001)  and PLEK (score < 0)  to encode proteins and 5) transcripts blasted (Blast2GO) to small RNAs (rRNA, tRNA, snRNA, snoRNA, miRNA, etc.) annotated in Rfam database. The remaining transcripts were considered as lncRNAs and were divided into known lncRNAs (based on the pig Ensembl database lncRNAs annotation) and novel lncRNAs. In addition, the obtained lncRNAs were classified depending on their genomic locations. The protein-coding sequences removed in step 1 of the identification pipeline were characterized and their genomic features were compared with those of lncRNAs (Welch’s t-test in R package, P < 0.05).
Differential expression analysis of lncRNAs
The combined Cufflinks and StringTie tools allowed us to analyze the expression profile of both known and unannotated lncRNAs. The identified lncRNA sequences were normalized to FPKM (fragments per kilobase of exon per million fragments mapped) values using the Cuffnorm (version 2.2.1) in the Cufflinks package . To identify differentially expressed lncRNAs (DELs), the expression levels of lncRNAs in TCDD-treated cells were compared to their respective expression levels in control cells. The corresponding P-values were determined for each incubation time by means of Cuffdiff . P-adjusted < 0.05 and log2 fold change (log2FC) ≥ 1.0 were set as a threshold for the significantly different expression.
Cis- and trans-regulated target gene prediction and Gene Ontology (GO) enrichment analysis
Cis-acting lncRNAs regulate the expression of genes that are positioned in the vicinity of their transcription sites, whereas trans-acting lncRNAs modulate the expression of genes being at independent loci . Protein-coding genes located within 10 kb from DELs were screened with the use of R package  and selected as DEL potential cis-regulated targets [44, 45]. To explore the trans-type interaction, the Pearson’s correlations coefficients (r > 0.7 or r < − 0.7) between DELs and protein-coding genes were analyzed. To analyze functions of the potential DEL target genes and their involvement in biological processes, the GO database was used (the established criteria: P-adjusted < 0.05). GO enrichment analysis was performed using g:Profiler software .
Alternative splicing events analysis
To recognize the differential expression of exons and the occurrence of splice junctions (i.e., differential usage of exons and junctions) in the identified lncRNAs, QoRTs  and JunctionSeq  tools were applied. First, QoRTs was used to generate the raw exon- and splice junction-counts. Then, JunctionSeq was employed to visualize the identified sites of differential usage of exons as well as splice junction in lncRNA loci of TCDD-treated cells in comparison to control cells (P-adjusted < 0.05).
To validate RNA-Seq results, three lncRNAs, expression of which was affected by TCDD in two (TCONS_00016901; TCONS_00031035) or three (TCONS_00034713) incubation times, were selected for real-time PCR (qRT-PCR). Complementary DNA was synthesized using RNA isolated from granulosa cells (control and TCDD-treated), 0.5 μmol/L oligo(dT)15 primer (Roche, Basel, Switzerland), 1 μmol/L hexanucleotide primers, 10 U RNase Out (Sigma Aldrich) and Omniscript RT Kit (Qiagen, Hilden, Germany). Reverse transcription reaction was performed at 37 °C for 1 h (Veriti Thermal Cycler, Thermofisher Scientific, Waltham, MA, USA). Specific primers for particular lncRNA were synthesized by Thermofisher Scientific Company (TCONS_00031035 – Assay ID: AIPAFN7; TCONS_00016901 – Assay ID: AIN1HH2; TCONS_00034713 – Assay ID: AIFAT9T). Glyceraldehyde 3-phosphate dehydrogenase (GAPDH; Assay ID: Ss03373286_u1) and β-actin (Assay ID: Ss03376563_uH) were used as reference genes. Real-time PCR was performed using TaqMan® Universal PCR Master Mix and TaqMan Gene Expression Assay (Thermofisher Scientific) in Applied Biosystems 7500 Fast Real-Time PCR System (Thermofisher Scientific). Thermal cycling conditions consisted of an initial denaturation step at 95 °C for 10 min and then 40 cycles of denaturation at 95 °C for 15 s followed by primer annealing at 60 °C for 1 min. The qRT-PCR for each sample was carried out in duplicate and non-template control was included in each run. To present the data as arbitrary units of the relative expression, lncRNA expression levels were normalized to the expression of GAPDH and β-actin. This was done by using comparative cycle threshold (Ct) method and the quantity based active schematic estimating (Q-BASE) model . Data were expressed as mean ± SEM. The difference in lncRNA expression level between control and TCDD-treated cells was evaluated using Student’s t-test (Statistica Software Inc., Tulsa, OH, USA). Differences with a probability of P < 0.05 were considered significant.
The sequencing of the porcine granulosa cell transcriptome provided 528.5 million reads (50–100 nt). After discarding adaptor sequences and low-quality sequences (Phred score Q < 20), the remaining 476.2 million reads were mapped to the annotated whole porcine genome (Sus_scrofa.Sscrofa10.2). The percentage of unique mapped reads ranged from 76.8% to 80.1%. As a result, 56,477 transcripts were collected (see Additional file 1).
Identification of lncRNAs in the porcine granulosa cell transcriptome
To identify lncRNAs from 56,477 mapped transcripts, the customized multi-step identification pipeline was applied (Fig. 1a). Twelve thousand one hundred seventy-three transcripts were obtained after filtering the protein-coding transcripts and transcripts that overlap with annotated protein-coding sequences. The removal of single exon transcripts and transcripts shorter than 200 nt yielded 4,160 transcripts and the evaluation of the protein-coding potential (Fig. 1b) produced 1,689 potentially non-coding long transcripts. The stringency of the selection process was additionally increased by discarding transcripts which blasted to small RNAs annotated in the Rfam database. Eventually, 1,666 RNA sequences were identified as lncRNAs, and 42 of them have already been annotated in databases. According to genomic localization of lncRNAs, the identified porcine novel lncRNAs were classified as intergenic lncRNAs (1,319 transcripts) and antisense lncRNAs (305 transcripts). The group of antisense lncRNAs (Fig. 1c) were located at the opposite strand of annotated protein-coding genes (281 lncRNAs), snoRNAs (12 lncRNAs), snRNAs (7 lncRNAs), miRNAs (3 lncRNAs), processed transcript (1 lncRNA) and pseudogene (1 lncRNA).
The comparison of lncRNAs and mRNAs of the granulosa cell transcriptome
The transcript length, exon length, exon number and expression level were compared between the identified lncRNAs (1,666 transcripts) and mRNAs (42,913 transcripts). The length of most of the lncRNAs ranged from 200 to 1,000 nt (Fig. 2a), while the length of most lncRNA exons ranged from 50 to 200 nt (Fig. 2b). Moreover, a majority of lncRNAs consisted of two exons (Fig. 2c) and the expression level of more than 50% of lncRNAs was lower than FPKM value of 2 (Fig. 2d).
The average length of lncRNAs (832.6 ± 33.06 nt) was shorter (P < 0.05) than that of mRNA (3284.4 ± 12.49 nt) (Fig. 3a) and the average exon length of lncRNA (318.2 ± 10.39 nt) was longer (P < 0.05) than that of mRNA (295.3 ± 1.03 nt) (Fig. 3b). In addition, the mean exon number of lncRNAs (2.6 ± 0.02) was lower (P < 0.05) than the mean exon number of mRNAs (11.09 ± 0.05) (Fig. 3c). The mean expression level (0.51 ± 0.003 vs. 0.49 ± 0.001 FPKM) did not differ (P > 0.05) between lncRNA and mRNA (Fig. 3d).
Differentially expressed lncRNAs in TCDD-treated porcine granulosa cells
Twenty two differentially expressed lncRNAs (P-adjusted < 0.05 and log2FC ≥ 1.0) were identified. The expression of 15, 4 and 7 lncRNAs were significantly altered after 3, 12 and 24 h of TCDD treatment, respectively. The log2FC value for DELs ranged from 3.77 (TCONS_00041714) to − 1.01 (TCONS_00048132) (Table 1). Among all DELs: 1) one lncRNA (TCONS_00034713) was up-regulated by TCDD in all incubation times, 2) two lncRNAs (TCONS_00016901 and TCONS_00031035) were up-regulated by the dioxin after two incubation times (3, 24 and 12, 24 h of TCDD treatment, respectively), 3) two lncRNAs (TCONS_00048132; TCONS_00030731) were down-regulated, and the downregulation was found after 24 h of TCDD treatment, 4) five lncRNAs (TCONS_00005658; TCONS_00016901; TCONS_00048979; TCONS_00060223; TCONS_00064401) were expressed only in TCDD-treated cells and 5) three lncRNAs (TCONS_00038918, TCONS_00030731, TCONS_00064964) up-regulated by TCDD were located in the antisense strand of protein-coding transcripts (MAF bZIP transcription factor F [MafF], transforming growth factor beta-induced [TGFβI], NHS actin remodeling regulator [NHS]). The expression profile of up- and down-regulated lncRNAs identified in granulosa cells after 3, 12 or 24 h of TCDD treatment is presented in Fig. 4.
The cis- and trans-target genes for DELs
The potential target genes for DELs acting in a cis- or trans-regulatory manner were predicted to investigate the possible significance of these lncRNAs in the porcine granulosal response to TCDD. In silico analysis produced 10 cis-target genes for DELs (see Additional file 2) and none of them were enriched (P > 0.05) in the GO classification.
To analyze the trans-type interaction between DELs (Table 1) and their target genes, the co-expression analysis was performed. As a result, a total of 916 negative (see Additional file 3) and 1,143 positive (see Additional file 4) correlations were detected. The negatively co-expressed genes were enriched in 28 GO terms (26 under “biological process” and two under “molecular function”) including “intracellular signal transduction” (GO:0035556), “response to stimulus” (GO:0050896) and “negative regulation of biological process” (GO:0048519). The negatively co-expressed trans-target genes involved, among others, AhR, TGFβI, endoplasmic reticulum to nucleus signaling 1 (ERN1), LIM domain kinase 2 (LIMK 2) and ephrin type-B receptor 4 (EPHB4) (see Additional file 5). The positively co-expressed genes were enriched in five GO terms (four in “biological process” and one in “molecular function”). Some of the genes were related to cellular response to xenobiotics, including “cellular response to chemical stimulus” (GO:0070887) and “regulation of signal transduction” (GO:0009966). This group of trans-target genes involved cytochrome P4501A1 (CYP1A1), MafF and NHS (see Additional file 6).
Alternative splicing events of lncRNAs in TCDD-treated porcine granulosa cells
We identified 33 events (sites) of differential usage of exons and splice junctions in lncRNAs loci of porcine granulosa cells after 3 h of TCDD treatment and only four events after 12 h of the treatment. No differential usage of exons and junctions was found after 24 h of TCDD treatment. Among the identified events, 28 and 2 were defined as differentially expressed exons, whereas 5 and 2 were described as splice junctions after 3 h and 12 h of TCDD treatment, respectively (see Additional file 7). The JunctionSet profile plot for the exemplary lncRNA XLOC (XLOC_020835) is presented in Fig. 5. Two lncRNAs - TCONS_00040606 and TCONS_00040607 – may be formed on the basis of this selected XLOC. Moreover, the expression of lncRNA TCONS_00040607 correlated negatively with the expression of AhR.
Validation of RNA-Seq data by real-time PCR
To validate the RNA-Seq results, three up-regulated lncRNAs, i.e., TCONS_00016901, TCONS_00031035 and TCONS_00034713 were chosen for real-time PCR. The expression of the selected lncRNAs entirely confirmed the results obtained by RNA-Seq (Fig. 6).
Recent advances in sequencing technologies revealed that the genome is widely transcribed, yielding a large number of ncRNAs. In the current study, we investigated the genome-wide lncRNA expression profile of porcine granulosa cells exposed to TCDD. The effects of TCDD on the porcine granulosa cell transcriptome have been previously examined, however the studies have been focused solely on the expression of genes (RNA-Seq ; real-time PCR ). In comparison to humans and rodents, information concerning the identification and expression of lncRNAs in pigs is limited [24, 50]. To the best of our knowledge, lncRNAs have not yet been identified in porcine granulosa cells and TCDD effects on the lncRNA expression profile have also not been examined.
A total of 1,666 lncRNAs and 42,913 mRNAs were identified in porcine granulosa cells in the present study. The majority of the identified lncRNAs were classified as intergenic lncRNAs. The identified lncRNAs were found to be shorter in comparison to those of mRNAs. In addition, the exon number of lncRNAs was lower and the exon length was longer than those of mRNAs. The features of the identified lncRNAs were similar to those found in other studies [24, 43, 51, 52].
The analysis of the expression profile of lncRNAs in TCDD-treated porcine granulosa cells revealed the presence of 22 DELs. Similarly to previously reported TCDD-induced changes in granulosal gene expression , the most of DELs (15 out of 22) were detected after 3 h of TCDD treatment. Since TCDD was shown to induce oxidative stress in the ovary [53, 54], the altered expression of lncRNAs after 3 h of TCDD treatment may be associated with the regulation of the early cellular stress response to TCDD.
Potential functions of DELs identified in porcine granulosa cells were indirectly predicted via searching their target cis- and trans-regulated protein-coding genes. Screening of genes which were located within 10 kb from DELs enabled us to select 10 cis-target genes, but none of them was enriched into GO classification. To identify the trans-target genes presumptively regulated by lncRNAs, we correlated the expression levels of DELs and protein-coding genes. A total of 916 negative and 1,143 positive correlations were detected. The GO results of negatively co-expressed genes suggested a possible involvement of DELs in a variety of biological processes such as “intracellular signal transduction”, “response to stimulus” and “negative regulation of biological process”. Interestingly, AhR was found among these negatively trans-regulated genes. Although the downregulation of AhR is a typical response to TCDD, recent reports concerning ncRNAs may provide more details on TCDD mechanism of action. Li et al.  demonstrated that miR-203 (microRNA) suppressed the expression of AhR in TCDD-treated human lung and hepatic cells. In the current study, we identified TCONS_00040607 i.e., lncRNA, the expression of which correlated negatively with that of AhR after 3 h of TCDD treatment. It is possible that TCONS_00040607 affects the expression of AhR and is involved in its negative regulation during the cellular response to TCDD. These results suggest that the regulation of AhR by lncRNAs may constitute part of cellular defense mechanism against dioxins. This hypothesis, however, needs additional experimental verification.
Some of the positively co-expressed genes were enriched in GO terms associated with “cellular response to chemical stimulus” and “regulation of signal transduction”. We found that the expression of CYP1A1, coding a protein playing an important role in the metabolism of xenobiotics , correlated positively with the expression of two DELs (TCONS_00034713 and TCONS_00031305). Due to the fact that TCDD treatment usually induces CYP1A1 expression, this enzyme is considered to be a molecular marker of TCDD action [57, 58]. In our previous study, the expression of CYP1A1 increased significantly in a time-dependent manner after 3, 12 and 24 h of porcine granulosa cell incubation with TCDD . We also demonstrated that the TCDD binding to the porcine CYP1A1 active site resulted in a rapid closure of the enzyme substrate channels. This phenomenon may partially explain TCDD’s high resistance to biodegradation . If the TCDD binding causes a continuous CYP1A1 blockage, the cellular response to TCDD may induce an extended synthesis of CYP1A1. The fact that the expression of two DELs: TCONS_00034713 and TCONS_00031305 correlated positively with the expression of CYP1A1 indicates their supportive role in the cellular reaction to TCDD.
Five DELs (TCONS_00005658; TCONS_00016901; TCONS_00048979; TCONS_00060223; TCONS_00064401) were found to be expressed only in TCDD-treated cells. The expression of two DELs (TCONS_00048979 and TCONS_00060223) correlated negatively with the expression of the same three genes: ERN1, LIMK2 and EPHB4. ERN1 is associated with endoplasmic reticulum stress and ER protein folding [60, 61]. EPHB4, in turn, is linked with the proliferation of ovarian carcinoma cells  and LIMK2 with actin microfilament disruption in porcine oocytes . The obtained data imply that TCONS_00048979 and TCONS_00060223 are mediators of TCDD action in porcine granulosa cells.
In addition, we identified three DELs located in the antisense strand of protein coding genes. These antisense strands were found to positively or negatively regulate the expression of their sense counterparts and, therefore, they attracted a lot of attention [64, 65]. In the current study, TCONS_00038918, TCONS_00030731 and TCONS_00064964 were found to be located in the respective antisense strands of MafF, TGFβI and NHS. MafF belongs to the small MAF family of basic-leucin zipper transcription factors, which affect genes encoding proteins responsible for xenobiotic metabolism and antioxidation [66,67,68]. TGFβ1 encodes an extracellular matrix (ECM) protein reported to interact with various matrix molecules (collagens, fibronectin and laminin), contributing to cell adhesion, proliferation and migration . ECM proteins were found to be affected by TCDD in marmosets . NHS products, in turn, were described to maintain cell morphology by remodeling the actin cytoskeleton . The in silico data concerning the possible relationships between lncRNAs and protein-coding genes in porcine granulosa cells treated with TCDD are supported by the results of our recent study . In this study, the abundance of heat shock proteins as well as cytoskeleton and ECM proteins were significantly affected by TCDD in porcine granulosa cells.
Similarly to protein-coding genes, lncRNAs also undergo alternative splicing, resulting in the formation of numerous lncRNA isoforms and, in consequence, extending their regulatory capabilities [73, 74]. It was demonstrated that xenobiotics may affect alternative splicing, modifying the process in a specific manner . In the current study, the events of differential usage of exons and splice junctions in lncRNAs loci were identified after 3 h and 12 h of cell incubation with TCDD. The majority of the identified events were defined as differentially expressed exons and only a few events were described as splice junctions. It is of interest that two alternatively spliced forms for XLOC_020835 were detected after TCDD treatment. In the current study, the expression of one of these forms i.e., TCONS_00040607, was found to negatively correlate with the expression of AhR. These facts implicate that TCONS_00040607 is involved in the regulation of the TCDD-affected AhR expression in porcine granulosa cells.
In the current study, we identified and characterized lncRNAs of porcine granulosa cells. We also examined the effects of TCDD on the lncRNA expression profile. The co-expression analysis revealed that the identified lncRNAs may influence the expression of numerous genes, including those involved in: 1) cellular response to TCDD (e.g., AhR), 2) dioxin metabolism (e.g., CYP1A1, MafF), 3) endoplasmic reticulum stress (e.g., ERN1) as well as 4) adhesion and proliferation of cells (TGFβ1, NHS, LIMK2). Our results suggest that lncRNAs constitute a part of the regulatory apparatus of TCDD action in porcine granulosa cells. This study may provide the foundation for future research focused on molecular effects exerted by TCDD in ovarian cells.
Aryl hydrocarbon receptor
- antisense lncRNAs:
LncRNAs originate from the opposite strand of DNA
AhR nuclear translocator
Fibroblast growth factor-basic human
The basic helix-loop-helix/PER-ARNT-SIM
Coding Potential Calculator
Cytochrome P450 1A1
Differentially expressed lncRNAs
Dulbecco’s modified Eagle’s medium
- E2 :
Ephrin type-B receptor 4
Endoplasmic reticulum to nucleus signaling 1
Fetal bovine serum
FlExible Extraction of LncRNAs
Gene Ontology database
- LIMK 2:
LIM domain kinase 2
LncRNAs originate from intergenic regions of the genome
LncRNAs originate from introns of protein-coding genes
Long non-coding RNAs
MAF bZIP transcription factor F
MEM non-essential amino acid solution
NHS actin remodeling regulator
- P4 :
Polycyclic aromatic hydrocarbon
RNA integrity number
Transforming growth factor beta-induced
Dempsey JL, Cui JY. Long non-coding RNAs: a novel paradigm for toxicology. Toxicol Sci. 2017;155:3–21.
Yu L, Tai L, Zhang L, Chu Y, Li Y, Zhou L. Comparative analyses of long non-coding RNA in lean and obese pigs. Oncotarget. 2017;8:41440–50.
Dykes IM, Emanueli C. Transcriptional and post-transcriptional gene regulation by long non-coding RNA. Genomics Proteomics Bioinformatics. 2017;5:177–86.
Liu D, Mewalal R, Hu R, Tuskan GA, Yang X. New technologies accelerate the exploration of non-coding RNAs in horticultural plants. Hortic Res. 2017;4:17031.
Prabhakar B, Zhong XB, Rasmussen TP. Exploiting long noncoding RNAs as pharmacological targets to modulate epigenetic diseases. Yale J Biol Med. 2017;90:73–86.
Chen X, Yan CC, Zhang X, You ZH. Long non-coding RNAs and complex diseases: from experimental results to computational models. Brief Bioinform. 2016;18:558–76.
Perry RB, Ulitsky I. The functions of long noncoding RNAs in development and stem cells. Development. 2016;143:3882–94.
Szcześniak MW, Makałowska I. lncRNA-RNA interactions across the human Transcriptome. PLoS One. 2016;11(3):e0150353.
Guo Q, Cheng Y, Liang T, He Y, Ren C, Sun L, et al. Comprehensive analysis of lncRNA-mRNA co-expression patterns identifies immune-associated lncRNA biomarkers in ovarian cancer malignant progression. Sci Rep. 2015;5:17683.
Basham KJ, Leonard CJ, Kieffer C, Shelton DN, McDowell ME, Bhonde VR, et al. Dioxin exposure blocks lactation through a direct effect on mammary epithelial cells mediated by the aryl hydrocarbon receptor repressor. Toxicol Sci. 2015;14:36–45.
ATSDR. Draft Update Toxicological Profile for Chlorinated Dibenzo-p-dioxins.Prepared by Research Triangle Institute for U.S. Department of Health and HumanServices, Agency for Toxic Substances Disease Registry (ATSDR), Atlanta, GA. 677 pp +appendices (1998).
Nicolopoulou-Stamati P, Pitsos MA. The impact of endocrine disruptors on the female reproductive system. Hum Reprod. 2001;7:323–30.
Pocar P, Fischer B, Klonisch T, Hombach-Klonisch S. Molecular interactions of the aryl hydrocarbon receptor and its biological and toxicological relevance for reproduction. Reproduction. 2005;129:379–89.
Piekło R, Grochowalski A, Gregoraszczuk EL. 2,3,7,8-tetrachlorodibenzo-pdioxin alters follicular steroidogenesis in time- and cell-specific manner. Exp Clin Endocr Diab. 2000;108:299–304.
Grochowalski A, Chrzaszcz R, Piekło R, Gregoraszczuk EL. Estrogenic and antiestrogenic effect of in vitro treatment of follicular cells with 2,3,7,8- tetrachlorodibenzo-p-dioxin. Chemosphere. 2001;43:823–7.
Petroff BK, Roby KF, Gao X, Son D, Williams S, Johnson D, et al. A review of mechanisms controlling ovulation with implications for the anovulatory effects of polychlorinated dibenzo-p-dioxin in rodents. Toxicology. 2001;158:91–107.
Moran FM, Vandevoort CA, Overstreet JW, Lasley BL, Conley AJ. Molecular target of endocrine disruption in human luteinizing granulosa cells by 2,3,7,8-tetrachlorodibenzo-p-dioxin: inhibition of estradiol secretion due to decrease 17a-hydroxylase/17,20-lyase cytochrome P450 expression. Endocrinology. 2003;144:467–73.
Mandal PK. Dioxin: a review of its environmental effects and its aryl hydrocarbon receptor biology. J Comp Physiol B. 2005;175:221–30.
Jablonska O, Piasecka J, Petroff BK, Nynca A, Siawrys G, Wąsowska B, et al. In vitro effects of 2,3,7,8-tetrachlorodibenzo-p-dioxin (TCDD) on ovarian, pituitary, and pineal function in pigs. Theriogenology. 2011;76:921–32.
Albertini DF, Combelles CM, Benecchi E, Carabatsos MJ. Cellular basis for paracrine regulation of ovarian follicle development. Reproduction. 2001;121:647–53.
Jablonska O, Ciereszko RE. The expression of aryl hydrocarbon receptor in porcine ovarian cells. Reprod Domest Anim. 2013;48:710–6.
Sadowska A, Nynca A, Korzeniewska M, Piasecka-Srader J, Jablonska M, Orlowska K, et al. Characterization of porcine Granulosa cell line AVG-16. Folia Biol (Praha). 2015;61:184–94.
Sadowska A, Nynca A, Ruszkowska M, Paukszto L, Myszczynski K, Orlowska K, et al. Transcriptional profiling of porcine granulosa cells exposed to 2,3,7,8-tetrachlorodibenzo-p-dioxin. Chemosphere. 2017;178:368–77.
Tang Z, Wu Y, Yang Y, Yang YT, Wang Z, Yuan J, et al. Comprehensive analysis of long non-coding RNAs highlights their spatio-temporal expression patterns and evolutional conservation in Sus scrofa. Sci Rep. 2017;7:43166.
Horisberger MA. A method for prolonged survival of primary cell lines. In Vitro Cell Dev Biol Anim. 2006;42:143–8.
Gregoraszczuk EL, Wojtowicz AK, Zabielny E, Grochowalski A. Dose-and time dependent effect of 2,3,7,8-tetrachlorodibenzo-p-dioxin (TCDD) on progesterone secretion by porcine luteal cells cultured in vitro. J Physiol Pharmacol. 2000;51:127–35.
Gregoraszczuk EL. Dioxin exposure and porcine reproductive hormonal activity. Cad Saude Publ. 2002;18:453–62.
Jablonska O, Piasecka-Srader J, Nynca A, Kołomycka A, Robak A, Wąsowska B, et al. 2,3,7,8-Tetrachlorodibenzo-p-dioxin alters steroid secretion but does not affect cell viability and the incidence of apoptosis in porcine luteinised granulosa cells. Acta Vet Hung. 2014;62:408–21.
Piasecka-Srader J, Sadowska A, Nynca A, Orlowska K, Jablonska M, Jablonska O, et al. The combined effects of 2,3,7,8-tetrachlorodibenzo-p-dioxin and the phytoestrogen genistein on steroid hormone secretion, AhR and ERβ expression and the incidence of apoptosis in granulosa cells of medium porcine follicles. J Reprod Dev. 2016;62:103–13.
Bolger AM, Lohse M, Usadel B. Trimmomatic: a flexible trimmer for Illumina sequence data. Bioinformatics. 2014;30:2114–20.
Dobin A, Gingeras TR. Mapping RNA-seq reads with STAR. Curr Protoc Bioinformatics. 2015;51:11.14.1–19.
Trapnell C, Roberts A, Goff L, Pertea G, Kim D, Kelley DR, et al. Differential gene and transcript expression analysis of RNA-seq experiments with TopHat and cufflinks. Nat Protoc. 2012;7:562–78.
Pertea M, Pertea GM, Antonescu CM, Chang TC, Mendell JT, Salzberg SL. StringTie enables improved reconstruction of a transcriptome from RNA-seq reads. Nat Biotechnol. 2015;33:290–5.
Pertea M, Kim D, Pertea GM, Leek JT, Salzberg SL. Transcript-level expression analysis of RNA-seq experiments with HISAT, StringTie and Ballgown. Nat Protoc. 2016;11:1650–67.
Zhan S, Dong Y, Zhao W, Guo J, Zhong T, Wang L, et al. Genome-wide identification and characterization of long non-coding RNAs in developmental skeletal muscle of fetal goat. BMC Genomics. 2016;7:666.
Han DX, Sun XL, Fu Y, Wang CJ, Liu JB, Jiang H, et al. Identification of long non-coding RNAs in the immature and mature rat anterior pituitary. Sci Rep. 2017;7:17780.
Kong L, Zhang Y, Ye ZQ, Liu XQ, Zhao SQ, Wei L, et al. CPC: assess the protein-coding potential of transcripts using sequence features and support vector machine. Nucleic Acids Res. 2007;35:W345–9.
Sun L, Luo H, Bu D, Zhao G, Yu K, Zhang C, et al. Utilizing sequence intrinsic composition to classify protein-coding and long non-coding transcripts. Nucleic Acids Res. 2013;41(17):e166.
Wucher V, Legeai F, Hédan B, Rizk G, Lagoutte L, Leeb T, et al. FEELnc: a tool for long non-coding RNA annotation and its application to the dog transcriptome. Nucleic Acids Res. 2017;45(8):e57.
Mistry J, Bateman A, Finn RD. Predicting active site residue annotations in the Pfam database. BMC Bioinformatics. 2007;8:298.
Li A, Zhang J, Zhou Z. PLEK: a tool for predicting long non-coding RNAs and messenger RNAs based on an improved k-mer scheme. BMC Bioinformatics. 2014;15:311.
Kung JT, Colognori D, Lee JT. Long noncoding RNAs: past, present, and future. Genetics. 2013;193:651–69.
Huber W, Carey VJ, Gentleman R, Anders S, Carlson M, Carvalho BS, et al. Orchestrating high-throughput genomic analysis with bioconductor. Nat Methods. 2015;12:115–21.
Ren H, Wang G, Chen L, Jiang J, Liu L, Li N, et al. Genome-wide analysis of long non-coding RNAs at early stage of skin pigmentation in goats (Capra hircus). BMC Genomics. 2016;17:67.
Zhang T, Zhang X, Han K, Zhang G, Wang J, Xie K, et al. Genome-wide analysis of lncRNA and mRNA expression during differentiation of abdominal Preadipocytes in the chicken. G3 (Bethesda). 2017;7:953–66.
Reimand J, Kull M, Peterson H, Hansen J, Vilo J. G:profiler—a web-based toolset for functional profiling of gene lists from large-scale experiments. Nucleic Acids Res. 2007;35:W193–200.
Hartley SW, Mullikin JC. QoRTs: a comprehensive toolset for quality control and data processing of RNA-Seq experiments. BMC Bioinformatics. 2015;16:224.
Hartley SW, Mullikin JC. Detection and visualization of differential splicing in RNA-Seq data with JunctionSeq. Nucleic Acids Res. 2016;44:e127.
Hellemans J, Mortier G, De Paepe A, Speleman F, Vandesompele J. qBase relative quantification framework and software for management and automated analysis of real-time quantitative PCR data. Genome Biol 2007;8(2):R19.
Zhao W, Mu Y, Ma L, Wang C, Tang Z, Yang S, et al. Systematic identification and characterization of long intergenic non-coding RNAs in fetal porcine skeletal muscle development. Sci Rep. 2015;5:8957.
Ran M, Chen B, Li Z, Wu M, Liu X, He C, et al. Systematic identification of long noncoding RNAs in immature and mature porcine testes. Biol Reprod. 2016;94:77.
Xia J, Xin L, Zhu W, Li L, Li C, Wang Y, et al. Characterization of long noncoding RNA transcriptome in highenergy diet induced non-alcoholic steatohepatitis minipigs. Sci Rep. 2016;6:30709.
Melekoglu R, Ciftci O, Cetin A, Basak N, Celik E. The beneficial effects of Montelukast against 2,3,7,8-tetrachlorodibenzo-p-dioxin toxicity in female reproductive system in rats. Acta Cir Bras. 2016;31:557–63.
Bhattacharya P, Keating AF. Impact of environmental exposures on ovarian function and role of xenobiotic metabolism during ovotoxicity. Toxicol Appl Pharmacol. 2017;261:227–35.
Li D, Liu C, Yu H, Zeng X, Xing X, Chenet L, et al. AhR is negatively regulated by miR-203 in response to TCDD or BaP treatment. Toxicol Res. 2014;3:142–51.
Androutsopoulos VP, Tsatsakis AM, Spandidos DA. Cytochrome P450 CYP1A1: wider roles in cancer progression and prevention. BMC Cancer. 2009;9:187.
Whitlock JP Jr. Induction of cytochrome P4501A1. Annu Rev Pharmacol Toxicol. 1999;39:103–25.
Kim S, Dere E, Burgoon LD, Chang CC, Zacharewski TR. Comparative analysis of AhR-mediated TCDD-elicited gene expression in human liver adult stem cells. Toxicol Sci. 2009;112:229–44.
Molcan T, Swigonska S, Orlowska K, Myszczynski K, Nynca A, Sadowska A, et al. Structural-functional adaptations of porcine CYP1A1 to metabolize polychlorinated dibenzo-p-dioxins. Chemosphere. 2017;168:205–16.
Minchenko OH. Inhibition of ERN1 signalling enzyme affects hypoxic regulation of the expression of E2F8, EPAS1, HOXC6, ATF3, TBX3 AND FOXF1 genes in U87 glioma cells. Ukr Biochem J. 2015;87:76–87.
Rashid HO, Yadav RK, Kim HR, Chae HJ. ER stress: autophagy induction, inhibition and selection. Autophagy. 2015;11:1956–77.
Kumar SR, Masood R, Spannuth WA, Singh J, Scehnet J, Kleiber G, et al. The receptor tyrosine kinase EphB4 is overexpressed in ovarian cancer, provides survival signals and predicts poor outcome. Br J Cancer. 2007;96:1083–91.
Jia RX, Duan X, Song SI, Su SC. LIMK1/2 inhibitor LIMKi 3 suppresses porcine oocyte maturation. PeerJ. 2016;4:e2553.
Khorkova O, Myers AJ, Hsiao J, Wahlestedt C. Natural antisense transcripts. Hum Mol Genet. 2014;23(R1):R54–63.
Villegas VE, Zaphiropoulos PG. Neighboring gene regulation by antisense long non-coding RNAs. Int J Mol Sci. 2015;16:3251–66.
Dhakshinamoorthy S, Jaiswal AK. Small maf (MafG and MafK) proteins negatively regulate antioxidant response element-mediated expression and antioxidant induction of the NAD(P)H:Quinone oxidoreductase1 gene. J Biol Chem. 2000;275:40134–41.
Katsuoka F, Motohashi H, Ishii T, Aburatani H, Engel JD, Yamamoto M. Genetic evidence that small maf proteins are essential for the activation of antioxidant response element-dependent genes. Mol Cell Biol. 2005;25:8044–51.
Massrieh W, Derjuga A, Doualla-Bell F, Ku CY, Sanborn BM, Blank V. Regulation of the MAFF transcription factor by proinflammatory cytokines in myometrial cells. Biol Reprod. 2006;74:699–705.
Kim JE. Molecular properties of wild-type and mutant betaIG-H3 proteins. Invest Ophthalmol Vis Sci. 2002;43:656–61.
Riecke K, Grimm D, Shakibaei M, Kossmehl P, Schulze-Tanzil G, Paul M, et al. Low doses of 2,3,7,8-tetrachlorodibenzo-p-dioxin increase transforming growth factor beta and cause myocardial fibrosis in marmosets (Callithrixjacchus). Arch Toxicol. 2002;76:360–6.
Zhang DD, Du JZ, Topolewski J, Wang XM. Review recent progress in identification and characterization of loci associated with sex-linked congenital cataract. Genet Mol Res. 2016;29, 15(3)
Orlowska K, Swigonska S, Sadowska A, Ruszkowska M, Nynca A, Molcan T, et al. The effects of 2,3,7,8-tetrachlorodibenzo-p-dioxin on the proteome of porcine granulosa cells. Chemosphere. 2018;212:170–81.
Bozgeyik E, Igci YZ, Sami Jacksi MF, Arman K, Gurses SA, Bozgeyik I, et al. A novel variable exonic region and differential expression of LINC00663 non-coding RNA in various cancer cell lines and normal human tissue samples. Tumour Biol. 2016;37:8791–8.
Sen R, Doose G, Stadler PT. Rare splice variants in long non-coding RNAs. Non-coding RNA. 2017;3(3):23.
Zaharieva E, Chipman JK, Soller M. Alternative splicing interference by xenobiotics. Toxicology. 2012;296:1–12.
We would like to thank Piotr Ciereszko for English language review.
This study was supported by National Science Centre (2012/05/B/NZ9/03333) and The Ministry of Science and Higher Education in Poland (UWM No. 528.0206.0806).
Availability of data and materials
The datasets analyzed during the current study are available in the NCBI BioProject database under accession number: PRJNA429720 (https://www.ncbi.nlm.nih.gov/bioproject/?term=PRJNA429720).
Ethics approval and consent to participate
Consent for publication
The authors declare that they have no competing interests.
Summary of the mapping of the RNA sequences to the reference genome. (XLSX 13 kb)
Cis-target genes predicted to be regulated by differentially expressed lncRNAs in porcine granulosa cells exposed to TCDD for 3, 12 and 24 h. (XLSX 50 kb)
Negatively co-expressed trans-target genes predicted to be regulated by differentially expressed lncRNAs identified in porcine granulosa cells exposed to TCDD for 3, 12 and 24 h. (XLSX 78 kb)
Positively co-expressed trans-target genes predicted to be regulated by differentially expressed lncRNAs identified in porcine granulosa cells exposed to TCDD for 3, 12 and 24 h. (XLSX 126 kb)
Functional enrichment analysis of the negatively co-expresed trans-target genes predicted to be regulated by differentially expressed lncRNAs identified in porcine granulosa cells exposed to TCDD for 3, 12 and 24 h. (XLSX 54 kb)
Functional enrichment analysis of the positively co-expresed trans-target genes predicted to be regulated by differentially expressed lncRNAs identified in porcine granulosa cells exposed to TCDD for 3, 12 and 24 h. (XLSX 36 kb)
Differential expression of exons and the occurence of splice junction in XLOCs of the lncRNAs identified in porcine granulosa cells exposed to TCDD for 3 and 12 h. (XLSX 16 kb)