Comparison of miRNA expression profiles in pituitary–adrenal axis between Beagle and Chinese Field dogs after chronic stress exposure

Ethics statements

In this experiment, all animals used were approved by the Animal Care Committee of Jinan University (Guangzhou, People’s Republic of China) approval number 20131018001, and strictly implemented in line with the experimental basic principles.

Chronic stress treatment

Chronic stress treatment was given briefly as follows: 6 unrelated purebred male Chinese Field dogs (CFD) and 6 unrelated purebred male beagles, similar in health, weight, and other aspects were selected randomly. From each of these two breeds, three dogs were selected for stress exposure via intermittent electrical stimulation and restraint stress, and the other 3 dogs of each breed, used as controls, were not exposed to the stress. Each morning for 10 days, dogs were restrained and electrically stimulated with a stable current of 10 mA for 6 s followed by a 6-s interval without stimulation, lasting for 20 min every day. Meanwhile, before and after each 20-min stress exposure session, 4 ml of blood was collected and isolated to examine cortisol level with a cortisol radiation immunoassay kit used.

Samples and RNA preparation

All 12 dogs were killed by air embolism on the 11th day, along with all parts of the adrenal cortex and pituitary tissues collected and fast frozen in liquid nitrogen until for miRNA-seq. Further, Trizol reagent (Invitrogen, USA) was used to isolate the RNAs of all the adrenal cortex and pituitary tissues collected above. The quantity and quality of RNA were examined by Agarose gel electrophoresis, NanoDrop 2000 (Thermo, USA), and Agilent 2100 (Agilent, USA). The RNA collected from each group, including 3 dogs of each, were pooled with equal mass respectively, thus obtaining 8 RNA pooled samples, including CFDP1 (CFD pituitary with stress exposure), CFDP2 (CFD pituitary with non-disposal), CFDAC1 (CFD adrenal cortex with stress exposure), CFDAC2 (CFD adrenal cortex with non-disposal), BP1 (beagle pituitary with stress exposure), BP2 (beagle pituitary with non-disposal), BAC1 (beagle adrenal cortex with stress exposure), and BAC2 (beagle adrenal cortex with non-disposal). In addition, tissues from the hippocampus were collected to perform slice analysis (Nissl staining), aiming to evaluate the influences of stress exposure on hippocampal region cells.

Small RNA library construction and sequencing

A total of 1 µg RNA was collected from each RNA pooled sample above. Then, specific 3′ and 5′ adaptors were added to the RNA ends of each group (Truseq TM Small RNA sample prep Kit, Illumina). After that, the RNA with adaptors were reverse transcribed into its 1st cDNA, with random primers added. PCR was performed to amplify the cDNA with 11–12 cycles run and then resulted in 8 cDNA libraries, including CFDP1, CFDP2, CFDAC1, CFDAC2, BP1, BP2, BAC1 and BAC2. After purified further by polyacrylamide gel electrophoresis (PAGE), each cDNA library was performed by 1*50 bp high-sequencing on Hiseq-2000 (Majorbio Inc, Shanghai, China).

Sequence data analysis

The raw reads obtained by high-sequencing were further performed with quality control to result in clean reads. In brief, the quality control was that of trimming low-quality reads (ambiguous N and length < 18 nt), 3′ adapters, 5′ adapters and poly (A) sequences, with SeqPrep (https://github.com/jstjohn/SeqPrep) and Sickle (https://github.com/najoshi/sickle) used. Thereby, the lengths of 18–32 bp reads, as the clean ones, were obtained and calculated for their length distributions, using Fastx-Toolkit (http://hannonlab.cshl.edu/fastx˙toolkit/). Meanwhile, the identical reads were collapsed to obtain the unique reads for the statistics of the small RNA’s species and abundance. Moreover, the clean reads were mapped to Rfam database (ftp://sanger.ac.uk/pub/databases/Rfam/) and GenBank noncoding RNA database (http://blast.ncbi.nlm.nih.gov/) to annotate the miscellaneous RNAs, thus filtering out the rRNA, scRNA, snoRNA, snRNA, tRNA and other non-coding RNA. The remaining sRNAs obtained above, including known miRNAs and unknown miRNAs, were then mapped to the CanFam 3.1 genome to analyze their distributions on genome and calculate their expressive quantity. In addition, the remaining sRNAs were mapped to the Canis Canine miRNA data of miRBase 21.0 to identify the known miRNAs. Meanwhile, the novel miRNAs were predicted using miRDeep2 (http://www.mdcberlin.de/en/research/research˙teams/systems˙biology˙of˙gene˙regulatory˙elements/projects/miRDeep/index.html) from unannotated sRNAs (Bonnet et al., 2004; Langmead et al., 2009; Friedländer et al., 2012). The data discussed here have been deposited in NCBI Sequence Read Archive and is accessible through GEO Series accession no. GSE72015 (http://www.ncbi.nlm.nih.gov/geo/query/acc.cgi?acc=GSE72015).

Analysis of differentially expressed miRNAs

To explore the differentially expressed miRNAs (DE-miRNAs) in pituitary and adrenal, on the basis of miRNA expression profiles, eight comparisons were made, including BP1_vs_BP2, CFDP1_vs_CFDP2, BP1_vs_CFDP1 and BP2_vs_CFDP2 in pituitary; BAC1_vs_BAC2, CFDAC1_vs_CFDAC2, BAC1_vs_CFDAC1 and BAC2_vs_CFDAC2 in adrenal cortex. Here, R package DEGseq using a Benjamini q-value of 0.001 as the cut-off (corrected p-value < 0.001) (http://www.bioconductor.org/packages/release/bioc/html/DEGseq.html) (Wang et al., 2010), combined with Value1/Value2 > 2 (Value indicated the normalized count of samples), were used to screen out the DE-miRNAs. Furthermore, the common miRNAs that differentially expressed in at least two comparisons from pituitary and adrenal cortex were selected for the following target prediction, respectively.

Target gene prediction of DE-miRNA and its GO, KEGG analysis

For the target gene prediction, we utilized Targetscan 7.0 (http://www.targetscan.org/vert˙70/). Importantly, target genes predicted with the cumulative weighted context + + score less than −0.2, were selected for further analysis. Then, the target genes were selected to perform GO and KEGG pathway analysis, using the KOBAS 2.0 Functional Annotation Tool (http://kobas.cbi.pku.edu.cn/program.inputForm.do?program=Annotate) (Corrected p-value < 0.05).

Conjoint analysis of DE-miRNAs and DEGs

To better identify the important miRNAs in regulating chronic stress response, the DEGs identified in our previous transcriptome study of the pituitary and adrenal cortex were compared with the target genes of DE-miRNAs in the pituitary and adrenal cortex separately, thereby obtaining the genes that were both DEGs and target genes in pituitary and adrenal cortex, respectively (as a matter of convenience, the genes that were both DEGs and target genes were named “DE-target genes”). Further, according to a negative correlation between miRNAs and their target genes in expression pattern, the DE-miRNAs and their corresponding DE-target genes that went against this pattern, were then removed, thus resulting in the final DE-miRNA and its corresponding DE-target gene (named DE-(miRNA-mRNA) pairs hereafter). Subsequently, combined with the expression level of DE-miRNAs and their corresponding DE-target genes’ relatedness to stress regulation, a number of DE-(miRNA-mRNA) pairs were selected to perform dual luciferse reporter assay.

Validation of DE-miRNAs by qPCR

To elucidate the validity of the miRNA-seq data, DE-miRNAs, including miR-30a, miR-124, miR-205 and miR-222, were further detected by qPCR. In each group, 1 µg of pooled RNA was reverse transcribed using ReverTra Ace qPCR RT Kit 101 (TOYOBO, Japan), along with bulge-loop RT primers added. The bulge-loop RT primers and primers for qPCR were both synthesized by RIBOBIO (Guangzhou, China). QPCR was performed on the Bio-Rad S1000 with Bestar SYBR Green RT-PCR Master Mix (DBI Bioscience, Germany). Besides, U6 and U48, as the reference genes, were used to normalize the miRNA expression on the basis of 2^−ΔΔCT method (Schmittgen & Livak, 2008). The detailed information of the primers and the annealing temperatures was presented in Table S2.

High-sequencing of small RNA

In this study, eight libraries including BAC2, BAC1, CFDAC2, CFDAC1, BP2, BP1, CFDP2 and CFDP1 were constructed for miRNA-seq. As a result, a number of raw reads ranging from 15,926,611 to 81,028,035 were obtained from these 8 groups (Table 1). After eliminating reads of low quality (ambiguous N and length < 18 nt), 3′ adapters, 5′ adapters and poly (A) sequences, a number of clean reads ranging from 11,900,933 to 74,445,014 were obtained from the 8 groups (Table 1). Further, the clean reads from each group were mapped to Rfam database (11.0, http://Rfam.sanger.ac.uk/) and GenBank noncoding RNA database (http://blast.ncbi.nlm.nih.gov/) to annotate miscellaneous RNAs. And the detailed annotations of small RNA in each library were presented in Fig. S1.

The statistics of size distributions showed that the majority of sRNAs were 21–23 nt in length across 8 samples, along with a significantly downregulated profile in post-stress group compared with that in pre-stress (Figs. 1A and 1B). Besides, when mapped to the reference genome of Canis Canine with the remaining sRNA reads (including known miRNAs and unknown miRNAs), Chromosomes 1, 3, 5, 20 and 31 were found to be mapped most abundantly (ratio > 5%) in dog (Fig. 2). Meanwhile, a bigger number of small RNA reads were showed in pituitary samples (especially in CFD pituitary samples) compared with that in adrenal cortex samples.

Identification of known miRNAs

When mapped to the Canis Canine data of miRBase 21.0 by the remaining sRNA reads, a total of 425 out of 453 known miRNAs (Table S3) were identified across eight samples. Of the 425 identified known miRNAs, BAC2, BAC1, CFDAC2, CFDAC1, BP2, BP1, CFDP2 and CFDP1 contributed to 332, 325, 332, 345, 327, 318, 352 and 355 miRNAs, respectively.

Besides, the importance attached to the abundance of the miRNA expression in these two breeds was conducive to identify the miRNAs of high activity in dog pituitary and dog adrenal cortex. In adrenal cortex samples, cfa-miR-99a was the most abundant known miRNA, followed with cfa-miR-30 family (including cfa-miR-30a, cfa-miR-30b, cfa-miR-30c and cfa-miR-30d), cfa-miR-26 family (including cfa-miR-26a and cfa-miR-26b) and cfa-miR-7 family (including cfa-miR-7 and cfa-miR-7g) (The top 20 known miRNAs of abundance in adrenal cortex were listed in Table S4). Meanwhile, cfa-miR-7 family (including cfa-miR-7 and cfa-miR-7g), cfa-miR-99a, cfa-miR-30 family (including cfa-miR-30a and cfa-miR-30d), and cfa-miR-125 family (including cfa-miR-125a and cfa-miR-125b) were found to be the most abundant known miRNAs in pituitary samples (The top 20 known miRNAs of abundance in pituitary were listed in Table S4). In previous high-throughput sequencing studies of other species, various isoforms of the miRNAs were always detected (Li et al., 2011; Visser et al., 2014). In this study, we found almost all known cfa-miRNAs were of isoforms, with a larger number of isoforms possessed by the miRNAs of higher abundance. Moreover, most of the isoforms showing the highest expression level were identical to the canonical forms in miRBase 21.0, but for some other miRNAs, like the cfa-miR-217, the isoforms that showed the highest expression level were not the canonical ones at all (Fig. S2), indicating the variability of miRNA expression pattern.

Identification of novel miRNAs

In this study, a total of 147 novel miRNAs were identified across eight samples (Table S5). Compared with the known miRNAs identified here, rather lower expression levels were showed by the identified novel miRNAs. Of the 147 novel miRNAs, there were just 18 that were greater than 1,000 reads in expression level (Table S6). MiR-24_20231 and miR-21_17074 as the top two of the 18 novel miRNAs above, their predicted secondary structures were presented in Fig. S3. Of the 147 novel miRNAs, 75 shared a seed sequence with other known miRNAs of other species (especially mammals).

Identification of differentially expressed known miRNAs

In this study, a total of 27 (6 up, 21 down), 37 (18 up, 19 down), 30 (22 up, 8 down), 21 (15 up, 6 down), 25 (19 up, 6 down), 58 (13 up, 45 down), 32 (7 up, 25 down) and 43 (32 up, 11 down) DE-miRNAs (P-value < 0.001, log₂(Fold_change) > 1) were detected in the comparisons of BAC1_vs_BAC2, BAC1_vs_CFDAC1, BAC2_ vs_CFDAC2, CFDAC1_vs_CFDAC2, BP1_vs_BP2, BP1_vs_CFDP1, BP2_vs_CFDP2 and CFDP1_vs_CFDP2 respectively (the DE-miRNAs identified from adrenal cortex were listed in Table S7, and that from pituitary were listed in Table S8). Furthermore, a total of 90 miRNAs expressing differentially in pituitary comparisons were obtained, with 51 common DE-miRNAs shared in at least 2 pituitary comparisons (common DE-miRNAs). Of the 51 common DE-miRNAs, cfa-miR-105a, cfa-miR-219-3p and cfa-miR-802 were the common DE-miRNAs of the four pituitary comparisons. Meanwhile, of the 62 DE-miRNAs detected in adrenal cortex comparisons, there were 39 common DE-miRNAs, with the three common miRNAs of cfa-miR-196a, cfa-miR-216b, and cfa-miR-514 shared in all these four adrenal cortex comparisons. In consideration of the functional impact exerted by miRNA expression level, we chose the DE-miRNAs whose total reads within pituitary or adrenal cortex were greater than 1,000 for subsequent analysis, thus resulting in a total of 39 and 18 common DE-miRNAs in pituitary and adrenal cortex separately (Fig. 3 and Table 2).

To validate the miRNA expression levels obtained by high-sequencing, 4 random DE-miRNAs, including cfa-miR-30a, cfa-miR-124, cfa-miR-205 and cfa-miR-222 were selected to perform qPCR (Fig. 4). The results showed cfa-miR-30a, cfa-miR-205 both were significantly downregulated in CFDP1_vs_CFDP2 (|Log₂(fold change)| > 1). Cfa-miR-124 was significantly upregulated but not significantly for cfa-miR-222 in CFDP1_ vs_CFDP2 (|Log₂(fold change)| > 1). These observations were basically consistent with those obtained by high-sequencing, indicating the high-throughput sequencing data was reliable.

miRNA target prediction, and GO and KEGG analysis

The common DE-miRNAs, including 33 and 22 in pituitary and adrenal cortex were used for miRNA target prediction, with the methods of Targetscan used. Consequently, 3,959 and 4,010 target genes were obtained from pituitary and adrenal cortex respectively (Table S9). Then 3,959 and 4,010 target genes were performed with GO and KEGG analysis separately. Results showed that neither of the target genes from pituitary and adrenal cortex could be clustered into a certain GO term or KEGG pathway, due to the corrected p-value > 0.05. Instead, IPA was employed to identify the canonical pathways enriched by the predicted target genes. In CFDP1_vs_CFDP2 (i.e., 3,959 predicted target genes), the most enriched pathway was the axon guidance signaling pathway (p = 1.34E–11).

DE-(miRNA-mRNA) pairs identification

In our previous transcriptome study (Luo et al., 2015), a total of 511 and 171 DEGs were identified in pituitary and adrenal cortex, respectively (Table S1). The DEGs identified in our transcriptome study were separately compared with the target genes in pituitary and adrenal cortex (i.e., 3,959 and 4,010 respectively), thus resulting in 105 and 10 common genes that were both included in DEGs and target genes (DE-target genes), respectively. Furthermore, the expression patterns of the 105 and 10 DE-target genes were compared with that of their corresponding miRNAs, without prejudice to the negative correlations between miRNAs and their target genes in both breeds. Consequently, 18 DE-(miRNA-mRNA) pairs in pituitary were obtained (Table 3). In consideration of the miRNA expression level and its potential target genes’ relatedness to stress regulation, we selected 6 DE-(miRNA-mRNA) pairs to perform dual luciferase reporter assays to elucidate their dependability, including cfa-miR-30a-GRIA2, cfa-miR-30a-SORCS3, cfa-miR-30a-SLC1A2, cfa-miR-205-MMD, cfa-miR-205-GLRB and cfa-miR-205-CHN1. The dual luciferase reporter assay results showed just the cfa-miR-205-MMD conformed with the predicted (Fig. 5A), suggesting an important role played by cfa-miR-205 in regulating MMD under stress exposure in pituitary. Further, MMD and cfa-miR-205 expression levels were examined in the CFD brain, cerebellum, hippocampus and hypothalamus tissues by qPCR, normalized with GAPDH, U6 and U48 respectively. As expected, MMD was found to be upregulated, while downregulated for cfa-miR-205 in CFDP1_vs_CFDP2 (Fig. 5B).