High levels of PIWI-interacting RNAs are present in the small RNA landscape of prostate epithelium from vitamin D clinical trial specimens.

BACKGROUND: Vitamin D, a hormone that acts through the nuclear vitamin D receptor (VDR), upregulates antitumorigenic microRNA in prostate epithelium. This may contribute to the lower levels of aggressive prostate cancer (PCa) observed in patients with high serum vitamin D. The small noncoding RNA (ncRNA) landscape includes many other RNA species that remain uncharacterized in prostate epithelium and their potential regulation by vitamin D is unknown.
METHODS: Laser capture microdissection (LCM) followed by small-RNA sequencing was used to identify ncRNAs in the prostate epithelium of tissues from a vitamin D-supplementation trial. VDR chromatin immunoprecipitation-sequencing was performed to identify vitamin D genomic targets in primary prostate epithelial cells.
RESULTS: Isolation of epithelium by LCM increased sample homogeneity and captured more diversity in ncRNA species compared with publicly available small-RNA sequencing data from benign whole prostate. An abundance of PIWI-interacting RNAs (piRNAs) was detected in normal prostate epithelium. The obligate binding partners of piRNAs, PIWI-like (PIWIL) proteins, were also detected in prostate epithelium. High prostatic vitamin D levels were associated with increased expression of piRNAs. VDR binding sites were located near several ncRNA biogenesis genes and genes regulating translation and differentiation.
CONCLUSIONS: Benign prostate epithelium expresses both piRNA and PIWIL proteins, suggesting that these small ncRNA may serve an unknown function in the prostate. Vitamin D may increase the expression of prostatic piRNAs. VDR binding sites in primary prostate epithelial cells are consistent with its reported antitumorigenic functions and a role in ncRNA biogenesis.

BACKGROUND

In the United States, 42% of adults have a vitamin D deficiency, a condition associated with increased morbidity from chronic diseases and cancers.1, 2, 3 Vitamin D is an essential hormone, which regulates gene expression through binding to the nuclear receptor and transcription factor vitamin D receptor (VDR).4 In addition to its well‐characterized systemic effects, vitamin D acts locally and VDR is ubiquitously expressed. Epidemiological studies have consistently shown an association between low vitamin D and increased aggressive pathology and mortality in prostate cancer (PCa) but not with PCa incidence.5, 6, 7, 8, 9 The recent Vitamin D and Omega‐3 Trail (VITAL) study did not show an effect of vitamin D on PCa incidence, but when the first 2 years of follow‐up were excluded, there was a significant 25% reduction in all‐cancer mortality in the vitamin D–supplemented patients.10 As the majority of PCa is indolent, the null relationship of vitamin D with PCa incidence is not as clinically relevant as its relationship with disease progression.

In the prostate, genomic targets of vitamin D have been suggested to decrease tumor aggressiveness by decreasing inflammation and proliferation and increasing apoptosis and differentiation.11 Studies of vitamin D in PCa cells in vitro and in animal models support antitumor activities.12, 13, 14, 15 Previous work from our lab and others have shown that vitamin D increases the expression of tumor‐suppressive microRNAs (miRNAs), small noncoding RNA (ncRNA) that regulate translation of messenger RNAs (mRNAs).16, 17 As miRNAs are generally downregulated in PCa, upregulation of miRNAs by vitamin D may play a role in its antitumorigenic effects in prostate.17 Very little is known about vitamin D regulation of other ncRNA species. A recent paper showed vitamin D‐regulation of melanoma‐associated long ncRNAs in keratinocytes.18 As roles are emerging for several non‐miRNA small ncRNAs in cancer, such as small nucleolar RNA (snoRNA), transfer RNA (tRNA) fragments, and PIWI‐interacting RNA (piRNA), it will be important to understand the effects of vitamin D on the expression of these molecules in the prostate.19

To investigate the normal small ncRNA landscape of the prostate epithelium, as well as the effect of high prostatic vitamin D concentrations on ncRNA expression, we isolated and sequenced small RNA from benign prostate epithelium from radical prostatectomy patients who participated in a vitamin D supplementation clinical trial. We examined vitamin D genomic regulation through VDR chromatin immunoprecipitation and sequencing in primary prostate epithelial cells. These studies revealed a previously unknown diverse population of small ncRNA in prostate epithelium and vitamin D‐regulation of piRNA expression.

MATERIALS AND METHODS

Patient specimens

Prostatectomy tissues were provided by Dr Theodorus Van der Kwast at the University of Toronto. The patients were participants in a double‐blind randomized clinical trial of vitamin D3 supplementation (Clinical trial CT00741364 http://www.clinical‐trials.gov) and deidentified.20 All patients in for sequencing were Caucasian and had a Gleason score of 3 + 3 or 3 + 4. For a 3 to 8 week period before surgery, patients in this study received 400 or 40 000 IU of liquid oral vitamin D3 daily. Prostatectomy tissue from the peripheral and transitional zones was freshly frozen and stored at −80°C. Serum 25D (cholecalciferol) and prostatic 1,25D (calcitriol) were measured by liquid chromatography‐tandem mass spectrometry and enzyme immunoassay, respectively, as part of the original trial (Table 1).20 Primary prostate epithelial cells (PrE) were isolated from deidentified benign prostatectomy tissue as determined by a pathologist and grown in PrEGM media (Lonza, Basel, Switzerland) for one passage on collagen‐coated plates17; the University of Illinois at Chicago (UIC) Institutional Review Board approved protocol #2011‐1138.

Table 1

Patient characteristics

Group	Patient sample ID	Age	BMI	Vitamin D status
				Serum 25D (nmol/L)	Prostatic 1,25D (pmol/kg)
Control	1	59	36.4	93	21.4
	2	57	22.6	62	22.9
	3	63	31.7	78	19.5
	4	50	24.3	81	21.6
	Average (StDev)	57.25 (5.34)	28.75 (6.45)	78.5 (12.76)	21.35 (1.40)
High Vitamin D	5	64	28.0	284	79
	6	64	26.6	336	40.2
	7	62	26.2	429	45.3
	8	58	24.5	431	35.4
	Average (StDev)	62 (2.82)	26.34 (1.44)	370 (72.47)	49.98 (19.76)

Abbreviations: BMI, body mass index; StDev, standard deviation.

Patients were enrolled in a clinical trial where they received either 400 IU (control) or 40 000 IU (high vitamin D) of vitamin D3 daily for 3 to 8 wk before radical prostatectomy. Vitamin D metabolites were measured as part of the original trial in prostatectomy tissues.20 Prostate 1,25D is reported here as the average measurement from the peripheral and transitional zones.

Laser capture microdissection RNA isolation and small‐RNA sequencing

Laser capture microdissection (LCM) RNA collection was performed as described.21 Approximately 1000 glands were collected from each patient. RNA was isolated using the RNAqueous‐Micro Kit (Ambion, Foster City, CA). Sequencing libraries were constructed using 500 ng of total RNA from each sample and the TruSeq Small RNA Sample Prep Kit (Illumina, San Diego, CA). Libraries were multiplexed and sequenced on one lane for 51 cycles on a HiSeq 2500 (Illumina) using TruSeq Rapid SBS sequencing chemistry v2 (Illumina) at the Core Genomics Facility (CGF) at UIC. Fastq files were generated with the bclfastq v1.88.4 (Illumina). PrE were treated with 50 nM 1,25D or ethanol for 24 hours before collection of total RNA. PrE libraries were constructed using 1 µg of total RNA and the TruSeq Small RNA Sample Prep Kit (Illumina). The libraries were run on one lane for 51 cycles on a HiSeq. 2500 (Illumina) in high output mode using TruSeq Rapid sequencing kit v3 (Illumina).

Adapter sequences were removed using cutadapt,22 and 3′ end trimming of bases with a Phred quality score less than 20. Reads less than 17 base pairs (bp) were discarded. For alignment, miRNA stem‐loop sequences were retrieved from miRBase (release 21),23 piRNA sequences from piRNABank,24 and ribosomal RNA (rRNA), small nuclear RNA (snRNA), snoRNA, and tRNA sequences from Ensembl version 90 (Hinxton, United Kingdom).25 Reads were aligned to small ncRNA sequences using SHRIMP 2.2.3 aligner in ungapped local alignment mode with seed size of 14 allowing two mismatches.26 The highest scoring match strata were reported. Counts were generated using a custom script so that reads mapping to multiple ncRNA species were designated in the hierarchical order: miRNA, rRNA, snoRNA, snRNA, tRNA, or piRNA. If one read sequence mapped to multiple of one species of small RNA, the read was split proportionally. Custom code is accessible at https://github.com/sgao30/smallRNAseq. Differential small RNA expression was calculated using EdgeR.27, 28 Small‐RNA sequencing and counts are deposited at NCBI GEO accession number GSE124589.

The Cancer Genome Atlas small‐RNA sequencing data

Access to The Cancer Genome Atlas (TCGA) data was obtained through dbGaP (project# 15240). We utilized small‐RNA sequencing files from the whole prostate in the prostate adenocarcinoma (PRAD) data set. We reanalyzed 11 benign small‐RNA sequencing samples from Roswell Park using our alignment pipeline. Adapter‐trimmed read sequences (both aligned and unaligned) were extracted using bedtools bam2fastq.29 Since TCGA small‐RNA sequencing reads are 15 to 30 bp, we only used reads 17 to 30 bp in our data set for comparisons.

piRNA sequence characterization

Consensus sequences of piRNAs with an average greater than 100 CPM (99 piRNAs and 92 piRNAs in control and high vitamin D groups, respectively) were created using the Berkley WebLogo Tool (https://weblogo.berkeley.edu/logo.cgi).30 The presence of a secondary biogenesis signature was assessed using the piPipes_smallRNA_intersect script from the piPipes package.31 Genomic locations (hg19) of piRNAs expressed greater than 100 CPM were determined using bedtools intersect29 with 1 bp overlap for the locations of piRNA sequences from piRBank,24 and UCSC genome annotations,32 RepeatMasker,33 or pooled/generic testis piRNA clusters from human testis in the piRNA cluster database (converted to hg19 with UCSC liftover).34, 35 piRNAs with multiple locations were designated in the following order: coding sequence, 3′‐untranslated region (3′‐UTR), 5′‐UTR, intron, piRNA cluster, and repeat DNA sequence.

Chromatin immunoprecipitation‐sequencing

PrE were grown in PrEGM media (Lonza) for one passage. Six 100 mm dishes per condition were treated for 2 hours with 50 nM 1,25‐vitamin D3 or ethyl alcohol (vehicle), cross‐linked in 1% formaldehyde for 10 minutes at reverse‐transcription (RT), then 125 mM glycine for 5 minutes. Plates were washed two times with ice‐cold phosphate‐buffered saline (PBS), incubated in TrypLE at 37°C for 3 minutes, and washed two times. Cells were lysed in Nuclear Isolation Buffer (Magna ChIP HiSens Chromatin Immunoprecipitation Kit, Millipore, Burlington, MA) on ice for 20 minutes with vortexing every 5 minutes. Nuclei were pelleted at 800g 5 minutes, then resuspended in 1 mL RIPA buffer with Protease Inhibitor Cocktail III (Millipore). Cells were sonicated on ice with a Sonic Dismembrator model 100 (Thermo Fisher Scientific, Waltham, MA). One milliliter (~1 µg) of sonicated DNA was used for chromatin immunoprecipitation (ChIP) and 10 µL was saved as input. Three hundred microlitre Dynabeads (Thermo Fisher Scientific) were washed three times and resuspended in 1000 µL of 1 mg/mL bovine serum albumin (BSA)/PBS/2 mM ethylenediaminetetraacetic acid and incubated with 24 µg of anti‐VDR C‐20x sc‐1008× (Santa Cruz Biotechnology, Dallas, TX) rabbit polyclonal for 5 hours, then washed three times. Antibody‐conjugated beads were incubated with sonicated DNA overnight rotating at 4°C. Beads were washed once with 5:1 RIPA:PBS containing 5 mg/mL BSA and Protease Inhibitor Cocktail III (Millipore), five times in 100 mM Tris‐HCl pH 7.5/ 500 mM LiCl/1% NP‐40 /1% sodium deoxycholate and once in TE buffer. DNA was eluted off the beads twice in 100 µL ChIP Elution Buffer (Millipore) containing 283 mM NaCl and 1 µL proteinase K and heated at 65°C for 1 hour with vortexing every 15 minutes. Eluted samples and input samples (added to 90 µL elution buffer) were incubated at 65°C overnight. DNA was purified using the Clontech NucleoSpin kit (Clontech, Mountain View, CA).

Libraries were prepared from 10 ng of ChIP DNA with the Library Construction Kit Kapa Biosystems (Wilmington, MA) with modifications: adapters were diluted 1:20 and adapter‐modified DNA were amplified for 10 cycles. The pooled libraries were sequenced on one lane each for 101 cycles on a HiSeq 2500 in high output mode using a TruSeq SBS sequencing kit version 3. Fastq files were generated and demultiplexed with the bcl2fastq v1.8.4 Conversion Software (Illumina).

ChIP‐seq validation and analysis

For polymerase chain reaction (PCR)‐validation, dynabeads conjugated to 8 µg of anti‐VDR C‐20x sc‐10008× (Santa Cruz Biotechnology) rabbit polyclonal, or 8 µg of anti‐Normal Rabbit IgG sc‐2027 (Santa Cruz Biotechnology) were incubated with 340 µL of the sonicated DNA as described above overnight at 4°C. In one condition, 40 µg of blocking peptide sc‐1008P (Santa Cruz Biotechnology) against anti‐VDR was added (a ratio of 5:1 blocking peptide:antibody). Bound DNA and inputs were eluted, reverse cross‐linked, and purified as above. Peak calling and differential peak analysis were performed using MACS236 with default parameters and EdgeR27, 28 at the Center for Research Informatics at UIC. In further analysis R statistical software v3.4.3 was used.37 A circos plot was generated with OmicCircos.38 Genes near peaks were identified with Bioconductor.39, 40 Gene set enrichment was performed with GSEA v.3.0,41 GO from MsigDB v6.142, 43 and KEGG v.6.1 pathways.44 Transcription factor enrichment was calculated using Homer v4.10,45 with default settings, including a range of 200 bp. Interactive Genomics Viewer (IGV) was used to visualize peaks.46 ChIP‐sequencing data and peak locations are deposited at NCBI GEO accession number GSE124576.

Immunostaining

Deidentified formalin‐fixed paraffin‐embedded (FFPE) tissue was from normal testis and benign prostate (UIC IRB protocol# 2012‐1033; UIC Biorepository (Chicago, IL). Immunofluorescence primary antibodies: PIWI‐like (PIWIL)‐1 (SAB4200365; Sigma, St Louis, MO) at 1:200 with Tris pH 9 heat‐induced epitope retrieval (HIER) 100°C, 10 minutes in a pressure cooker; KRT5 (Poly19055; BioLegend, San Diego, CA) at 1:1000; PIWIL‐2 (SAB3500749; Sigma) at 1:200 with sodium citrate pH 6 HIER 100°C, 5 minutes. Secondary antibodies: Alexa‐Fluor anti‐rabbit 555 and anti‐mouse 488 at 1:200; 4′,6‐diamidino‐2‐phenylindole (DAPI) nuclear counterstain. For immunohistochemistry, PIWIL‐4 (HPA057508; Sigma) at 1:300 with citric acid pH 6 HIER 125°C, 4 minutes. Rabbit IHC Kit (Abcam, Cambridge, UK) was used with hematoxylin counterstain. Imaged on the Vectra Automated Multispectral Imaging System (Perkin Elmer, Waltham, MA) and EVOS Cell Imaging System (Thermo Fisher Scientific). Immunocytochemistry primary antibodies (1:200‐1:300): PIWIL‐1 (15659‐1‐AP; Proteintech, Rosemont, IL); PIWIL‐2 (SAB3500749; Sigma), and PIWIL‐4 (HPA057508; Sigma); VDR (D‐6; sc‐13133; Santa Cruz Biotechnology). Secondary antibodies: 1:200 Alexa‐Fluor anti‐rabbit 488 or Alexa‐Fluor anti‐mouse 555; Phalloidin actin stain and DAPI nuclear counterstain. Imaged on the Cell Insight CX7 High‐Content Screening Platform (Thermo Fisher Scientific). For western, 30 µg protein lysate was used with the NuPage System (Invitrogen, Carlsbad, CA), and Odyssey CLx imager (LI‐COR, Lincoln, NE) with Odyssey Blocking Buffer. Primary antibodies at 1:1000: TDRKH (13528‐1‐AP; Proteintech); GAPDH (MAB 374; Chemicon, Temecula, CA); VDR (C‐20X, sc‐10008x and D‐6, sc‐13133; Santa Cruz Biotechnology); H3K4me3 (rabbit polyclonal; Millipore). Secondary antibodies at 1:10 000: LI‐COR anti‐rabbit 555 and anti‐mouse 488.

Tissues, cell culture, and quantitative RT‐PCR

Deidentified prostate tissue and deidentified normal testis tissue were obtained through the UIC Biorepository (UIC IRB 2011‐1138). PrE were grown in PrEGM media (Lonza), for one passage. Primary stromal cells (PrS) were isolated under the same IRB as PrE as described,47 grown in MCDB media (Sigma) for 1 to 3 passage. RWPE1, PC‐3, LNCaP, 22RV1, and LAPC4 were obtained through ATCC (Manassas, VA) and grown less than 20 passages. RWPE1 were grown in KSFM media (Gibco, Dublin, Ireland), PC‐3 in DMEM (Gibco) supplemented with 10% fetal bovine serum (FBS) (Gibco), LNCaP, 22RV1, and LAPC4 in RPMI (Gibco) with 10% FBS. Complementary DNA (cDNAs) were made with High‐Capacity cDNA RT Kit (Invitrogen). Primers are in Table 2. Quantitative PCR (qPCR) used SYBR green (Roche, Basel, Switzerland) and QuantStudio 6 (Thermo Fisher Scientific). PCR settings: 95°C 10 minutes (×1), 95°C 15 seconds, 58°C 30 seconds, 72°C 30 seconds (×40). Relative quantity was calculated by the −ΔΔC t method, using GAPDH or HPRT1 for normalization.

Table 2

List of primers used

Target	Strand	Sequences (5′‐3′)	Product length	Experiment
PIWIL1	F	CAGAATCCACTTCCTGGAACA	223	RT‐qPCR
	R	GAGCAGGAACACGAATGACA
PIWIL2	F	CCCAGGTTGTCAATGTTCG	278	RT‐qPCR
	R	CAGGCTGTCCACAATCTCC
PIWIL3	F	GCCACTAAAAGAGCGGAGAG	231	RT‐qPCR
	R	TCCAAACTGGTACCATGACG
PIWIL4	F	TCACCCAGACGTGTGAGAAG	249	RT‐qPCR
	R	GCTGAGCCTCACTGTTGTCA
HPRT1	F	TGCTGACCTGCTGGATTACA	208	RT‐qPCR
	R	CTGCATTGTTTTGCCAGTGT
GAPDH	F	AAGGTCGGAGTCAACGGATTTGGT	196	RT‐qPCR
	R	TGATGACAAGCTTCCCGTTCTCAG
LNA piR‐hsa‐9286	N/A	TCCTCTTGAGGTCTCTCTTCA	N/A	In situ hybridization
LNA piR‐hsa‐15063	N/A	AACACCGTCGTTCATTCATTCAT	N/A	In situ hybridization
CYP24A1/1	F	AGCACACCCGGTGAACTC	140	ChIP validation, known VDR binding site
	R	TGGAAGGAGGATGGAGTCAG
CYP24A1/2	F	TTCAAGAGGTCCCCAGACAC	203	ChIP validation, negative control
	R	AGTCGGGGCTTAACGATTCT
GAPDH	F	TACTAGCGGTTTTACGGGCG	166	ChIP validation, negative control
	R	TCGAACAGGAGGAGCAGAGAGCGA
TDRKH	F	CCAGGAGGCTGTGAAACTCAT	87	RT‐qPCR
	R	TCCACATCAATCCGAGCACC

Abbreviations: ChIP, chromatin immunoprecipitation; RT‐qPCR, reverse‐transcription quantitative polymerase chain reaction; VDR, vitamin D receptor.

In situ hybridization

The miRCURY LNA microRNA ISH Optimization Kit (Exiqon, Vedbaek, Denmark) was followed with modifications. FFPE sections (5 µm) of a deidentified prostate tissue microarray were baked at 60°C, deparaffinized, and incubated for 20 minutes at 37°C with 15 µg/mL Proteinase K. Digoxigenin (DIG)‐labeled piR‐hsa‐9286 and piR‐hsa‐15063 custom LNA‐probes (sequences in Table 2) were used 80 nM, incubated at 48°C for 60 minutes followed by stepwise 5 minutes washes in SSC at 42°C 1 × −0.5 × −0.2×, then SSC 0.2× at room temperature. Slides were blocked and incubated for 60 minutes with alkaline phosphatase (AP)‐conjugated anti‐DIG antibody (Roche) at 1:200. AP was visualized with Vector Red alkaline phosphatase substrate (Vector Laboratories, Burlingame, CA).

Statistical analysis

All statistics not yet mentioned were performed in GraphPad Prism 8.0.0. A value of P < 0.05 is considered significant.

RESULTS

LCM small‐RNA sequencing of prostate epithelium

The patients in the current study were from the control (400 IU/day) and high‐dose vitamin D (40 000 IU/day) groups of a clinical trial of cholecalciferol (25(OH)2D3) supplementation for 3 to 8 weeks before radical prostatectomy.20 This inactive prohormone is the form of vitamin D that is predominant in circulation and permitted patients to achieve supra‐physiological levels of prostatic vitamin D without adverse effects on clinical parameters of calcium or phosphate homeostasis.20 Characteristics of the patients, including vitamin D metabolite levels in serum and prostate tissue at the time of radical prostatectomy, are shown in Table 1. LCM was used to isolate benign epithelium from the prostatectomy tissues followed by total RNA isolation and sequencing.

Small‐RNA sequencing reads were trimmed and aligned to a reference library of annotated human pre‐miRNA, rRNA, snoRNA, snRNA, tRNA, and piRNA sequences. In total, 68% to 80% of reads per sample mapped to small ncRNA (Table SI), and in the small percentage of cases (<0.01%) where a read mapped equally well to multiple small ncRNA classes, the class was assigned hierarchically in the order of small ncRNA listed above. The majority of the mapped reads represented miRNAs, but surprisingly piRNAs, which are generally thought to be testis‐specific and re‐expressed in cancer,48 comprised approximately 20% of the mapped reads in the prostate epithelium (Figure 1A). The length distributions of reads were as anticipated for miRNAs (22 bp) and piRNAs (24‐32 bp) (Figure 1B). The longer rRNAs, snoRNAs, snRNAs, and tRNAs spanned a range of read lengths, suggesting detection of fragmented species. Of note, tRNA fragment lengths in the high vitamin D prostates compressed to a peak at the length of tRNA 5′‐halves (30‐34 bp), suggesting increased active tRNA processing (Figures 1B and S1). Prostatic vitamin D concentrations had a more pronounced effect on piRNAs than any other small ncRNA. Patients in the high vitamin D group expressed significantly more piRNAs (Figure 2A), and 22 individual piRNAs were significantly differentially expressed between the control and high vitamin D groups at an false discovery rate (FDR) of 20% (Figure 2B and Table SI). Seven miRNAs, including miR‐146a and miR‐595, and one snRNA were also differentially expressed (Table SIII).

Figure 1

Noncoding RNA in prostate epithelium. Small RNA in prostate epithelium was laser capture microdisssected and sequenced. A, Noncoding RNA composition of all mapped reads in each patient sample in the control (right) and high vitamin D (left) groups. B, Length distribution of reads mapped to small noncoding RNA species and unmapped reads in the control (right) and high vitamin D (left) groups (total reads = total reads after filtering low‐quality reads). High Vit. D, high vitamin D; miRNA, microRNA; piRNA, PIWI‐interacting RNA; rRNA, ribosomal RNA; snoRNA, small nucleolar RNA; snRNA, small nuclear RNA; tRNA, transfer RNA. [Color figure can be viewed at wileyonlinelibrary.com]

Figure 2

piRNAs in prostate epithelium and their regulation by vitamin D. A, The percent of total‐mapped reads mapping to piRNAs is higher in the high vitamin D group than the control group by a two‐tailed unpaired t test, *P = 0.037. B, The expression of the two most differentially expressed piRNAs in CPM‐mapped piRNA. C, Consensus sequence of piRNAs expressed at a level greater than 100 CPM (~100 piRNA) in control and high vitamin D groups. D, Secondary piRNA signature (prevalence of complementary 5′ 10 bp overlaps between reads) was assessed. Reads mapping to piRNAs and unmapped reads were included in the analysis for each sample. The average relative frequency of antisense 5′ ends at different lengths of overlap is shown for the control and high vitamin D groups. E, The locations of piRNA sequences in the genome for the piRNAs expressed greater than 100 CPM. One base pair overlap was required for annotation. The sequences from repetitive elements are further broken down in the bar graphs. CDS, coding sequence; CPM, counts per million; High Vit. D, high vitamin D group; piRNAs, PIWI‐interacting RNA; repetitive, repetitive element; rRNA, ribosomal RNA; tRNA, transfer RNA; UTR, untranslated region. [Color figure can be viewed at wileyonlinelibrary.com]

Characterization of piRNAs and PIWIL proteins in prostate

Due to the unexpectedly high expression of piRNAs in benign prostate and their association with vitamin D, we investigated prostatic piRNAs further. Although vitamin D increased overall piRNA abundance, it did not significantly alter characteristics of the most highly expressed piRNAs in prostate epithelium. These piRNAs display the characteristic 5′ preference for uracil, but no preference for adenosine at position 10, which is a signature of the secondary biogenesis pathway31 (Figure 2C). Analysis of the 5′‐complementarity of reads mapping to piRNAs and unmapped reads, representing potentially unannotated piRNAs, showed no selection for 10 bp offsets, suggesting the absence of secondary piRNAs (Figure 2D). Many piRNA sequences match several genomic locations, so the genomic context was assigned hierarchically giving preference to intragenic matches. Highly expressed piRNA sequences in prostate epithelium, primarily matched within introns and 5′‐UTR, rather than uniquely from piRNA clusters or intergenic transposable elements (Figure 2E). The majority of piRNAs originating uniquely from repetitive elements partially overlapped with tRNA genes or pseudogenes, not transposable elements as observed in gametes. However, most of the piRNAs within introns (74% control and 77% high vitamin D) lie within fragments of repetitive sequences within the intron, and/or repetitive sequence elsewhere in the genome.

We tested if the class of Argonaut proteins that bind piRNAs, PIWIL proteins, were also present in the prostate. In humans, there are four PIWIL proteins. PIWIL‐1, ‐2, ‐4 were present in prostate tissue by immunostaining (Figure 3A). PIWIL‐3 was not tested due to unsuitable antibodies. Staining for all PIWIL proteins was stronger in prostate epithelium than in stroma. Primary prostate epithelial cells (PrE) were found to express PIWIL‐1, ‐2, and ‐4 in nuclear and perinuclear regions (Figure 3B). Transcript expression of PIWIL proteins in PrE cells was detectable but lower than in normal human testis tissue (positive control) (Figure 3C). Fresh frozen benign whole prostate tissue expressed similar levels of PIWIL transcripts as PrE cells (Figure 3D). PIWIL2, PIWIL4, and PIWIL1 were expressed, but PIWIL3 was only detected at low levels in some samples. Unlike PrE cells, primary prostate stromal cells (PrS) expressed PIWIL2, PIWIL3, and PIWIL4, but not PIWIL1. Immortalized prostate cell lines on average expressed lower levels of PIWIL transcripts than primary cell lines, regardless if the immortalized cell line was benign (RWPE1) or from PCa (PC‐3, LNCaP, 22RV1, and LAPC4) (Figure 3D). These results suggest that PIWIL proteins are expressed in prostate tissue, and therefore, the piRNAs expressed in prostate tissue may be functional.

Figure 3

PIWIL expression in prostate tissues and cell lines. A, Human testis tissue (left) and benign human prostate tissue (right). First panel, H&E stain; second panel, immunofluorescence for PIWIL‐1 (green) on adjacent sections of tissue. Staining for KRT5 (red) was included to identify prostate basal epithelium. DAPI (blue) stains nuclei. Third panel, immunofluorescence for PIWIL‐2 (green) on adjacent sections of tissue. Forth panel, immunohistochemistry for PIWIL‐4 (brown) on the same tissues. B, PIWIL‐1, ‐2, and ‐4 (green) are detected in PrE cells. DAPI (blue) stains nuclei, and phalloidin (red) stains actin fibers. C, PIWIL transcript expression relative to GAPDH expression in PrE cells and normal testis tissue. GAPDH was used because it is expressed at similar levels in testis and prostate. Four technical replicates (error bars = standard deviation). D, PIWIL transcript expression relative to HPRT in benign prostate tissue (n = 1), PrE (n = 7), PrS (n = 5), and immortalized prostate cell lines RWPE1 (n = 6), PC‐3 (n = 6), LNCaP (n = 5), 22RV1 (n = 5), and LAPC4 (n = 3). Each biological replicate contained two technical replicates (error bars = standard error). DAPI, 4′,6‐diamidino‐2‐phenylindole; GAPDH, glyceraldehyde 3‐phosphate dehydrogenase; H&E, haemotoxylin and eosin; KRT5, keratin 5; PIWIL, PIWI‐like; PrE, prostate epithelial cell; PrS, primary stromal cell. [Color figure can be viewed at wileyonlinelibrary.com]

Comparison of small‐RNA sequencing of LCM epithelium to the whole prostate in TCGA

We compared the small ncRNA profiles from LCM benign prostate epithelium to the small‐RNA sequencing profiles of benign whole prostate available through The Cancer Genome Atlas Prostate Adenocarcinoma cohort (TCGA PRAD). This analysis served two purposes: to validate our findings and to evaluate the utility of LCM in minimizing contamination from the ubiquitous prostate stroma present in benign prostate. TCGA small‐RNA sequencing focused on miRNAs, but other small ncRNAs have been analyzed in TCGA data, including piRNAs.49, 50, 51 The TCGA tissue processing and sequencing platform were similar to those used in this study, however, TCGA selected for RNA fragments below 30 bp51 (Figure S2A), whereas this study selected for below 50 bp. Therefore, to facilitate comparison, only the 17 to 30 bp reads from the LCM‐collected tissue were compared with TCGA reads. TCGA reads from the 11 TCGA benign prostate specimens collected at Roswell Park were reanalyzed through our analysis pipeline.

Hierarchical clustering of LCM small ncRNA profiles with TCGA small ncRNA profiles showed a clear distinction between the groups. The LCM small ncRNA profiles were more highly internally correlated than the TCGA small ncRNA profiles, suggesting that LCM improved sample homogeneity by removing stromal bias (Figures 4A and 4B and S2B). The miRNA profiles of LCM prostate epithelium and whole prostate tissue were fairly concordant, with the same miRNAs highly expressed in both groups (Figure 4A). In contrast, there were marked differences in the profiles of another small ncRNA between the data sets (Figures 4B and S2B). This was apparent for piRNAs, suggesting that piRNA expression profiles may be strongly affected by the stromal bias in TCGA data (Figure 4B). To test for epithelial expression, we selected two piRNAs (<30 bp) which were highly expressed in our LCM samples (>1000 CPM‐mapped piRNA) and poorly expressed in TCGA samples (<10 CPM‐mapped piRNA) for in situ hybridization on a tissue microarray containing prostatectomy tissue from 32 patients. These piRNAs are only expressed from one genomic location; piR‐has‐15603 is unannotated and piR‐has‐9286 is intronic. Both piRNAs were detected in prostate and consistently showed higher expression in prostate epithelium than surrounding stroma (Figure 4C). These results demonstrate the value of performing LCM in studies of prostate tissue specimens.

Figure 4

Comparison to TCGA PRAD benign whole tissue small‐RNA sequencing data. Heatmaps of (A) miRNA and (B) piRNA expression in our LCM‐collected and in TCGA benign prostate whole tissue small‐RNA sequencing data. The log (CPM) for individual small RNAs within the samples (CPM reads mapped to each small RNA type) was calculated. Hierarchical clustering used Spearman correlation as the distance metric. Small RNA are ordered by average log (CPM) expression level in the LCM‐collected epithelium. (C) In situ hybridization for two piRNAs highly expressed in our LCM‐collected data and expressed at very low levels in TCGA data. H&E stain of an adjacent section is included for reference. No probe and scrambled probe controls are shown on the left. CPM, counts per million; DAPI, 4′,6‐diamidino‐2‐phenylindole; H&E, haemotoxylin and eosin; LCM, laser capture microdissection; miRNA, microRNA; piRNA, PIWI‐interacting RNA; TCGA, The Cancer Genome Atlas. [Color figure can be viewed at wileyonlinelibrary.com]

Comparison of small‐RNA sequencing of LCM epithelium to small‐RNA sequencing of PrE cells

Small‐RNA sequencing data was also collected from benign PrE that had been isolated from histologically benign areas from a radical prostatectomy tissue and grown in culture for a single passage then treated with 50 nM of 1,25D or vehicle for 24 hours. Compared with LCM‐collected prostate epithelium, PrE expressed predominately miRNAs and only a very small proportion of piRNAs (2.4% control, 3% vitamin D‐treated) (Figure 5A and 5B). This suggests that cell culture conditions cause dysregulation of small ncRNA in PrE. The most highly expressed ncRNA from prostate epithelium were compared individually to their expression in PrE cells. There was a significant correlation between LCM epithelium and PrE expression levels for all ncRNA, but the relationship was only weak to moderate for miRNAs (Figures 5C and S3). Of note, despite the decreased abundance of piRNAs in PrE cells, their expression in PrE was strongly correlated with LCM epithelium (Figure 5D). This finding, along with the presence of PIWIL proteins in PrE, supports the use of these cells as a model to study prostatic piRNAs.

Figure 5

Small noncoding RNA in prostate epithelial cells. PrE were treated with vehicle or vitamin D for 24 hours. A, Noncoding RNA composition of all mapped reads in control and vitamin D‐treated PrE. B, Length distribution of reads mapped to small noncoding RNA species and unmapped reads in PrE. C, The expression of the top expressed (avg. CPM > 100) miRNAs by Log(CPM) (Log2 CPM‐mapped miRNA) in control prostate epithelium was compared with their expression in Log(CPM) in control PrE cells. The Pearson correlation coefficient was calculated. D, The expression of the top expressed (avg. CPM > 100) piRNAs by Log(CPM) (Log2 CPM‐mapped piRNA) in control prostate epithelium was compared with their expression in Log(CPM) in control PrE cells. The Pearson correlation coefficient is shown. CPM, counts per million; miRNA, microRNA; piRNA, PIWI‐interacting RNA; PrE, prostate epithelial cell; rRNA, ribosomal RNA; snoRNA, small nucleolar RNA; snRNA, small nuclear RNA; tRNA, transfer RNA. [Color figure can be viewed at wileyonlinelibrary.com]

VDR‐ChIP‐sequencing in primary prostate epithelial cells

In tandem to the small‐RNA sequencing experiments, we investigated the genomic activity of vitamin D in the same PrE cells through VDR‐ChIP‐sequencing. PrE cells treated with 50 nM 1,25D showed nuclear localization the VDR within 1 hour, persisting until 4 hours, and a 2 hours treatment was selected for ChIP‐seq (Figure 6A). Rigorous testing of the ChIP conditions was performed. The VDR antibody efficiently enriched VDR by immunoprecipitation (Figure 6B). VDR ChIP‐PCR was used for validation and showed the interaction of VDR with the characterized VDR‐binding site in the CYP24A1 promoter, but not a CYP24A1 region ~300 bp away or the GAPDH promoter (Figure 6C). A blocking peptide toward the antibody blocked this interaction (Figure S4). The expected VDR binding site upstream of CYP24A1 is shown in Figure 6D. Nearly 5000 additional VDR binding sites were identified, with 690 significantly increased, and 807 significantly decreased VDR binding after 2 hours treatment with 50 nM 1,25D (Table SIV). Although more binding sites were significantly decreased following 1,25D treatment, the fold change (FC) of these sites (avg. log2(FC) = 1.48) was less than at sites with significant increases in VDR binding (avg. log2(FC) = 4.80).

Figure 6

VDR ChIP‐seq in PrE cells. VDR ChIP was performed in PrE treated with 50 nM 1,25D for 2 hours. A, Immunofluorescence for VDR in PrE cells treated for 0‐6 hours with 50 nM 1,25D. Nuclear localization was observed for 1 and 4 hours treatments. B, Immunoprecipitation of VDR in PrE cell lysates with VDR antibody and IgG control. GAPDH and H3M3K4 are used as negative and loading controls. C, VDR ChIP‐PCR in cells treated with 50 nM 1,25D for 2 hours shows antibody binding to the CYP24A1 VDR binding site (CYP24A1/1). A proximal region of CYP24A1 promoter (CYP24A1/2) and GAPDH are included as negative controls. D, Interactive Genome Viewer tracks for the read density of VDR ChIP (1,25D‐EtOH), RWPE1 DNaseI hypersensitive sites (ENCODE), H3K4me3 sites (NCBI GEO, GSE63094), and H3K27ac sites (NCBI GEO, GSE63094) across the CYP24A1 gene. E, A circos plot depicting chromosomal locations of significant VDR ChIP peaks as a heatmap. The inner rings show chromosomal locations of differentially expressed (P < 0.05) ncRNA from the small‐RNA sequencing data. ChIP, chromatin immunoprecipitation; EtOH, ethyl alcohol; GAPDH, glyceraldehyde 3‐phosphate dehydrogenase; IgG, Immunoglobulin G; ncRNA, noncoding RNA; PrE, prostate epithelial cell; VDR, vitamin D receptor. [Color figure can be viewed at wileyonlinelibrary.com]

The genomic localization and log2(FC) of VDR binding peaks are depicted in a circos plot (Figure 6E). VDR binding was predominately on autosomal chromosomes. The location of differentially expressed (pre‐FDR correction) small ncRNA from the LCM small‐RNA sequencing data are shown in the inner circles. There was no significant enrichment of VDR binding sites in proximity to these small ncRNA, suggesting a lack of direct regulation by VDR. However, several small ncRNA biogenesis genes were within 5 kb of a VDR‐ChIP peak that significantly changed following vitamin D treatment (Table SV). TDRKH (Tudor and KH domain‐containing protein) is a PIWIL adapter protein with an intronic VDR peak (Figure S5A). Four individual patient‐derived PrE cell lines showed elevated expression of TDRKH mRNA when treated with vitamin D (Figure S5B). An increase in TDRKH protein was observed in vitamin D‐treated RWPE1 prostate epithelial cells (Figure S5C). This suggests that TDRKH is a vitamin D‐regulated gene in prostate epithelium.

There was a high percentage of VDR‐ChIP peaks that showed decreased binding after 1,25D treatment. We compared these genomic locations to those where VDR binding increased. Gene set enrichment analysis using GO and Kegg terms were performed on genes with transcription start sites within 5 kb of a peak. Genes near peaks that were higher in control cells were enriched in pathways related to translation, oxidative phosphorylation, ncRNA processing, and p53 signaling (Figure 7A). Genes near peaks that were higher after 1,25D treatment were enriched in pathways related to morphogenesis and motility (Figure 7B). The central 200 bp of VDR‐ChIP peaks were evaluated for the enrichment of transcription factor binding sites. The consensus binding site for VDR was enriched under both conditions. In the peaks that were significantly higher in control cells, ETS transcription factor binding sites, such as Elk4 and ETS, were the most highly enriched (Figure 7C). Conversely, in peaks that were significantly higher under 1,25D treatment, bZIP transcription factor binding sites, such as Atf3 and AP‐1, were the most highly enriched (Figure 7D).

Figure 7

VDR binding site analysis. A, GSEA for GO and Kegg terms of genes with transcription start sites within 5 kb of a VDR‐ChIP peak that decreased after 1,25D treatment. B, GSEA for GO and Kegg terms of genes with transcription start sites within 5 kb of a VDR‐ChIP peak that increased after 1,25D treatment. C, Top 10 significantly enriched annotated transcription factor binding sites within 200 bp of the center of VDR‐ChIP peaks, in peaks that decreased after 1,25D treatment. D, Top 10 significantly enriched annotated transcription factor binding sites by HOMER within 200 bp of the center of VDR‐ChIP peaks, in peaks that increased after 1,25D treatment. ChIP, chromatin immunoprecipitation; GO, Gene Ontology; GSEA, Gene set enrichment analysis; Kegg, Kyoto Encyclopedia of Genes and Genomes; VDR, vitamin D receptor

DISCUSSION

We investigated the associations between vitamin D and ncRNA in prostate epithelium using small‐RNA sequencing of LCM‐collected tissue from a clinical trial of high‐dose vitamin D supplementation, and VDR ChIP‐sequencing of primary prostate epithelial cells. We previously reported that vitamin D upregulates miR‐100 and miR‐125b in prostate epithelium, which may contribute to its antitumor effects.17 In prostate stroma, vitamin D upregulated miRs, including miR‐21, and miR biogenesis gene DICER.47 Vitamin D has also been shown to also modulate miRNA expression in other tissues and serum.16 However, the effect of vitamin D on the majority of small ncRNA species is unknown. Our results suggest that piRNAs are expressed in normal benign prostate epithelium, and associated with vitamin D.

We first characterized the small RNA transcriptome of the normal prostate epithelium. Due to the high stromal content of prostate tissue, characterization of RNA in the prostate epithelium is challenging. Our study supports that using LCM to isolate prostate epithelium produces more consistent data than whole prostate sequencing, and more biologically relevant data than PrE cell sequencing. In benign prostate epithelium, we found an unexpectedly high expression of annotated piRNAs. PiRNAs were discovered and characterized in sperm, where they silence transposable elements.52, 53, 54 In comparison to miRNAs, piRNAs are longer (26‐31 bp), more numerous, and expressed as single‐stranded precursors independently of Dicer.54 By definition, piRNAs interact with PIWIL proteins, which we found in prostate.52, 53 Recent small‐RNA sequencing studies of human somatic tissues have also detected hundreds of annotated or predicted piRNAs in ovary, liver, brain, and cell lines.55, 56, 57, 58, 59 The piRNAs identified in other human tissues share characteristics with the piRNAs we identified in prostate epithelium, such as a large portion mapping to introns and a weak or absent secondary piRNA biogenesis signature. In contrast, piRNAs in the germline are expressed from piRNA clusters and transposable elements and have a strong secondary biogenesis signature.54 This may suggest that somatic piRNAs have different functions than in the germline.

In comparison to control patients, prostate epithelium from patients with high concentrations of vitamin D expressed significantly more piRNAs. Overall abundance of other small ncRNA did not change, although the distribution of tRNA read lengths was altered to suggest tRNA halves, or tiRNAs. tiRNAs may reduce translational activity, which is consistent with reported negative effects of vitamin D on proliferation.60 Furthermore, tiRNAs are reported to be dysregulated in cancer and may have functional roles.61 The gene for Angiogenin (Ang), which produces tiRNA from tRNA, was located within 5 kb of a VDR‐ChIP peak in PrE cells, suggesting potential VDR regulation of this gene. Interestingly, tiRNAs have been shown to interact with PIWIL‐4 in a human somatic cell line.62 Prostates with high vitamin D also showed differential expression of several individual ncRNAs. Nearly two dozen piRNAs were differentially expressed seven miRNAs, and one snRNA. The miRNAs included anti‐inflammatory miR‐146a, which in this study was negatively associated with vitamin D, in contrast to studies in colon tissue and plasma.16 The differentially expressed miRNAs also did not include another miRNA previously identified as regulated by vitamin D.16 This inconsistency may be due to limitations of this study, including insufficient sample size, and the absence of patients meeting the clinical threshold of deficiency in vitamin D (<75 nmol/L serum 25D). All control patients receiving 400 IU of vitamin D had clinically sufficient vitamin D serum levels. As most of the United States population is deficient in vitamin D, it will be important to validate the regulation of ncRNA by vitamin D in a deficient setting.1

The potential function of piRNAs in vitamin D signaling is not known, but intriguingly the piRNA/PIWIL system promotes differentiation of somatic tissues in nonmammalian model organisms.63 It is tempting to hypothesize that piRNAs play a role in vitamin D‐induced differentiation of prostate epithelium, which is considered to contribute to its antitumor effects. In contrast, a large discrepancy exists between vitamin d‐ and reported piRNA‐associated functions, as studies indicate that piRNAS and PIWIL proteins promote tumorigenesis. Elevated expression of piRNAs and PIWIL proteins has been reported in several cancers, including increased expression of PIWIL‐2 in PCa.48, 64 A large reanalysis of TCGA small‐RNA sequencing for piRNA expression demonstrated that the expression profiles of piRNAs can differentiate normal and tumor tissues by unsupervised hierarchical clustering.50 However, while the association of higher piRNA expression in cancer observed in this study is supported by many others,48 differences in our LCM data suggest that analysis of piRNAs in TCGA PRAD should be interpreted with caution. Very few studies have focused on piRNAs in PCa specifically. Consistent with a protective function, studies in LNCaP PCa cells have shown that PIWIL‐1 overexpression is protective against genomic rearrangement of the TMPRSS:ERG gene caused by prolonged testosterone exposure.65 In predicting whether vitamin D‐induced changes in prostatic piRNAs are antitumorigenic, it is important to consider that some studies suggest that both low and high serum levels of vitamin D increase the risk of developing PCa.66 In the clinical trial that our LCM tissues were obtained from, there was decreased Ki67 proliferation marker expression in the prostates with higher vitamin D, suggesting reduced tumor proliferation.20 However, serum vitamin D levels similar to those in the high vitamin D group have not been investigated in epidemiological studies of PCa, as they are not common in natural populations.

We investigated gene regulation by vitamin D through VDR ChIP‐sequencing in PrE cells. Several genes involved in ncRNA biogenesis and function were located within 5 kb of a VDR peak that was significantly altered under 50 nM 1,25D conditions. VDR binding decreased near RNA polymerase II and III components, suggesting overall decreased transcription that likely affects both mRNA and ncRNA expression. In addition, many aminoacyl tRNA synthetase genes were near peaks that decreased after vitamin D treatment. This is consistent with the phenotype of decreased proliferation observed in vitamin D‐treated prostate cell lines.12 Unexpectedly, a VDR binding site near miRNA processing gene DCGR8 was decreased following vitamin D exposure. This could suggest decreased miRNA processing, in contrast to our previous finding of vitamin D upregulation of miRNAs in PrE cells under similar conditions.17 Patient variability is a major limitation here since only one primary cell line was used for ChIP‐sequencing and small‐RNA sequencing, and this patient did show lower miRNA expression following vitamin D treatment. In addition, the functional outcome of VDR binding in all these cases is not known. We validated one piRNA biogenesis gene, TDRKH, as upregulated by vitamin D in prostate epithelial cell lines. In mouse spermatocytes, the TDRKH homolog binds to PIWIL‐1 and ‐4 homologs and promotes primary piRNA biogenesis, PIWIL‐4 nuclear localization, and LINE1 DNA methylation.67 It is unclear if TDRKH contributes to the differential expression of piRNAs in high vitamin D prostates.

Other interesting findings arose from the VDR ChIP‐sequencing data. Under 1,25D‐stimulated conditions, VDR bound near genes enriched in morphogenesis and differentiation pathways. In the adult prostate, vitamin D may act in concert with an androgen to promote the secretory phenotype of glandular cells,68 which is lost in PCa. Intriguingly, in stimulated PrE, VDR bound near genes regulating neural differentiation. This finding, along with vitamin D‐downregulation of interleukin 6 in prostate,14 may be relevant to neuroendocrine transdifferentiation, a development in PCa that is associated with poor prognosis.69 Transcription factor site enrichment in VDR peaks was very different between unstimulated and stimulated peaks. Unstimulated peaks were highly enriched for ETS transcription factor binding sites. Activation of ETS genes through gene fusions, such as TMPRSS2:ERG, is common in PCa and drives proliferation.70 Negative cross‐talk between vitamin D and ETS transcription factors has been previously reported; the TMPRSS2:ERG fusion increases CYP24A1 expression thereby reducing cellular 1,25D levels, and vitamin D treatment of VCaP PCa cells negatively regulated gene sets associated with ERG signaling.71 However, others have shown that Ets‐1 can also directly bind VDR and promote its ligand‐independent activation.72 Consistent with these studies, our data suggest that in PrE cells, unliganded VDR colocalizes with ERG transcription factors. Upon 1,25D stimulation, VDR binding at these sites decreases, which may negatively affect ERG‐mediated gene expression.

CONCLUSIONS

In summary, prostate epithelium contains a diverse population of small ncRNA, and an unexpected abundance of annotated piRNAs and their partner PIWIL proteins. Supplementation with high‐dose calcitriol in PCa patients is associated with changes in the expression profiles of prostatic ncRNAs, including a significant upregulation of piRNAs. VDR binding sites in PrE cells were located near genes consistent with reported functions of vitamin D in prostate proliferation, apoptosis, and differentiation. Vitamin D may also regulate ncRNA expression through the many ncRNA biogenesis genes located near VDR‐ChIP binding peaks. Limitations of this study include the low sample numbers and the lack of data from vitamin D‐deficient patients. Analysis of annotated ncRNAs also precluded identification of novel ncRNAs from this data set. Future studies are needed to characterize a functional piRNA/PIWIL system in the prostate, and to determine if the functions of vitamin D/VDR in prostate epithelium identified in this study contribute to its antitumorigenic effects in the prostate.

CONFLICT OF INTERESTS

The authors declare that there are no conflict of interests.

Supporting information

Click here for additional data file.

Supporting information

Click here for additional data file.