Microarray data and DGE profiling
Two gene expression profiles (GSE35988, GSE32269) were obtained from Gene Expression Omnibus (GEO, https://www.ncbi.nlm.nih.gov/geo). GSE35988 contained 27 metastatic CPRC samples and 47 localized PCa samples, and GSE32269 was composed of 29 metastatic CPRC samples and 22 primary PCa samples.
Animals and castration procedure
Twenty healthy and mature male Sprague-Dawley (SD) rats (12 weeks old and the body weight reached 350 ~ 400 g) were obtained from the Animal Research Centre of the Lv Ye Pharmaceutical Company (Yantai, China) and were housed at an environmental temperature of 20 ± 2°C, receiving 12 h of light and 12 h of dark, with food and water provided ad libitum. Rats were randomly divided into two groups by using a random number table (10 animals of each group), namely sham-operated control group (Con-EP) and castrated group (Cas-EP). Before castration treatment, rats were anesthetized by intraperitoneally injection of sodium pentobarbital (30 mg/kg body weight). Thereafter, the efferent ducts and the testicular arteries and veins were ligated, then the testes were removed from the ligation point, leaving the intact epididymides in the scrotal area. In the meantime, the sham operation was also taken for the control group. Analgesia was provided by administration of meloxicam Q 24hrs (1 mg/kg PO or SC) (Metacam, Boehringer Ingelheim, St Joseph, MO) the day prior, the day of and the day following surgery. No obvious behavioral changes were observed in both control and the operation group. At the conclusion of the experiment, all the surviving rats were intraperitoneally anesthetized with sodium pentobarbital (30 mg/kg body weight) and their necks were dislocated for euthanasia. The whole epididymides were taken out from the castrated and sham-operated control rats after 7 days’ recovery and were immediately frozen in liquid nitrogen and stored at -80 ℃ before use. Protocols for the use of animals in these experiments were approved by the Research Animal Care and Use Committee of Yantai University (the approval reference number: YD.No20190916S0350210[315]) and were carried out in strict accordance with Guidelines for Ethical Review of Laboratory Animal Welfare of China. Ethical approval was obtained before the initiation of the research work and all efforts were made to minimize suffering.
Library preparation and sequencing
Total RNA was isolated from the epididymides from the control group and castrated group, respectively, by using RNAiso Plus reagent (TaKaRa, Dalian, China). The purity and quality of total RNA were evaluated by Agilent 2000 and were considered sufficient for subsequent library construction and sequencing. Tag libraries for the samples were prepared using the Illumina gene expression sample, Prep Kit. Each sample of 6 µg of the total RNA was extracted and purified by oligo (dT) magnetic bead adsorption. Primed by oligo (dT), mRNA bound to the bead was used as the template to synthesize a double-stranded cDNA. The 5' ends of tags were digested with restriction enzyme NlaIII at CATG sites. The cDNA fragments with 3' ends connected to Oligo(dT) beads were purified and the Illumina adaptor 1 was ligated to the 5' ends. The junction of Illumina adaptor 1 and CATG site is the recognition site of MmeI, which is a type of endonuclease with separated recognition sites and digestion sites. It cuts at 17 bp downstream of the CATG site, producing tags with adaptor 1. After removing 3' fragments with magnetic beads precipitation, Illumina adaptor 2 is ligated to the 3' ends of tags, acquiring tags with different adaptors of both ends to form a tag library. After 15 cycles of linear PCR amplification, 105 bp fragments were purified by 6% TBE PAGE Gel electrophoresis. After denaturation, the single-chain molecules were fixed onto the Illumina Sequencing Chip (flowcell). Each molecule was grown into a single-molecule cluster sequencing template through Situ amplification. Then, they were sequenced with the method of sequencing by synthesis (SBS). Each tunnel generated millions of raw reads with a sequencing length of 49 bp.
Tag annotation and gene expression level quantification
Five steps were needed to transform raw tags into clean tags: (1) removal of the 3' adaptor to preserve the 21 nt long sequences; (2) removal of the empty reads (reads with only adaptor sequences but no tags); (3) removal of the low-quality tags (tags with unknown sequences); (4) removal of too long or too short tags, keeping the 21 nt tags; (5) removal of tags with only one copy number (probably because of sequencing error). To convert digital profiles to gene expression, the tag sequences were mapped to the reference genes of rat from NCBI, with no more than 1 nucleotide mismatch. Clean tags mapped to multiple loci were filtered, only unambiguous mapped tags were kept. To identify gene expression, the number of unambiguous tags for each gene was calculated and then normalized to the number of transcripts per million tags (TPM). Gene ontology was assigned to all mapped genes to annotate their possible functions.
Real-time quantitative PCR analysis of genes
Real-time quantitative PCR (qPCR) was performed to examine the relative gene expression levels. Briefly, by using ReverTra Ace reverse transcriptase (Toyobo Co., Osaka, Japan), 1 ng~1 μg total RNA was reverse-transcribed into cDNA with oligo(dT)18, and then qPCR was performed on the QIAGEN’s Roter-Gene Q instrument by using Platinum SYBR Green qPCR SuperMix-UDG (Life Technologies, Cat. no.: 11733-038, CA, USA). In each reaction, 20 μL reaction mixtures containing 1 μL cDNA were incubated at 95℃ for 5 min, followed by 40 cycles of 95℃ for 10 s and 60℃ for 45 s. GAPDH was taken as an endogenous reference. 2-△△Ct method was used to calculate the differences of the expression level of genes in samples examined [11]. All experiments were run in triplicate and results were shown as the mean±SD (n=10). The primer sequences were available in Table 1.
Table 1
Primer sequences of qPCR
Gene
|
Primer sequence (5’-3’)
|
Nfkbie
|
F: GTGGACTGGATGGAGATTCTTG
R: TTTCCTGGTGGCTGGTAATG
|
Plscr4
|
F: CAGCTTGGGACACTAGGTTATT
R: GGGAACTAAGGGCGTCATTT
|
Apoc1
|
F: ACAAGGACAGGGTAGAGAAGA
R: ACAGGAAGTGCGATGAAGAG
|
Lcn8
|
F: GGGTAGAAGGCTTGTTCCTTAC
R: CTCTTTCTGAACCCACTGATCTT
|
Prdx1
|
F: TGTGTCCCACGGAGATCATTGCTT
R: TGTTCATGGGTCCCAATCCTCCTT
|
Prdx3
|
F: TGGTGTCATCAAGCACCTGAGTGT
R: AAGCTGTTGGACTTGGCTTGATCG
|
GAPDH
|
F: TGGGTGTGAACCACGAGAA
R: GGCATGGACTGTGGTCATGA
|
Identification of DEGs of digital profiling data
To identify the differentially expressed genes in the epididymis of rats after bilateral castration, a modified algorithm according to Audic S, et al. was used [12]. Suppose the number of clean tags corresponding to gene A was x, the expression level of each gene was only a small fraction of all genes expression level, therefore, p(x) follows the Poisson distribution. Furthermore, suppose the number of clean tags of sample 1 and 2 was N1 and N2, respectively, while the clean tags of gene A in sample 1 and 2 were x and y, respectively. The probability of gene A expression was equal in two samples could be calculated as the formula:
The acquired p-value was controlled with false discovery rate (FDR) using Benjamini-Hochberg method [13]. Differentially expressed genes with FDR≤0.001 and log2 fold-change≥1 were extracted. The correlation between the two libraries was measured by the Pearson correlation coefficient.
Identification of DEGs of microarray data
Online analysis tool GEO2R (https://www.ncbi.nlm.nih.gov/geo/geo2r) and R software (v4.0.5; https://www.r-project.org) were used for the identification of DEGs. The raw data from GSE35988 on GPL6480 and GSE32269 on GPL96 were normalized, transformed into expression values with GEO2R. Then empirical Bayes methods were applied to identify DEGs between CRPC and primary PCa with the limma package [14]. DEGs were defined with false discovery rate (FDR)<0.001 and |log2(fold change)|>1[15].
Functional enrichment and PPI network analysis of DEGs
The Database for Annotation, Visualization, and Integrated Discovery (DAVID; http://www.david.niaid.nih.gov) [16] and the R package clusterProfiler were used to perform GO and KEGG pathway analyses separately on CRPC-castration and CRPC-specific DEGs. Search Tool for the Retrieval of Interacting Genes/ Proteins (STRING; https://string-db.org) was used to construct PPI networks [17, 18]. Protein interactions with combined scores of > 0.15 were selected for the PPI network construction. Further, the Cytoscape software (v3.8.2) was utilized to calculate the node degree by MCODE and CytoHubba apps [19].