MicroRNAs (miRNAs) are short non-coding RNAs that are widely expressed in the genome of humans, animals and plants [1]. It has also been demonstrated that miRNA is present in serum and plasma and are relatively stable and protected from degradation despite the presence of high ribonuclease activities [2-4]. There is evidence that these small non-coding RNA molecules play a role in many cell activities and that their expression level can be an indicator for specific patho-physiological conditions [5].
In the last decade, numerous studies have focused on miRNAs expression levels as biomarkers for pathophysiological disorders in humans and several domestic animal species. Circulating miRNAs expression levels may thus serve as biomarkers for physiological conditions, such as pregnancy, and pathological conditions during pregnancy [6-8]. Because of the high sensitivity and specificity and relative cost, the quantitative reverse transcription polymerase chain reaction (RT-qPCR) is currently the preferred technique for measuring the miRNAs expression levels [9]. However, the expression level analysis requires endogenous normalizers for normalization of the data. Thus, the selection of species-specific and tissue specific endogenous normalizers is critical for any study on miRNA expression [10-17]. The aim of normalization is to remove as much as possible the variation between samples, except for the difference that is a consequence of the condition itself. Variations between samples, might come from differences in RNA quality and integrity, or variations between the amounts of RNA which are loaded for each samples [17-19].
Normalizing transcripts have been well characterized for messenger RNA in virtually every tissue [20]. However, for miRNA there are not many normalizers or reference microRNAs that have been characterized [18, 21, 22]. Initially, small, noncoding RNAs such as 5s RNA and U6 sRNA were most commonly used as reference genes for miRNA quantification [22, 23]. Recently, there is growing evidence that the above-mentioned small RNAs are degraded in some serum samples and are highly variable among samples. These characteristics can cause significant bias when normalizing the miRNA expression and may lead to misleading results. Therefore, researchers are starting to look for more consistently expressed miRNA to be used as a normalizer in their studies [18]. It has been proposed that the best way to approach the analysis of miRNA expression data is to use a set of reference miRNAs that is tissue and species specific [24]. For instance, the miRNAs miR-23a and miR-191 are used for normalization in profiling studies of human cervical tissues [25], and let-7a and miR-16 have been selected as reference genes in human breast cancer tissues [26]. Also for quantification of circulating miRNAs in humans and mice, several studies have described the use of endogenous normalizers [27-30]. However, in livestock, few studies have addressed these issues. Studies using porcine tissues have reported that miR-93, miR-25, miR-106a, miR-17-5p, and miR-26a can be used as stable reference miRNAs studies [31, 32]. Recently, a combination of miR-93 and miR-127 has been proposed as an endogenous normalizers for circulating miRNAs in cattle [33].
Presently, no reference miRNAs for equine serum or chorioallantoic membrane tissue have been reported in the literature. The objective of the present study was to select a set of suitable reference miRNAs, which are consistently expressed in chorioallantoic membranes (CAM) throughout gestation, and separate set of reference miRNAs, which exhibit stable expression in serum of both pregnant and non-pregnant mares and can be used as normalization factors for miRNA qPCR analysis. Our hope was to identify stable miRNAs which can be used to normalize miRNAs expression level within both CAM and serum samples.
Material and Methods:
Serum and tissue collection and preparation
Serum from diestrus mares (n=3) and geldings (n=3) as well as from pregnant mares at 4, 6 and 10 months of gestation (n=3/stage) were collected. Samples were centrifuged at 1000 rpm for 10 min at 4°C. Supernatant was removed and stored at -20°C. Samples of chorioallantoic membrane were collected from pregnant mares at 4 mo (n=7), 6 mo (n=4) and 10 mo (n=7), as well as from the chorioallantoic membrane of spontaneously expelled placenta after normal parturition (n=3). The uterus of pregnant mares was recovered immediately after euthanasia and the chorioallantoic membrane was gently peeled off the endometrium. Full thickness samples were collected from this part of the placenta. CAM samples were stored overnight at 4°C in RNALater (Applied Biosystems, Carlsbad, CA) and then moved to −80°C until further processing.
RNA Extraction and cDNA Synthesis
Serum
Samples were thawed on ice and RNA was extracted starting from 400μL of serum using a miRNeasy Micro Kit (Qiagen) according to manufacturer’s instructions with small modification. TRIzol LS (Life Technologies, USA) was used as lysis reagent; 200 μL of chloroform was used instead of 140 µL; and the final elution of RNA was performed using 20 µL RNase-free water instead of 14 µL. Concentration of the RNA samples was measured using (1) the nanoDrop DP-1000 spectrophotometer (Agilent Technologies, Palo Alto, CA), (2) the Qubit® 2.0 Fluorometer (Life Technologies) and (3) the Bioanalyzer (Agilent Technologies, Santa Clara, CA, USA) 6000 Pico Kit.
Complementary DNA synthesis was done with the miScript II RT Kit (Qiagen, Valencia, CA, USA) using fixed volume of extracted RNA according to the manufacturer instructions. The 20 μL reverse transcription reaction mixture consisted of 10 μL of extracted RNA, 4 μL of miScript HiSpec buffer, 2 μL of miScript Nucleics mix, 2 μL of RNase-free water, and 2 μL of miScript Reverse Transcriptase mix. The reaction mixture was incubated for 60 min at 37 °C and then 5 min at 95 °C, using a thermal cycler. Three microliter of each sample was added to a tube to make the pooled serum sample.
Chorioallantoic membrane
Samples were thawed on ice and RNA was extracted using 75 mg of CAM sample in 1mL of Trizolâ„¢ Reagent (Invitrogen, Carlsbad, CA, USA) according to the manufacturer recommendation. RNA quantity was measured by NanoDrop DP-1000 spectrophotometer (Agilent Technologies, Palo Alto, CA).
Complementary DNA synthesis was carried out using miScript II RT Kit (Qiagen, Valencia, CA, USA). The 20 μL reverse transcription reaction mixture consisted of 2 μL of extracted RNA (containing 400ng of extracted RNA), 4 μL of miScript HiSpec buffer, 2 μL of miScript Nucleics mix, 10 μL of RNase-free water, and 2 μL of miScript Reverse Transcriptase mix. The reaction mix was incubated for 60 min at 37 °C and then 5 min at 95 °C, using a thermal cycler. Three microliter of each sample was added to a tube to make the pooled CAM sample.
RT- qPCR
Mature miRNAs was detected using RT-qPCR. The RT-qPCR was performed using miScript SYBR Green PCR kits (Qiagen), containing a QuantiTect SYBR Green PCR Master Mix and the miScript Universal Primer along with the miRNA-specific primer according to the manufactures guidelines using the following cycling conditions: 95°C for 15 min; 40 cycles of 94°C for 15 sec, 55°C for 30 sec; 70°C for 30 sec. After the PCR cycles, a DNA melting curve was generated in order to discriminate between specific and non-specific amplification products.
Initially, serial dilutions of pooled cDNA from serum were used for each primer set to adjust the assay. The different volumes of cDNA (0.06, 0.12, 0.25 and 0.5 μL) were used to determine the optimal volume of cDNA for qPCR assay. Geometric mean of CT values for each dilution was calculated. Optimal volume of cDNA (0.25μL) was selected for serum samples according to the lowest geometric mean of CT values.
RT-qPCR was performed in triplicate for serum and duplicate for chorion samples. PCR efficiency was calculated using LinRegPCR (version 2012.0) (http://www.hartfaalcentrum.nl) to ensure that all primers resulted in PCR efficiencies of 1.8 – 2.1.
Selection of Candidate Reference miRNAs
Potential candidate miRNAs were selected from a previously generated microRNA sequencing dataset, from chorioallantoic membrane and serum from pregnant (4 and 10 mo), postpartum and diestrus mare serum (unpublished data). Read counts were normalized to average miRNA expression for each miRNA for serum and CAM separately. The selected potential reference miRNAs were selected using NormFinder software (V 0.953) from list of 1185 miRNAs for CAM and list of 653 miRNAs for serum. NormFinder was chosen because of the ability of its algorithms to differentiate intragroup variation from intergroup variation. This makes NormFinder a suitable tool for identifying candidate miRNAs in bigger sample sizes compared to geNorm [34].
The top 20 candidates were selected for serum and CAM separately, and primers were designed using miRprimer software (version 2.0) [35]. 5S rRNA and U6 snRNA, historically used miRNA normalizers, were also used. Primers that were designed for serum were first tested on pooled samples made from serum. Designed primers for CAM were also first tested on a pooled sample made from all CAM samples. Primers, which were found to be reliable and did not produce primer-dimers, were used for further analyses.
Data analysis for miRNA stability
The assessment of the putative reference miRNA for normalization was initially evaluated using geNorm through the R-based SLqPCR package (Version 1.0.0) [21, 36] and NormFinder (version0.953 ) [37]. The correlation between the “stability valueâ€, generated by NormFinder and the “M value “, generated by geNorm was calculated by Pearson correlation coefficient. Spearman`s rank-order correlation coefficient was also performed on the miRNAs ranking with geNorm and NormFinder using the function “cor.test†with the method set as “spearman†in the “Hmisc†package. In case of a discrepancy between the results of geNorm and NormFinder, BestKeeper (version 1.0), was used to find the potential reference miRNA [38].
The geometric means of the Cq values run in triplicate (in serum) and duplicates (in CAM) were further used for the comparison of the stability of each candidate reference miRNA. For the NormFinder, the Cqs were transformed to relative quantities (RQ) using the following formula, in which RQ is the relative quantity of the miRNA of interest, E is the mean PCR efficiency, “minCq†is the lowest estimated Cq of the miRNA and “SampleCq†is the Cq of each samples [39, 40]. RQ= E -(minCq-sampleCq). The resulting RQs were used as input data in the Normfinder, to calculate stability values in the samples under analysis. The lowest stability value represents the lowest variation and the highest stability in expression
In R-based SLqPCR package (version 1.0.0), Cq values were transformed to relative quantities (RQ) via function relQuantPCR, in which (the lowest relative quantity for each gene was set to 1.0) then used as a input data in function selectHKgenes to generates a ranking of genes according to their M values, resulting in the identification of the genes with the most stable expression in the samples under analysis; the lower the M value, the higher the gene’s expression stability [36].
The BestKeeper algorithm using raw Cq values for each sample for pairwise correlation, and calculates the coefficient of variance of the Cq values and establishes the correlation coefficients to estimate a reference miRNAs index [38, 41].The final ranking of the BestKeeper software is usually performed by assessing the correlation coefficients of each individual gene with the geometric mean of Cq values of all candidate reference genes (the BestKeeper Index) [34]. Hence, the correlation coefficient is used to assess the most stable miRNAs [38].
Each miRNA was ranked on the basis of the results obtained from the analysis with different software, which a rank of 1 indicating the most stable gene, and the arithmetic average of each gene rank was calculated and reported as the average ranking stability [42].
Statistical analysis
In order to look at the correlation in the miRNA concentration measured with different method, Pearson correlation coefficients were calculated on continuous variables using the function “cor†of the “Hmisc†package [43].
Results
MiRNA concentration and quality
The concentration of extracted RNA from CAM samples was within the manufacturer’s stated detection limit for the Nanodrop; ranging between 471 to 9133 ng/μL with median of 1605ng/μL,
The concentrations of extracted RNA from serum samples were below the detection limit of the Nanodrop. There was a high correlation in serum miRNA concentration (r= 0.928, P <0.01) between Qubit and Bioanalyzer, however measured concentration was still too low to monitor yield and ensure consistent miRNA input across all samples (concentration ranged between 0.14 ng/μL to 1.08ng/μL measured with Qubit and between 21 pg/μL to 1.189 pg/μL measured with Bioanalyzer). For this reason, equal volume of extracted miRNA for all samples were used for qPCR.
Selection of Candidate Reference miRNAs and Primers
Initially, all the primers were evaluated by testing them on pooled samples of either CAM or serum. Primers were deemed reliable if they had a Cq < 37 [44], and did not produce primer dimers (Cq=NA in water sample). Primers meeting these criteria were used for further experimentation; these included nine serum primers and 10 CAM primers (Table 1).
Expression levels of candidate reference miRNAs
As shown in Figure 1, the quantification cycle (Cq) values of candidate miRNAs ranged from 16 to 30 for CAM samples and from 22 to 35 for serum samples. The range of Cq values indicates the stability of reference miRNAs expression across the samples.
In serum, the range of Cq values for eca-miR-125-2-5p, eca-miR-10b-5p and eca-miR-10a-5p was less than 6 cycle. Eca-miR-128-1-3p had the least stability across serum samples with a range of 10.5 cycles. Percentage of coefficient of variation (CV%) between the samples for each candidate was calculated. Eca-miR-125-2-5p, eca-miR-10b-5p and eca-miR-10a-5p had the lowest CV% and eca-let-7g had the highest CV% (3%, 4.4%, 5.8% and 9.4% for eca-miR-125-2-5p, eca-miR-10b-5p, eca-miR-10a-5p and eca-let-7g, respectively).
In CAM, the range of Cq values for eca-miR-7-3-5p, 5S-rRNA and eca-miR-500-1-3p was less than 6 cycle. Eca-miR-19b-2-3p had the widest range of Cqs between all the potential candidates for CAM reference miRNA. Eca-miR-7-3-5p, 5S-rRNA and eca-miR-8908a-1-5p had the lowest coefficient of variation percentage (CV%) between all the CAM samples (5.5%, 5.8% and 6.4% for eca-miR-7-3-5p, 5S-rRNA and eca-miR-8908a-1-5p, respectively). Eca-miR-19 had the highest CV%, (13.7%) among all the candidates.
Selection of suitable endogenous normalizer
The expression stability of all candidate miRNAs were initially evaluated by NormFinder and geNorm.
NormFinder
In serum, the three most stable miRNAs were eca-miR-21-5p, eca-miR-10a-5p and eca-let-7a. In CAM, the three most stable miRNAs were eca-106a-5p, eca-miR-500-1-3p and eca-miR-8908a-1-5p (Figure 2).
geNorm
In serum, the top three-ranked miRNA were eca-miR-21-5p, eca-let-7a and eca-miR-10a-5p and in CAM the top three-ranked miRNA were eca-miR-500-1-3p, eca-miR-130a-3p and eca-miR-8908a-1-5p (Figure 3).
Correlation between NormFinder and geNorm
MicroRNAs were ranked according to the results in NormFinder and geNorm, in each ranking a rank of 1 indicates the most stable gene (Table 2).
Spearman correlation demonstrated high correlation between miRNA ranking of geNorma and NormFinder in serum, (rho=0.97). In contrast, there was no correlation between the ranking of miRNAs between geNorm and NormFinder for CAM (rho = 0.56).
Similar results were obtained by using Pearson correlation coefficient between geNorm average expression stability value “M†and NormFinder “Stability valueâ€; r=0.967 with p-value< 0.01 for serum and r=0.43 with the p-value=0.1 for CAM candidates.
Preparation of list of CAM`s potential candidate miRNA for BestKeeper
Since we were not able to demonstrate a correlation between selected miRNA with geNorm and NormFinder in CAM, BestKeeper was used to come up with a third ranking. BestKeeper has the ability of computing only 10 candidates, therefore the three least stable miRNAs identified by geNorm (eca-miR-7-3-5p, eca-miR-377-3p, U6-snRNA) and NormFinder (eca-miR-19b-2-3p, eca-miR-377-3p, 5S-rRNA) were excluded and the remaining 8 candidates were subjected to the BestKeeper analysis.
After excluding the top three least stable miRNAs identified by geNorm and NormFinder, the expression stability for the remaining 8 candidate miRNAs was also re-evaluated by geNorm and NormFinder (Figure 4).
BestKeeper
The three most stable miRNA in CAM samples were eca-miR-8908a-1-5p, eca-miR-106a-5p and eca-miR-369-5p (Table 3).
Correlation between three methods
Each CAM putative reference miRNA was ranked again on the basis of the results obtained from the new analysis with geNorm, NormFinder, and BestKeeper, with a rank of 1 indicating the most stable gene [42]. The average of each gene rank was also calculated and reported as the average ranking stability (Table 4). Spearman’s rank-order correlation coefficient was performed on the miRNAs rankings obtained in the three programs. Average ranking stability was also included in the computation. There was good correlation between the ranking of miRNAs candidate in all three methods and the average ranking, (r= 0.96, 0.9 and 0.9; between average ranking and geNorm, NormFinder and BestKeeper, respectively; Figure 5).
Discussion
MicroRNAs are important regulators of cell and tissue behavior. In the reproductive system, different miRNAs have been implicated in the development of ovarian follicles and corpus luteum, uterine cyclicity, establishment of pregnancy, embryonic development and placental development [6, 45-49]. Evaluating the expression of miRNAs during various conditions is the most common way to study miRNA dynamics. While RT-qPCR is widely applied to evaluate changes in miRNA expression levels, the accuracy of this method is easily influenced by cDNA quality, RNA integrity and PCR efficiency. Use of reference miRNA as normalizers for relative quantification of miRNA expression is therefor indispensable to avoid bias in expression results and to make comparison of miRNA expression data more reliable and meaningful. Due to the fact that a universal housekeeping gene/miRNA does not exist, one or more endogenous normalizers need to be validated independently for each set of tissues [21, 50]. The aim of this study was to find accurate, reliable reference miRNAs for chorioallantoic membrane and serum at different stage of pregnancy in the mare.
In the present study, RNA was extracted from equine serum samples and chorioallantoic membrane. The RNA concentration of CAM samples was successfully measured by Nano Drop. However, in serum samples the extracted, RNA concentrations were under the detection limits. It has been reported that accurate measurement of RNA concentration by spectrophotometric methods for serum samples is impossible [51]. To measure the concentration of the extracted RNA from serum, Qubit® and Bioanalyzer were used. There was a high correlation between the measured concentrations between two methods. However, measured concentration was too low to ensure consistent miRNA input across all samples. This finding was previously described by various authors [10, 52, 53] yet, there are also some reports of measuring serum miRNA concentration with Qubit and NanoDrop [54, 55]. To overcome the limitation of unknown starting concentration, we used fixed volume of the extraction for RT reaction as described previously [10, 48].
Initially all the primers for putative candidates were tested on the pool samples. Primers that had a Cq greater than 37 were excluded from the list. Also those primers, which formed primer dimers, were excluded from the list. Consequently, thirteen primers for CAM and nine primers for serum samples were selected for subsequent assessment. Series of algorithms have been developed to for comparison of the stability of candidate reference genes [34]. Because there are no universally accepted algorithm and use of a single algorithm may introduce bias, we used the two most widely used methods namely NormFinder and geNorm, to compare the stability of normalizers in the present study.
In serum miRNAs candidates, there was a high correlation between geNorm and NormFinder finding, and both programs indicated that eca-miR-21-5p was the most stable candidate. Rankings for all candidates followed the same order in both program except slight variation in the position of eca-let-7a-5p and eca-miR-10a-5p. Eca-miR-125-2-5p showed the lowest stability with both programs, and surprisingly, had the least CV% and standard deviation between all the candidates when looking at the raw data. In the geNorm software, the gene expression stability threshold value (M) is set at <1.5, and candidate reference genes above this threshold should be considered as non-suitable as reference gene [56]. In the present study, the M-value of eca-miR-125-2-5p is proximate to the threshold (M value = 1.4). According to both programs, eca-miR-21-5p, eca-let-7a-5p and eca-miR-10a-5p were the most stable candidates and most appropriate to serve as normalizer.
A recent sequencing study of serum from three gelding horses showed that miR-21-5p and eca-let-7a-5p are among the 10 miRNA species exhibiting the highest plasma expression. This makes these miRNAs good candidates for reference miRNAs in serum [57]. MiR-21-5p was found to be potentially involved in mice trophectoderm and trophoblast development and function [49, 58]. There are also some reports regarding the function of miR-21-5p in promoting cell proliferation, migration, and invasion in head and neck carcinoma, osteosarcoma and ovarian carcinomas [59]. Likewise miR-21 but to a smaller extend let-7a and miR-10a are playing roles in many reproductive tract functions; let-7a-5p expressed in bovine endometrium across the estrous cycles [60]. Let-7a-5p was also shown to be elevated in the endometrium of 20-day pregnant mouse [61]. Let-7a-5p also participates in endometrial decidualization during the window of implantation and affects the implantation competency of the activated blastocysts [61-63]. Moreover, miR-10a-5p is associated with successful embryo implantation in women [64].
In the CAM samples there was insufficient correlation between NormFinder and geNorm and therefore, BestKeeper was used as the third algorithm to detect appropriate reference miRNAs. The BestKeeper algorithm can only compare 10 candidate genes. Hence, the three miRNA identified as the least stable with geNorm and NormFinder were excluded from the candidate list. After excluding the least stables candidates, geNorm, NormFinder and BestKeeper algorithms were used again to build new candidate lists. The BestKeeper software calculates the correlation coefficients of each individual gene with the geometric mean of all genes. This different basis for determining gene stability results in a substantially different gene ranking compared with the other two algorithms [34]. Despite the differences among the three algorithms, the rankings of the most stable candidates were similar. All the proposed candidates had an M value less than 1 and standard deviation less than 1.5, (measured in geNorm and Bestkeeper, respectively) which are less than proposed gene expression stability threshold value (M<1 and SD <1.5) [34] and was only exceeded by eca-miR-8908e-5p. When combining the output of the three algorithms by calculating the average ranking, it was determined that eca-miR-8908a-1-5p was the most stable miRNAs followed by eca-miR-369-5p and eca-miR-106a-5p.
Similarly, in serum it has been reported in other species that miR-369-5p and miR-106a-5p expression is altered in specific pathologic conditions. Expression of miR-369-5p is highly elevated in urine of human patients with intrahepatic cholestasis of pregnancy and in rats with pancreatic ductal adenocarcinoma [65, 66]. Also, miR-369-5p is involved in mesenchymal stem cell phenotype and epithelial-mesenchymal transition in endometrial carcinosarcomas [67, 68]. Mir-106a-5p is among the most abundant miRNAs detected in the plasma of healthy women during the menstrual cycle [69] and miR-106a-5p is the most abundant miRNAs in placental exosomes [70].
Previously, numerous RNA species, including rRNA and snRNA have been used as normalizers. However there is increasing proof that these small RNAs are highly variable among samples and are not good candidate normalizer for miRNA expression [18, 71-78]. The result of this study is in agreement with this finding. In serum samples, 5S rRNA and U6 snRNA were not expressed within acceptable cycle threshold. In CAM samples even though 5S rRNA expression had one of the lowest coefficient of variation and standard deviation between samples, but it was not identified as a suitable candidate with both algorithms.
The most stable miRNAs identified in the present study were detected in samples collected during the physiological status and cannot be applied in other studies using similar tissues collected in pathological condition without prior revalidation. Consequently, the provided list is a set of putative candidates for studying pathological conditions in CAM and serum samples during the equine pregnancy in equine species.
In conclusion, we were able to identify a set of reference miRNAs for chorioallantoic membranes and serum during different stage of equine pregnancy that can be used as normalizing factors for miRNA RT-qPCR experimental analysis