Chemerin Impact on Alternative mRNA Transcription in the Porcine Luteal Cells

Chemerin participates in the regulation of processes related to physiological and disorder mechanisms in mammals, including metabolism, obesity, inflammation, and reproduction. In this study, we have investigated chemerin influence on alternative mRNA transcription within the porcine luteal cell transcriptome, such as differential expression of long non-coding RNAs (DELs) and their interactions with differentially expressed genes (DEGs), differences in alternative splicing of transcripts (DASs), and allele-specific expression (ASEs) related to the single nucleotide variants (SNVs) frequency. Luteal cells were collected from gilts during the mid-luteal phase of the oestrous cycle. After in vitro culture of cells un-/treated with chemerin, the total RNA was isolated and sequenced using the high-throughput method. The in silico analyses revealed 24 DELs cis interacting with 6 DEGs and trans-correlated with 300 DEGs, 137 DASs events, and 18 ASEs. The results enabled us to analyse metabolic and signalling pathways in detail, providing new insights into the effects of chemerin on the corpus luteum functions related to inflammatory response, leukocyte infiltration, the occurrence of luteotropic and luteolytic signals (leading to apoptosis and/or necroptosis). Validation of the results using qPCR confirmed the predicted expression changes. Chemerin at physiological concentrations significantly modifies the transcription processes in the porcine luteal cells.


Introduction
Chemerin (CHEM) is a hormone discovered in 1997 during studies analysing the pathogenesis of psoriasis as a product of the RARRES2 gene (originally known as TIG2; short names of genes and macromolecules are explained in the Abbreviations part), whose expression was increased in response to retinoid substances [1]. In 2007, due to its metabolic functions and high secretion by the adipose tissue, CHEM was classified as an adipokine [2,3]. More recently, it was also categorized as a cytokine due to its involvement in immune processes, such as tissue-dependent pro-or anti-inflammatory activity and chemotactic activity on leukocytes, particularly macrophages, dendritic and NK cells [3,4]. Three membrane receptors capable of binding CHEM molecules have been identified in mammals-CMKLR1 (also referred to as ChemR23), GPR1 and CCRL2 [5]. The first two receptors bind the N-terminus of CHEM with similar high affinity and, through binding to may exert different biological functions [38]. Alternatively spliced variants may contain premature stop codons. Protein-coding transcripts carrying a nonsense codon upstream of the terminal exon junction may be targeted for degradation. This mechanism serves to modulate the amount of protein products of a gene without altering its expression in the mammalian cells [39]. Additionally, certain clipped proteins may directly antagonize the full-length proteins activity by retain domains capable of interacting with protein complexes or substrates [40]. Significantly, expression of alternatively splicing transcript variants is often tissue specific [41]. Alternative splicing within genes which products are involved in luteolysis, maternal recognition of pregnancy, and vascular growth and regression is observed in CLs [42,43]. In this publication, we performed a detailed identification of differential alternative splicing events (DASs) appearing in CHEM-treated and native LCs.
Most genes are expressed equally from both alleles; however, some genes are differentially expressed in an allele-dependent manner. Under the influence of regulatory factors such as DNA methylation, histone modifications, and non-coding RNA activity, the production of a protein encoded by one of the alleles of the same gene can be promoted in heterozygotes [44]. This phenomenon is an important genetic factor leading to phenotypic variation induced by, among others, variable environmental conditions [45]. The intensity of angiogenesis in humans has been linked to the occurrence of specific single nucleotide polymorphisms directly related to plasma CHEM levels [46]. This leads to the conclusion that the hormone analysed in the present study is capable of promoting the expression of one of the allelic variants in the affected cells. Allele-specific expression variants (ASEs) were detected out of single nucleotide variants (SNVs) showing imbalance ratios within DEGs.
Despite extensive research analysing the general mechanisms of CHEM action, to date, the details of its influence at the molecular level on specific cells resulting in the modulation of transcription and translation processes remains poorly described. This manuscript is the extension of our previously published analyses of transcriptional profiles on the same research model [19]. Moreover, we have analysed the results obtained by our team using other laboratory techniques on identical research models [17,20,21]. The purpose of this study was to investigate in detail the molecular mechanisms by which CHEM can influence and regulate the transcriptional processes in porcine LCs during the mid-luteal phase of the oestrous cycle. In pursuit of this objective, we focused on alternative mRNA transcription mechanisms, such as DELs interactions with DEGs, as well as DASs within genes, and ASEs within DEGs, which were analysed using high-throughput RNA sequencing data.

Collection of Samples, In Vitro Cell Culture and Chemerin Administration
The experiment was performed on five mature Large White × Polish Landrace female pigs, aged 7-8 months, with a body weight of 120-130 kg, obtained from a private breeding farm. The females, during the mid-luteal phase of the oestrous cycle (days 10-12 of the cycle), were used for the research. Gilts were daily observed for oestrus behaviour in the presence of a boar. The day when the symptoms of the second oestrus appeared was designated as day 0 of the oestrous cycle. Additionally, ovarian morphology was verified to confirm the phase of the cycle [47]. Immediately after humane slaughter, ovaries were collected and transported to the laboratory in ice-cold PBS supplemented with 100 IU/mL penicillin and 100 µg/mL streptomycin.
Luteal cells were isolated from the slaughter material according to the procedure described by Kaminski et al. [48]. Briefly, CLs were dissected from the ovaries, divided into smaller fragments, and dispersed in Ham's F-12 nutrient mixture (Sigma-Aldrich, St. Louis, MO, USA) containing 1% of bovine serum albumin fraction V (BSA; Sigma-Aldrich, St. Louis, MO, USA) and a mix of antibiotics (Sigma-Aldrich, St. Louis, MO, USA). Swine CLs were enzymatically decomposed 4-6 times using 0.125% trypsin (Sigma-Aldrich, St. Louis, MO, USA) solution at 38 • C, centrifuged for 10 min at 300× g at 21 • C, and washed three times. Luteal cells were filtered from trypsinized tissue residues with a Cells 2022, 11, 715 4 of 29 75 µm mesh nylon filter and resuspended in fresh Ham's F-12 mixture. The number of cells was estimated with a haemocytometer and their viability of ≈90% was determined with 0.4% trypan blue dye (Sigma-Aldrich, St. Louis, MO, USA) exclusion.
Luteal cells were resuspended at a concentration of 2 × 10 6 /2 mL in Ham's F-12 medium supplemented with 20% foetal bovine serum (Sigma-Aldrich, St. Louis, MO, USA), 1% BSA, and antibiotics. The cells were pre-incubated in a humidified incubator with 95% air and 5% CO 2 atmosphere for 48 h. The serum-containing culture medium was rejected, and LCs were flushed with serum-free Ham's F-12 nutrient mixture. Porcine LCs were cultured for 24 h in F-12 mixture with 1% BSA, antibiotics, and with or without the recombinant human CHEM (cat. no. 268-10032; RayBiotech, Peachtree Corners, GA, USA) at the concentration of 200 ng/mL for the experimental and control group (n = 5), respectively. The human recombinant CHEM was applied due to the unavailability of porcine CHEM in the onset of laboratory phase of the experiment. In accordance to Du and Leung's study [49], the porcine and human CHEM amino acid homologs show 84% identity. Furthermore, Luangsay et al. [50] have shown that the C-terminus fragment of the mature protein (YFPGQFAFS) is highly conserved in all mammalian species and is highly important for CHEM biological activity. Rytelewska et al. [21] have shown that the identity of the mentioned short amino acid fragment in mammals equals 89% and suggested that CHEM interaction with the cognate receptors is well conserved across mammalian species. The concentration of CHEM was selected on the basis of its blood plasma level in women [46] and gilts [14]. The whole process of performed analysis is shown in Figure 1.

RNA Isolation, Library Preparation and High-Throughput Sequencing Procedure
From in vitro cultures of LCs treated and untreated with CHEM, total RNA was extracted with the use of RNeasy Mini Kit (Qiagen, Germantown, MD, USA) with DNase included in RNase-free DNase Set (Qiagen, Germantown, MD, USA), in accordance with the manufacturer's procedure. The purity (A 260 /A 280 ) and quantity (wavelength 260 nm, A 260 ) of the obtained RNA was estimated spectrophotometrically with an Infinite M200 Pro (Tecan, Männedorf, Switzerland). RNA integrity was tested using the Bioanalyzer 2100 (Agilent Technology, St. Clara, CA, USA). RNA integrity number in a range between 8 and 10 qualified isolated RNA for the transcriptome high-throughput sequencing (RNA-Seq) analysis and the quantitative real-time polymerase chain reaction (qPCR) validations. RNAs extracted from all samples were stored at −80 • C. The strand-specific sequencing libraries were prepared separately from each sample RNA using the TruSeq Stranded mRNA Library Prep Kit (Illumina, San Diego, CA, USA), according to Shen et al. [51] recommendations. The RNA-Seq was performed on the NovaSeq 6000 instrument (Illumina, San Diego, CA, USA). The assumed minimum depth of sequencing was 100 million reads per sample.

In Silico Analyses
The mRNA transcription mechanisms, such as production of lncRNA by cells, alternative splicing of mRNA molecules, and preference for production of mRNA encoded by one of the alleles available in the genetic material, were analysed in order to explore the regulatory mechanisms induced by CHEM in the porcine LCs. The bioinformatic analyses were performed according to Paukszto et al. [52,53], with modifications.

Raw Reads Pre-Processing and Mapping to a Reference Genome
Quality of generated 2 × 150 nt raw paired end reads during sequencing were controlled using FastQC software v. 0.11.8 [54]. Adapters and low-quality regions of raw reads were clipped with the use of Trimmomatic tool v. 0.38 [55] (Q Phred score at 5 and/or 3 ends < 20; average Q Phred score < 30). Reads were trimmed to equal length of 90 nt. Processed reads were rechecked for adapter content and quality with FastQC. The remaining sequences were mapped to the Sus scrofa reference genome with Ensembl annotations v. 11.1.99 [56] with the use of STAR mapper v. 2.7.3a [57] with the operating parameters Cells 2022, 11, 715 5 of 29 recommended by Jakobi [58]. StringTie software v. 1.3.5 [59] with enabled "fr-firststrand" parameter, was applied to annotate and estimate the expression of genes and transcripts. Counts per gene and per transcript were computed using the prepDE Python script provided by the StringTie's authors [60]. The DEGs analysis process was conducted using the Ballgown tool v. 2.18.0 [61], following the procedure described by Makowczenko et al. [19], with the operating parameters: q-value < 0.05 and |log2FC| ≥ 0.5.

Figure 1.
A stage-by-stage pipeline of the research performed within the scope of this manuscript, including the use of the results of differential gene expression analysis performed on the same research model (doi:10.3390/genes11060651). Abbreviations: lncRNAs-long noncoding RNAs, DELs-differentially expressed lncRNAs, ASEs-allele-specific expression variants, DASs-differential alternative splicing events. Identification of lncRNAs in the porcine LCs was performed using multistage workflow. Firstly, low-expressed and protein-coding transcripts were removed from the processed dataset. Secondly, one-exonic and short (length < 200 nt) transcripts were excluded. The coding potential of the remaining transcripts was estimated using four tools: the coding potential assessment tool (CPAT) v. 1.2.4 [62], the coding potential calculator 2 (CPC2) v. 2.0.1 [63], the flexible extraction of lncRNAs (Feelnc) [64], and the predictor of lncRNAs and mRNAs based on an improved k-mer scheme (PLEK) v. 0.2 [65]. Simultaneously, transcripts' sequences were aligned into the protein families (Pfam) database v. 33.1 [66] using "hidden Markov model"-based HMMER software v. 3.3.2 [67]. The transcripts found to be devoid of coding potential (CPAT score < 0.364; CPC2 score < 0; Feelnc coding potential < 0.558; PLEK label = "noncoding") by three of the mentioned tools and with low similarity to Pfam records (e-value > 10 −3 ) were used in subsequent stages of analysis. The differences in detected lncRNAs expression between CHEM-treated and control samples were computed with the Ballgown tool v. 2.18.0 [61] using binominal test, with a cut-off q-value < 0.05 and |log 2 (fold change)| > 0.5. Moreover, the sequences of the obtained DELs were aligned to small RNA model records of the RNA families (Rfam) database v 14.4 [68] with the use of Infernal cmscan v 1.1.3 tool [69].
A further analysis was conducted in R environment v. 4.0.3 [70], with an integrated high-throughput genomic data analysis toolkit-Bioconductor v. 3.12 [71]. It was performed to discover the relationships between the obtained DELs, and mRNAs encoded by DEGs. The relatively close proximity of DELs-DEGs located on the same chromosome were described as cis-interactions. The pairs located in different parts of the genome were also characterized according to transcriptional profiles (trans-action). Analysis of cis-actions was performed within two searching regions on the porcine genome: overlapping, where lncRNA and mRNA-coding genes were located in direct vicinity, and distant, where lncRNA genes were located upstream or downstream of the mRNA-coding genes at a distance of up to 10,000 bp. Trans-actions between DELs and DEGs transcriptional profiles were examined using Pearson's correlation coefficient, with a cut-off |r| > 0.9 and p-value < 0.05. In addition, the binding strengths of lncRNAs to potential mRNA targets were tested using LncTar v. 1.0 software [72] with the normalized free energy (ndG) calculation. Moreover, the probability of hydrogen-bonding propensities, and van der Waal's interaction between the secondary structures of the detected lncRNAs and potential protein targets were calculated using lncPro v. 1.0 tool [73]. In subsequent analyses, only those trans-actions for which LncTar indicated an ndG < −0.1 or lncPro showed a probability > 0.9 were included, providing more detailed information on the DEL-DEG trans-actions detected in the previous steps of the analysis. The identified interactions between DELs and DEGs were visualized using Cytoscape software v. 3.8.3 [74].

Differential Alternative Splicing Events Analysis
Alternative splicing events were predicted with the super-fast pipeline for alternative splicing analysis 2 (SUPPA2) tool v. 2.3 [75]. Briefly, equal length processed reads uniquely mapped against the reference genome were retrieved, and, subsequently, were re-mapped against S. scrofa reference transcriptome using Salmon mapper v. 1.1.0 [76]. Differential alternative splicing events between experimental and control groups were statistically tested and the percentage of splicing inclusions (PSI) for all splicing events was calculated. Detected DASs were considered as statistically significant with p-value < 0.05, and |∆PSI| > 0.1. All discovered DASs were classified into seven categories by SUPPA2: alternative 5 splice site (A5), alternative 3 splice site (A3), mutually exclusive exons (MX), retained intron (RI), skipping exon (SE), alternative first exon (AF) and alternative last exon (AL). Significant DASs were visualized using the ggsashimi Python tool v. 0.5.0 [77].

Single Nucleotide Variants Identification and Allele-Specific Expression Variants Analysis
The single nucleotide variants in transcripts were discovered by aligning RNA-Seq reads to the genome reference sequence. Calling of SNVs with multi-sample data was conducted using the pipeline consisting of Picard software v. 2.6.0 [78], the replicate multivariate analysis of transcript splicing and discovery of differential variants in RNA (rMATS-DVR) tool v. 1.0.0 [79], and the gold standard genome analysis toolkit (GATK) v. 3.6.0 [80]. Known genomic positions of porcine SNVs were obtained from Ensembl database v. 11.1.99 [56]. All previously created BAM files were recalibrated with the Picard tool. The individual rMATS-DVR modules firstly identified the potential occurrence of SNVs in all sample replicates, and then detected the differences in allele frequencies between CHEM-treated and control sample groups. Low-quality and disrupted SNVs were filtered out on the basis of GATK standard parameters: total depth of base coverage (>10), root mean square mapping quality (>40), quality by depth (>2), mapping quality rank sum (>−12.5), rank sum test for relative positioning of reference versus alternative alleles within reads (>−8). Single nucleotide variants with alternative allele fraction (AAF) > 0 within at least half of the RNA-Seq samples were selected for further analyses. Allelic expression bias of discovered SNVs between experimental and control groups were considered as statistically significant with |∆AAF| > 0.1 and false discovery rate (FDR) < 0.05. An allelic imbalance ratio of ASEs candidates was confirmed using the Differential alternative splicing events between experimental and control groups were statistically tested and the percentage of splicing inclusions (PSI) for all splicing events was calculated. Detected DASs were considered as statistically significant with p-value < 0.05, and |ΔPSI| > 0.1. All discovered DASs were classified into seven categories by SUPPA2: alternative 5′ splice site (A5), alternative 3′ splice site (A3), mutually exclusive exons (MX), retained intron (RI), skipping exon (SE), alternative first exon (AF) and alternative last exon (AL). Significant DASs were visualized using the ggsashimi Python tool v. 0.5.0 [77].

Single Nucleotide Variants Identification and Allele-Specific Expression Variants Analysis
The single nucleotide variants in transcripts were discovered by aligning RNA-Seq reads to the genome reference sequence. Calling of SNVs with multi-sample data was conducted using the pipeline consisting of Picard software v. 2.6.0 [78], the replicate multivariate analysis of transcript splicing and discovery of differential variants in RNA (rMATS-DVR) tool v. 1.0.0 [79], and the gold standard genome analysis toolkit (GATK) v. 3.6.0 [80]. Known genomic positions of porcine SNVs were obtained from Ensembl database v. 11.1.99 [56]. All previously created BAM files were recalibrated with the Picard tool. The individual rMATS-DVR modules firstly identified the potential occurrence of SNVs in all sample replicates, and then detected the differences in allele frequencies between CHEM-treated and control sample groups. Low-quality and disrupted SNVs were filtered out on the basis of GATK standard parameters: total depth of base coverage (>10), root mean square mapping quality (>40), quality by depth (>2), mapping quality rank sum (>−12.5), rank sum test for relative positioning of reference versus alternative alleles within reads (>−8). Single nucleotide variants with alternative allele fraction (AAF) > 0 within at least half of the RNA-Seq samples were selected for further analyses. Allelic expression bias of discovered SNVs between experimental and control groups were considered as statistically significant with |ΔAAF| > 0.1 and false discovery rate (FDR) < 0.05. An allelic imbalance ratio of ASEs candidates was confirmed using the ꭕ 2 goodness-2 goodness-of-fit test (p-value < 0.05). The exact location of identified ASEs within genes' regions and the effect of single nucleotide mutations on the transcription process of specific proteins were determined using the Ensembl variant effect predictor (VEP) web tool [81].

Functional Annotation of Target Genes
To analyse the functions and involvement of potential target DEGs of the identified DELs, DASs and ASEs in biological processes, the Ko-based annotation system (KOBAS) web tool v. 3.0 [82] searching gene ontology (GO) [83,84], the Reactome Knowledgebase [85], and the Kyoto encyclopaedia of genes and genomes (KEGG) [86] databases was implemented. The gene enrichment was performed for every analysed phenomenon individually, with the use of Fisher's exact test with Benjamini and Hochberg correction (FDR < 0.05). Furthermore, KEGG analysis of target genes were clustered to generate new molecular network interaction maps using a Pathview v. 1.30.1 R package [87]. The statistical significance values were recalculated for the regrouped data (FDR < 0.05).

Quantitative Real-Time PCR Validation
To validate the results discovered during DELs analysis, the qPCR method was applied. The 500 ng of each RNA sample used in the RNA-Seq were rewritten to cDNA, yielding a total volume of 10 µL of each mixture, with the use of Omniscript RT Kit (Qiagen, Germantown, MD, USA), and 0.5 µg of oligo(dT) (Roche, Penzberg, Germany). The reaction of reverse transcription was conducted at 37 • C for 1 h and was terminated by incubation at 93 • C for 5 min.
The quantitative real-time PCR was performed using an AriaMx Real-Time PCR System (Agilent Technology, St. Clara, CA, USA). The constitutively expressed ACTB and GAPDH as reference genes were used. Primer sequences for target transcripts (CL.9638.3, CL.12742.3, ENSSSCT00000075362 and ENSSSCT00000078829) were developed using Primer Express software 3 (Applied Biosystems, Waltham, MA, USA). The primer sequences of reference and target genes are listed in Table 1. Rection mixtures with a final volume of 20 µL consisted of 30 ng prepared cDNA, 200 nM of the forward and reverse primers, 12.5 µL of the Power SYBR Green PCR Master Mix (Applied Biosystems, Waltham, MA, USA) and RNase-free water. Quantitative real-time PCR was performed according to the procedure: preliminary cDNA denaturation and enzymes activation at 95 • C for 3 min, 40 cycles of denaturation at 95 • C for 30 s, annealing at 66 • C for 1 min, and elongation at 72 • C for 45 s. For ACTB primers the annealing temperature was lowered to 61 • C, and for EN-SSSCT00000075362 primers preliminary denaturation was extended to 10 min and the annealing temperature was lowered to 60 • C. Negative controls were prepared by replacing the cDNA with water. The reactions were performed in technical duplicate for each sample. The amplification specificity was examined at the end of qPCR by melting-curve analysis.
The relative expression levels of validated transcripts were calculated using the comparative cycle threshold (∆∆CT) method and normalized using the geometrical means of the reference genes expression [88]. The results of qPCR were statistically processed using Student's t-test with a significant p-value < 0.05 in the R environment v. 4.0.3 [70] and were presented as mean values ± standard error of the mean (SEM).

Overall Statistics of RNA-Seq Data Mapping
The total number of raw reads obtained across all samples by the deep sequencing of transcriptome was 1,154,494,264, and the total number of reads obtained after the pre-processing step was 1,024,385,352. The total number of processed reads mapped to the reference porcine genome (S. scrofa v. 11.1.99) was 1,021,940,610 (Table 2). Of them, 96.25% were mapped uniquely and 3.18% were mapped to multiple loci. Among all RNA-Seq libraries used in this study, the average number of reads' bases mapped to the coding DNA sequences (CDS) were 59.52%, untranslated regions (UTR) were 23.32%, intronic regions were 3.35%, and intergenic locations were 13.81%. The detailed statistics of the sequencing and pre-processing data, principal component analysis and sample-to-sample distance matrix were summarized by Makowczenko et al. [19].

Long Non-Coding RNA Identification and Cis-/Trans-Acting on Protein-Coding Genes
The first steps of filtration involved removing from the dataset sequences designated as "protein-coding biotype" in the Ensembl database, as well as single-exon RNAs and RNAs shorter than 200 nt. This process resulted in the detection of 25,252 potential lncRNA candidates in the porcine LCs. Of the surviving transcripts, 9340 were labelled as "lncRNA biotype" in the Ensembl database. The remaining unknown sequences in the data scope were scanned using four tools-CPAT, CPC2, Feelnc and PLEK, and in consequence reducing the number of novel lncRNAs to 1759. Of these, 107 survived alignments to the Pfam database using HMMER tool ( Figure 2). Analysis of genes expression differences between CHEM-treated and control samples revealed that 24 of the detected lncRNAs (encoded by 12 genes) undergo expression modifications in LCs under the influence of the CHEM, including 20 previously known and 4 unidentified lncRNAs (Figure 3). Among the identified DELs, 16 were down-regulated while 8 were up-regulated. None of the detected transcripts were matched by Infernal cmscan to small RNA models in the Rfam database. The detected DELs are summarized in Table 3. Cis-interaction analysis revealed 11 cases of DELs affecting 6 DEGs based on their colocalization in the porcine genome. Of these, DELs overlapped target gene loci in 3 interactions, 5 interactions occurred in DELs located within 10,000 bp downstream from DEGs, and 3 interactions occurred up to 10,000 bp upstream from the protein-coding genes ( Table 4). The analysis of trans-actions between DELs and DEGs revealed the 9295 interaction events occurring with 469 DEGs (assuming a high level of Pearson's correlation coefficient). Additional analysis of mRNA-lncRNA and protein-lncRNA bindings revealed, respectively, 1393 and 115 cases, restricting the number of previously discovered trans-interactions to 1440 performed with 300 DEGs (Table S1). All detected cis-and trans-interactions between DELs and DEGs are summarized in Figure 4.

Differential Alternative Splicing Events of Differentially Expressed Genes
The adopted procedure, incorporating the SUPPA2 tool, allowed the detection of 43,918 alternative splicing events, including 137 DASs resulting from the comparison of CHEM-treated and control samples. All detected DASs are localized within 92 DEGs. Among all detected DASs, 13 were classified as A5, 12 as A3, 2 as MX, 13 as RI, 21 as SE, 60 as AF and 16 as AL ( Figure 5). Selected events of alternative splicing occurring within the CXCL12, MCL1, FHL1 and SLA-DQB1 genes are visualized in Figure 6, while all identified cases are summarized in Table S2.

Single Nucleotide Variant Calling and Allele-Specific Expression Variants
During the in silico analyses of ten RNA-Seq libraries, 151,603 SNVs were called for the transcriptome of porcine mid-luteal LCs. After the multi-step filtering procedure, based on the standard GATK parameters and the AAF occurrence in at least half of the samples, 37 SNVs were obtained and directed to the next steps of analysis. Among all remaining SNVs, 18 showed statistically significant allelic imbalance occurring at the same loci (|∆AAF| > 0.1 and FDR < 0.05) between CHEM-treated and control samples. In accordance with the VEP annotations, 13 ASEs is novel and 5 is previously known: rs1110933145, rs322079801, rs332564681, rs346401196, and rs55619031. Detected ASEs overlapped 38 transcripts encoded by 16 genes. Localization consequences of individual ASEs within the genes and transcripts components were summarized in Figure 7. As annotated by VEP, 11 ASEs affected non-coding regions of transcripts encoded by 11 genes with consequences difficult to predict. Four ASEs affected CDS regions of transcripts encoded EIF3F, OSBP, PIK3C2A and SLC25A24 genes with low probability to change coded protein. Four ASEs affected transcripts encoded by EIF3F, HPS5 and OSBP genes non-disruptively with chance to change protein effectiveness. Impact of 1 ASE on transcripts encoded by CXCL2 gene is disruptive and causes loss of protein function by modifying the splice acceptor of second intron. According to the results of the sorting intolerant from tolerant (SIFT) module [91] implemented in VEP, of the 4 missense ASEs that can affect the protein product conformation, 2 have a deleterious effect (on HPS5 and OSBP genes) and 2 could be tolerated (EIF3F and OSBP genes). The detailed features of discovered ASEs, their biological impact on genes and protein translation process were summarized in Table S3.

Functional Annotation of Target Protein-Coding Genes
In order to summarize and combine the results obtained by analysing three different biological processes occurring at the transcription RNA level, an additional research step was carried out. Initially, the KOBAS tool was used to enrich terms provided in GO, the Reactome Knowledgebase and KEGG databases, separately for each of the three biological phenomena analysed. Gene ontology enrichment analysis of genes cis interacting with DELs, genes trans-interacting with DELs, genes containing DASs and/or ASEs showed a significant contribution of, respectively, 5, 273, 39 and 10 genes encompassed to 30, 159, 2 and 39 biological process (BP) terms, 5, 34, 1 and 17 molecular function (MF) terms, and 4, 24, 3 and 12 cellular component (CC) terms (Table S4). The Reactome Knowledgebase analysis demonstrated the significant enrichment of 3 molecular pathways by genes cis interacting with DELs, 20 pathways by genes trans-interacting with DELs, and 1 pathway by genes containing ASEs (Table S4). The Reactome Knowledgebase enrichment analysis did not identify any statistically significant molecular pathways in which the genes containing DASs were involved.
To predict the detailed biological effect of the observations described in the previous sections, the protein-coding genes related to DELs, DASs and ASEs were distinguished in maps of metabolic or signal transduction KEGG pathways. The 84 enriched pathways indicated the statistical significance. Due to the unrelatedness of the studied tissue and the purpose of this study, the pathways associated with pathological processes, such as autoimmune diseases, carcinogenesis, and infections initiated by bacteria, prions, protists and viruses were not analysed. Details of the remaining 35 pathways associated with physiological mechanisms are summarized in Table 5 and Figures S1-S18. Additionally, important for CLs functioning, "Arachidonic acid metabolism" pathway that did not reach statistical significance during KEGG enrichment was examined ( Figure S19). Table 5. Summary of KEGG pathways that are statistically implicated by protein-coding genes associated with DELs, DASs, and ASEs occurring in the chemerin-treated porcine luteal cells during days 10-12 of the oestrous cycle. The 'Visualization' column contains the numbers of figures available in the supplemental materials presenting molecular pathways included in the discussion.

Pathway Name Pathway ID Input (Background) Number FDR Visualization
Cytokine-cytokine receptor interaction

Discussion
This study represented the first ever attempt to describe the CHEM impact on alternative mRNA transcription mechanisms viz. DELs, DASs, and ASEs in the mammalian CLs. The aforementioned analyses were performed on data derived from high-throughput sequencing of the in vitro cultured mid-luteal porcine LCs transcriptomes. Inference of the physiological effects of the molecular processes modified by CHEM was carried out on the basis of interaction or direct association of individual events with protein-coding genes, whose transcriptional profiles were described by Makowczenko et al. [19]. Briefly, we identified 24 DELs that underwent 11 cis interactions with 6 DEGs and 1440 trans-interactions with 300 DEGs, 137 DASs across 92 genes, and 18 ASEs present on transcripts encoded by 16 genes.
One of the three CHEM membrane receptors is CCRL2, often described as an "atypical" receptor due to the lack of chemical signal transduction [92]. Monnier et al. [8] reported a direct association between activation of the MAPK (including NF-κB transcription factors) and Jak/STAT signalling pathways ( Figures S5 and S9, respectively) by pro-inflammatory factors and the up-regulation of CCRL2 gene expression in the human and mouse endothelial cells. In a recent study, we demonstrated the presence of CCRL2 protein in the porcine LCs on days 10-12 of the oestrous cycle [17] and postulated activation of both mentioned signal transduction pathways under the influence of CHEM [19]. Admittedly, we did not note statistically significant differences in CCRL2 mRNA abundance in CHEM-treated LCs compared to the control group [19]. However, in the current study, we can observe significant differences in the process of alternative splicing, classified as AF with (∆PSI = −0.48), resulting in increased proportion of CCRL2 transcript variants containing first exon of 252 nt length in comparison to those containing first exon of 56 nt length. Both alternative sequences contain a fragment of 5 UTR. It is highly likely that the observed alteration of the first exon may affect the stability of transcriptional variants and/or the translation process, affecting the concentration of CCRL2 proteins on the surface of LCs under the CHEM influence. The effect may promote a chemotactic gradient of CHEM through its binding by CCRL2 proteins, increasing the migration of leukocytes into the CLs, as described in the case of other cell types, such as endothelial and epithelial cells [7,8,93,94]. The above observation requires further detailed studies to elucidate the reported phenomenon.
Data on the relationship of immune cells and produced by them inflammatory cytokines to the luteolysis management still remain limited in pigs. Nevertheless, according to the hypothesis by Pate and Keyes [95], immune cells are actively involved in the loss of steroidogenic potential and the destruction of CLs. Standaert et al. [96] have shown that in the porcine CLs, immune cells are present in addition to the natively present cells and their number changes depend on the ongoing stage of the oestrous cycle. Ziecik et al. [97] have shown that genes related to lymphocyte infiltration are activated in CLs with the onset of functional luteal regression in gilts. Moreover, increased monocytes/macrophages abundance during the late-luteal phase is directly related to the presence of pro-inflammatory cytokines that regulate CLs apoptosis and structural luteolysis. In our previous study, we have suggested that CHEM may indirectly regulate the process of luteolysis by influencing the immune system cells motility as the chemotactic factor and agent that modifies the levels of pro-and anti-inflammatory cytokines [17]. Among the data analysed, we noted a significant effect of DELs and ASEs on genes involved in TNF signalling pathway ( Figure S2), the activation of which results in increased production of the cytokines, chemokines, and adhesion factors analysed in the following paragraphs. The CHEM effects on the genes transcription and processes affecting the inflammatory status of porcine CLs during the mid-luteal phase of the oestrous cycle are summarized in Figure 9.
The chemokines CCL2 and CCL4 are known as chemoattractants and activators of immune cells, such as monocytes, T cells, and NK cells [98,99]. The main function of CCL5 is to activate diverse lymphocyte populations, but it can also stimulate eosinophils, basophils, and dendritic cells [100]. Witek et al. [101] in a recent study on the porcine CLs found that expression of genes encoding CCL2, CCL4, and CCL5 increases during the late-luteal phase of the oestrous cycle compared to the mid-luteal phase. During our previous study [19], we found that the concentration of CCL2, CCL4, and CCL5 mRNAs increased in response to CHEM. The current study revealed the negative correlations of these genes' transcriptional profiles with DELs' profiles. We have also identified this DEGs' intermolecular bonding with DELs at the mRNA-lncRNA level (one for CCL2 and CCL5 each), and protein-lncRNA level (2 for CCL2 and 11 for CCL4). The presence of CHEM in the environment, via effects on genes encoding the described chemokines, may contribute to leukocyte infiltration into CLs. In situ, based on additional signals from LCs, immune cells may initiate luteal regression. Figure 9. The summary of the CHEM effects on the genes transcription and processes affecting the inflammatory status of porcine CLs during the mid-luteal phase of the oestrous cycle. The arrows next to the processes names or genes abbreviations symbolize the effect of the CHEM influence: ↑ (red)-stimulatory, ↓ (green)-inhibitory. Details are described in the text. Abbreviations: CCL 2/4/5-C-C motif chemokine ligand 2/4/5 genes, CL-corpus luteum, CTSS-cathepsin S gene, ICAM1-intercellular adhesion molecule 1 gene, IL 1A/6/8-interleukin 1A/6/8 genes, MHCIImajor histocompatibility complex II, NF-κB-nuclear factor κ-light-chain-enhancer of activated B cells, PSMB 8/9/10-proteasome 20S subunit β 8/9/10 genes, TGFB3-transforming growth factor β3 gene, TRAF3IP2-TNF receptor associated factor 3 interacting protein 2 gene. This thesis is also supported by the increase in the levels of IL1A and TRAF3IP2 mRNAs [19] and their interactions at the protein-lncRNA level with DELs (6 actions per gene). Both genes are involved in the expression induction of, among others, IL1A, IL6, IL8 (also named CXCL8) and TGFB3, as the result of the NF-κB signal transduction pathway activation ( Figure S15). Under CHEM influence, the mRNAs abundance of these interleukins was strongly up-regulated, whereas the transcription of TGFB3 downregulated [19]. We identified protein-lncRNA interactions with DELs for each of the four genes (6, 3, 3, and 2 trans-actions respectively). The expression profiles of these factors are distinctive of cellular senescence [102,103]. The senescent cells can potentially attract and activate different immune cells subsets (such as NK cells, B cells, T cells, neutrophils, dendritic cells or monocytes/macrophages), which in turn, can be engaged in the immune response and removal of senescent cells [104]. Interestingly, two research teams reported ICAM1 expression up-regulation in the senescent cells [105,106], which may facilitate communication with immune cells and induce luteal regression, as described by Olson et al. in rats [107]. In our previous studies, ICAM1 mRNA concentration was up-regulated [19]. Furthermore, 6 trans-interactions with DELs at protein-lncRNA level were identified.
Among all chemokines' gene transcriptional profiles analysed in this and previous work [19], only the transcription of CXCL12 (encoding SDF-1 chemokine protein) decreased, with a simultaneous increase in the mRNA abundance of its receptor-CXCR4. Within CXCL12, we identified DAS event classified as A3 (∆PSI = −0.37), resulting in the formation of transcripts lacking 2597 nt fragment within the penultimate exon containing CDS-ending sequence and partial 3 UTR sequence ( Figure 6A). Additionally, in CHEM-treated cells we identified an ASE site numbered rs346401196 in the Ensembl database (∆AAF = 0.187), resulting in an increase in the frequency of nucleotide C, at site A within the last exon containing a partial 3 UTR sequence (or intron in some transcript variants; Figure 7). Expression of CXCL12 and CXCR4 genes was found in the human pre-ovulatory GCs [108] as well as in the sheep LCs [109,110]. Available reports focus on SDF-1 and CXCR4 proteins contribution to chemotactic activity of immune cells and the resulting consequences for folliculogenesis [108,111,112]. Although little is known about these proteins' involvement in CLs regulation, McIntosh et al. [110] observed the increased expression of both genes in LCs during early stages of pregnancy. The researchers indicated the luteotropic effect of SDF-1 conductive to the establishment and maintenance of intrauterine environment suitable for conceptus implantation. Furthermore, Nishigaki et al. hypothesized that SDF-1 may facilitate luteinization of human preovulatory follicles [113]. The potential changes in the amount of SDF-1 and CXCR4 proteins in CHEM-treated LCs may suggest CHEM influence on the chemokine action and, indirectly, leukocytes infiltration of the luteal tissue. This phenomenon may be related to the discovery of de Poorter et al. [114] regarding the possibility of heteromer formation between the CHEM receptor CMKLR1 and CXCR4. This may have important implications for the transduction of signals by both types of receptors, mainly cross-inhibition [3]. Furthermore, the same research team found that CHEM co-incubated with SDF-1 decreases their specific binding to mouse bone marrowderived dendritic cells [114]. The SDF-1 and CXCR4 protein pair is involved in a number of molecular processes, such as "Cytokine-cytokine receptor interaction" (Figure S1), "NFkappa B signalling pathway" (Figure S5), "Regulation of actin cytoskeleton" (Figure S13), and 'Leukocyte transendothelial migration' (Figure S14).
Noteworthy are the results obtained during our previous research [19], in which the abundance of the HLA-DRA, SLA-DMA, and SLA-DQB1 mRNAs decreased under CHEM treatment in the swine mid-luteal LCs. Furthermore, we detected 4 (2 positively and 2 negatively corelated), 1 (positively corelated), and 8 (7 positively and 1 negatively corelated) trans-interactions of mentioned genes with DELs, respectively, and DAS event qualified as SE within the SLA-DQB1 gene, resulting in the absence of exon 4 in the final transcripts (∆PSI = 0.26; Figure 6D). The proteins encoded by these genes are antigens belonging to the MHC II group [115]. Members of this group are present on the small and large LCs, and their expression increases significantly after the onset of luteal regression in humans, cattle and ewes [116][117][118]. MHC II class peptides are recognized by CD4+ T cells ( Figure S10), activation, of which results in their clonal expansion and leads to the release of cytokines such as TNFα and IFNγ [95]. Therefore, when the expression of these genes is inhibited and/or the formation of shorter length transcript variants is promoted, the amount of MHC II membrane proteins is also reduced. This may lead to suppressed CD4+ T lymphocytes activation and inhibition of secretion of luteolytic cytokines into the CLs microenvironment.
On the other hand, all three genes-PSMB9, PSMB10, PSMB8, encoding the β1i, β2i, and β5i particles of the immunoproteasome 20S catalytic subunit, respectively, and PSME1 and PSME2 genes encoding two particles, PA28α and PA28β, building the immunoproteasome 19S regulatory subunit in cytosol [119], were all up-regulated in CHEM-treated swine LCs [19]. Additionally, for PSMB10 (1 positively correlated), PSMB8 (4 negatively correlated), and PSME2 (1 positively and 3 negatively correlated), trans-interactions with DELs were detected ( Figure S17). Immunoproteasomes degrade ubiquitin-tagged proteins by preparing the hydrophobic C-ends of the proteins to fit into the grooves of MHC I complexes ( Figure S10), allowing the prepared fragments to be exposed to CD8+ T lymphocytes which may lead to cell destruction (apoptosis) through a cytotoxic effect [119]. The expression of proteins building immunoproteasomes can significantly increase in cells exposed to pro-inflammatory cytokines, or oxidative stressors [119,120]. Most notably, Luo and Wiltbank [121] reported increased expression of PSMB9, PSMB8, PSMB10, and CTSS genes in the luteinized bovine GCs co-cultured with T cells. Interestingly, Cash et al. [122] demonstrated the conversion of CHEM molecules into shorter fragments with potent anti-inflammatory effects by a number of cysteine proteases, including CTSS. During the current study, we identified two cases of DASs in the CTSS gene, classified as an AF (∆PSI = −0.30 and −0.18) which may result in modification of the mature protein. The presented modifications in the expression of genes which products build immunoproteasomes, as well as interactions with DELs and the occurrence of DASs in the porcine mid-luteal LCs, would suggest the involvement of CHEM in cellular response that can lead to degradation of CLs.
Among the obtained functional analysis results, we observe the high proportion of signalling pathways and metabolic processes directly related to immune responses, such as "IL-17 signalling pathway" (Figure S3), "NOD-like receptor signalling pathway" (Figure S4), "Toll-like receptor signalling pathway" (Figure S7), "Cell adhesion molecules (CAMs)" (Figure S11), and 'Phagosome' ( Figure S18). This is further evidence of the CHEM influence on the intensification of processes related to immune mechanisms in the microenvironment of the porcine mid-luteal CLs. The limitation of the number of immune cells in in vitro cell culture research models makes it impossible to observe all the nuances of the reported molecular processes that occur in living organisms, and further studies to shed new light on the mentioned phenomena are worthwhile.
Within FHL1, the gene linked to the Jak/STAT signal transduction pathway, we found five forms of statistically significant DASs categorized as AF, whose occurrence was caused by CHEM (∆PSI ∈ {−0.22} ∪ [0.13, 0.30]; Figure 6B). The biological function of the protein encoded by this gene is poorly understood in the processes unrelated to carcinogenesis. However, Fimia et al. [123] demonstrated expression of FHL1 in the mammalian ovarian cells. Matulis & Mayo [124] linked the knock-down of this gene to the downregulation of the NR5A transcription factor in the mouse GCs, which directly affects the production of proteins important for steroidogenesis, such as CYP11A1 (also referred to as P450scc) catalysing the conversion of cholesterol to pregnenolone, the direct substrate of progesterone. Also, accordingly to Ziecik et al. [97], expression of NR5A1 gene was down-regulated on day 14 vs. day 12 (late-luteal vs. mid-luteal phase) of the estrous cycle in porcine CLs, which may be associated with the acquisition of luteolytic sensitivity. The appearance of DASs within FHL1 affecting the decreased activity of the proteins produced would explain the decrease in CYP11A1 mRNA concentration and its protein production under the CHEM influence in ovarian cells [10].
We can infer a positive effect of CHEM on angiogenesis in the mid-luteal CLs of gilts based on observations both at the transcriptomic level-increase in VEGFA mRNA level [19], negative correlation with seven DELs showing six protein-lncRNA interactions and one mRNA-lncRNA interaction, and at the protein level-increase in VEGFA secretion [21]. Vascular endothelial growth factor A through its signal transduction pathway ( Figure S16), plays a major role in endothelial cell survival, vascular stabilization and integrity [125][126][127]. This finding seems to be in line with CHEM effect on PGE 2 production. By treating the porcine mid-luteal LCs with CHEM, we observed an increase in mRNA abundance of PTGS2 and PTGES creating the pathway from the primary substrate-arachidonic acid to the production of luteotropic factor-PGE 2 [19]. Moreover, we have found the protein-lncRNA trans-interactions associated with negative correlation of this DEGs' transcriptional profiles with 3 and 2 DELs, respectively ( Figure S19). According to the "two signal switch" hypothesis, PGE 2 , acting on LCs around day 12 of the oestrous cycle in pigs, prevents the acquisition of luteolytic sensitivity by activating the so-called "rescue switch" that results in the maintenance of angiogenesis, steroidogenesis, and LCs' survival [97].
In previous studies, we have reported a stimulating effect of CHEM on concentration of CASP3 mRNA encoding the major execution protease of the apoptosis process in the porcine LCs [19]. During the present study, we identified the lncRNA molecule CL.9638.3. Its expression was inhibited by CHEM and was negatively correlated with that of CASP3 ( Figure S12). For both molecules, we additionally detected the probability of combining protein-lncRNA spatial structures of 94.66%. However, it should be noted that the recent study by our team, in which the ELISA method was implemented, showed no effect of CHEM on the amount of CASP3 protein produced by the mid-luteal LCs of gilts [21]. This allows us to assume that in spite of CHEM influence on CASP3 generation at the transcriptomic level, this effect is neutralized during the next steps of its production. However, it cannot be ruled out that CHEM may modify CASP3 protein changing its activity. The evidence of the CHEM influence on the process of apoptosis was the increased expression of two lncRNAs-CL.9638.3 and CL.2666.5-showing high probability of binding to MCL1 proteins (95.59% and 90.56%, respectively) and exhibiting significant correlation values of transcription levels (r = −0.97 and 0.96). Furthermore, within the MCL1 gene, we identified DAS event, classified as AF (∆PSI = −0.18; Figure 6C). In the previous publication [19] we demonstrated the MCL1 mRNA abundance under the CHEM-treatment increased. Interestingly, Shee-Uan et al. [128] found that the increase in MCL1 expression in LCs prevents apoptosis and increases cell viability, and thus may play an important role in regulating the lifespan of CLs. It seems, therefore, that CHEM effect on MCL1 production and apoptosis process in porcine CLs is ambiguous. It is also worth paying attention to the slight increase in the percentage of splicing inclusions classified as RI event within DUSP4 (∆PSI = 0.11), whose transcription enhancement under the influence of CHEM we previously observed [19]. The aforementioned CHEM effect consisted in generation suppression of transcript variants containing the last exon with the terminal CDS fragment and 3 UTR at the advantage of the variant having two exons encoding analogous regions of the transcript. The DUSP4 encodes the nuclear phosphatase involved in ERK and p38 dephosphorylation ( Figure S6), and thereby controlling the activity of the MAPK-ERK1/2-p38 signal transduction pathway [129,130]. Increased expression of DUSP4 results in providing an anti-apoptotic effect and promotes cell survival [131,132].
An interesting observation that requires further detailed study seems to be the significant effect of CHEM on genes involved in necroptosis (a caspase-independent form of programmed cell death reminiscent of necrosis) in the mid-luteal LCs of gilts ( Figure S8). Hojo et al. [133], studying bovine CLs, concluded that luteolysis is an acute tissue damage process, whereas the mechanism of apoptosis itself is not acute. It seems that an efficient regression of the CLs must consist of a number of other processes, including necroptosis. Markers of necroptosis are the proteins RIPK1 and RIPK3, whose activation increases during luteal regression and which, forming homo-and heteromers, phosphorylate a number of proteins that are necroptosis executors [133,134]. No differences in the mRNA abundance and occurrence of transcription-associated phenomena (DELs interactions, DASs and ASEs) were found for any gene encoding the aforementioned proteins [19]. Nonetheless, we observed the increase in the transcription of some genes and an intensification of transcription-associated mechanisms in genes involved in the signal transduction pathways, that are also involved in both ubiquitination (TNF pathway; Figure S2) and induction of RIPK1 expression (Jak/STAT pathway; Figure S9). Moreover, CHEM increased TNFAIP3 mRNA concentration [19] and revealed 6 negatively correlated interactions with DELs at the protein-lncRNA level. The protein encoded by this gene shows the ability to deubiquitinate the other pro-necroptotic factor, RIPK3. These finding indicate the possibility of CHEM-induced specific preparation of cells for necroptosis which could be triggered by additional factors.

Conclusions
Performing in silico analyses of alternative mRNA transcription related to CHEM effects on the porcine LCs enabled a more detailed insight and description of the sophisticated and subtle molecular mechanisms controlling cells functions. We have gained further arguments indicating that CHEM affects the MAPK (including NF-κB transcription factors) and Jak/STAT signalling pathways, local generation and action of cytokines and chemokines, as well as essential for CLs processes such as angiogenesis, prostaglandins production and steroidogenesis. The phenomena analysed at the molecular level indicate activation of both pro-apoptotic (or even pro-necroptotic) and pro-survival mechanisms, supporting CHEM engagement in the control of CLs lifespan.
Supplementary Materials: The following are available online at https://www.mdpi.com/article/ 10.3390/cells11040715/s1, Table S1: Overview of trans-interactions between DELs and DEGs, Table  S2: Overview of detected DASs, Table S3: Overview of detected ASEs, Table S4: Results of GO and Reactome functional analysis of DEGs acting with DELs or within which DASs or ASEs were identified, Figure S1: Cytokine-cytokine receptor interaction with DELs' interactions, DASs and ASEs mapped, Figure S2: TNF signalling pathway with DELs' interactions and ASEs mapped, Figure S3: IL-17 signalling pathway with DELs' interactions, DASs and ASEs mapped, Figure S4: NOD-like receptor signalling pathway with DELs' interactions and ASEs mapped, Figure S5: NF-kappa B signalling pathway with DELs' interactions, DASs and ASEs mapped, Figure S6: MAPK signalling pathway with DELs' interactions and DASs mapped, Figure S7: Toll-like receptor signalling pathway with DELs' interactions mapped, Figure S8: Necroptosis with DELs' interactions and ASEs mapped, Figure S9: Jak-STAT signalling pathway with DELs' interactions and DASs mapped, Figure S10: Antigen processing and presentation with DELs' interactions and DASs mapped, Figure S11: Cell adhesion molecules (CAMs) with DELs' interactions and DASs mapped, Figure S12: Apoptosis with DELs' interactions and DASs mapped, Figure S13: Regulation of actin cytoskeleton with DELs' interactions, DASs and ASEs mapped, Figure S14: Leukocyte transendothelial migration with DELs' interactions, DASs and ASEs mapped, Figure S15: Cellular senescence with DELs' interactions mapped, Figure S16: VEGF signalling pathway with DELs' interactions mapped, Figure  S17: Proteasome with DELs' interactions mapped, Figure S18: Phagosome with DELs' interactions, DASs and ASEs mapped, Figure S19: Arachidonic acid metabolism with DELs' interactions mapped.

Acknowledgments:
The authors thank Dorota Makowczenko for assistance with the graphics processing.

Conflicts of Interest:
The authors declare no conflict of interest regarding the publication of this paper. The funders had no role in the design of the study; in the collection, analyses, or interpretation of data; in the writing of the manuscript, or in the decision to publish the results.