Abstract
MicroRNAs (miRNAs) are highly conserved 22-mer RNA molecules, encoded by plants and animals that regulate the expression of genes binding to the -UTR of specific target mRNAs. The amount of miRNAs in a total RNA sample depends on the recovery efficiency that may be significantly affected by the different purification methods employed. Traditional approaches may be inefficient at recovering small RNAs, and common spectrophotometric determination is not adequate to quantify selectively these low molecular weight (LMW) species from total RNA samples. Here, we describe the use of qualitative and quantitative lab-on-a-chip tools for the analysis of these LMW RNA species. Our data emphasize the close correlation of LMW RNAs with the expression levels of some miRNAs. We therefore applied our result to the comparison of some miRNA expression profiles in different tissues. Finally, the methods we used in this paper allowed us to analyze the efficiency of extraction protocols, to study the small (but significant) differences among various preparations and to allow a proper comparison of some miRNA expression profiles in various tissues.
1. Introduction
In the last few years, a new class of highly conserved -mer non-coding RNAs, microRNAs (miRNAs), has emerged as an important player in post-transcriptional gene expression control in different physiological and pathological conditions: metabolism, proliferation, cell death, differentiation and development, viral infection, and cancer [1–3]. They specifically bind the -untranslated regions (-UTR) of target mRNAs, promoting either mRNA degradation or translation arrest [4–6]. Each miRNA may control the activity of many genes, and almost 30% of the genome could be regulated in such a way, which renders these small molecules as important as the transcription factors [7, 8].
Owing to their extreme importance as regulators, the isolation and the precise quantification of these tiny molecules are therefore fundamental. As purity and integrity are essential requirements not only for total RNA but also for these small species, RNA extraction protocols must also account for the maximizing of their recovery. Sometimes, the handling and recovery of small RNAs are not straightforward. Nowadays, it is well established that the traditional glass-fiber total RNA extraction protocol may be inefficient at recovering small RNAs. In addition, the common spectrophotometric determination of total RNA is not adequate to quantify low molecular weight (LMW) species selectively.
We and other authors found previously that the recovery of LMW RNA species is significantly affected by the specific purification process [3, 9]. Our results showed that different extraction strategies lead to significantly different recovery of LMW species including miRNAs. We also demonstrated that using the same amount of total RNA (from different tissues), different amounts of miRNAs may be obtained. Moreover, the amount of LMW RNA species does not perfectly parallel that of miRNAs: even with the same extraction protocol, the concentration of miRNAs may differ significantly among various tissues. Nevertheless, this may also be because of a different global expression of miRNA by different tissues. All of these considerations are extremely important when one compares the relative expression of many miRNAs in different tissues. Commonly, some small non-coding RNAs (i.e., U6, snoZ30, and others) are used as normalization controls to compare miRNAs expression in different tissues. The choice of a proper control should be conditioned to its real “endogenous behavior”, that is, it should have a constant and equal expression in all tissues. If the amount of that control small non-coding RNA is different, a relative quantification is not feasible. In this latter case, only the comparison between various miRNAs and the control gene in the same specimen (tissue) is methodologically acceptable. Therefore, a method allowing a reliable estimation of the goodness of an endogenous control is highly desirable.
In order to address all of these issues we used a lab-on-a-chip technology to characterize the LMW RNA fractions obtained with different RNA extraction protocols and evaluated miRNAs recovery with specific quantitative real-time PCR (qRT-PCR).
In this paper we report and discuss the methods that allowed us to analyze the efficiency of extraction protocols, to study the small (but significant) differences between various preparations and to obtain a correct comparison of some miRNA expression profiles in various tissues.
2. Materials and Methods
2.1. Total RNA Extraction and Small RNAs Enrichment Protocols
Total RNA was extracted using three different methods: an acid phenol/guanidine isothiocyanate solution (TRIzol Reagent, Invitrogen), a glass-fiber filtration protocol (MirVana miRNA Isolation Kit, Ambion) that provides also a procedure to isolate and enrich low molecular weight (LMW) RNAs from higher molecular weight (HMW) RNAs, and another common glass-fiber purification protocol (RNEasy Mini Kit, Qiagen). All extractions were performed according to manufacturer’s instructions.
2.2. Cell Lines
Three different cell lines (HeLa, COS-1 and a lymphoblastoid cell line (LCL) obtained following an already reported method [10]) were cultured using standard procedures, trypsinized (if necessary), and pelleted by centrifugation. Approximately cells for each extraction were resuspended in the appropriate lysis solution contained in the RNA extraction kit and treated according to manufacturer’s instructions.
2.3. Electrophoresis and RNA Elution from Gel
The integrity of RNA samples was checked by gel electrophoresis (agarose 1%) stained with ethidium bromide. Gel images were acquired and analyzed with the Quantity One (software Ver.2.0, Biorad). Two micrograms of each RNA samples were resolved into a 15% polyacrylamide TBE-Urea gel and the image acquired. A second polyacrylamide gel, performed with the same experimental conditions, was used for RNA bands extraction. The visible bands were separately excised using a modified “crush-and-soak” method [11]. Briefly, each band was gently disrupted with a pin and incubated in a soaking solution (0.5 M ammonium acetate, 2 mM EDTA, 0.1% SDS) for 3 hours at C. The mixture was centrifuged at 14000 g for 5 minutes, and the supernatant was collected in a new tube. The solution was centrifuged at 12000 g for 10 minutes through a YM-50 spin column, washed with water , and centrifuged at 12000 g for 5 minutes. Purified fractions were collected in new tubes by inverting the columns and spinning at 1000 g for 3 minutes. These RNA samples were then separately run on the Agilent Bioanalyzer.
2.4. Agilent 2100 Bioanalyzer Chip Preparation
Total RNA samples were analyzed with the Total RNA 6000 Nano Kit (Vers. II), specifically optimized for total RNA analysis with the Agilent 2100 Bioanalyzer. RNA samples obtained after extraction from polyacrylamide gel were run with the same kit. For miRNA quantification we used the dedicated Small RNA kit. The instrument uses fluorescence detection, monitoring the emission between 670 nm and 700 nm. The run was performed according to manufacturer’s instructions. Electropherograms were analyzed using the Agilent 2100 Expert B.02.06 software that includes data collection, presentation, and interpretation functions.
2.5. Real-Time PCR
For the small RNAs recovery efficiency evaluation of different extraction procedures, hsa-miR-21 and small nucleolar Z30 (snoZ30) were assayed by Real-Time PCR (TaqMan, Applied Biosystems), according to manufacturer’s instructions [12]. For the analysis of miRNA expression in different tissues, three miRNAs (hsa-miR-26a, hsa-miR-26b, and hsa-miR-134) and two controls (U6 and snoZ30) were arbitrarily chosen as model miRNAs. A panel of five tissues (brain, skeletal muscle, heart, liver, and uterus) was chosen for evaluation, and the corresponding total RNA (Clontech, BD Biosciences) was analyzed with Agilent 2100 Bioanalyzer. Briefly, from five to ten nanograms of each RNA samples were retrotranscribed to cDNA using the High-Capacity cDNA Archive Kit (Applied Biosystems) with specific primers. The principle of TaqMan MicroRNA Assays is a specific stem-loop reverse transcription (RT) primer. The short length of mature miRNAs ( nt) prohibits conventional design of a random-primed RT step followed by a specific real-time assay. In the former case, the resulting RT amplicon is a suitable template for standard real-time PCR with TaqMan assays. Reactions were performed incubating samples for 30 minutes at C, 30 minutes at C, 5 minutes at C, and finally cooled on ice. PCR products were assayed with specific probes using the TaqMan Universal PCR Master Mix, No AmpErase UNG (Applied Biosystems) following the manufacturer’s protocol. PCR reactions were performed incubating samples for 10 minutes at C, then for 15 seconds at C and 60 seconds at C for 45 cycles by means of ABI PRISM 7900HT Sequence Detection System. Data were analyzed using the SDS software (Ver. 2.1).
2.6. Deconvolution with PeakFit Software
The early electropherogram region (from 23 to 29 seconds), obtained from LCL RNA sample extracted with TRIzol reagent, was exported in ASCII format with the Agilent 2100 Expert B.02.06 and imported in PeakFit V.4.12. The curve was automatically fitted using the least-squares method. Data were further adjusted using the Savitzky-Golay smoothing algorithm. The overall model was fitted with chromatographic Gaussian curves of variable amplitudes. The fitting procedure was iterated until a constant value was obtained.
3. Results
3.1. The Recovery of Low Molecular Weight (LMW) RNAs Is Affected by Different Extraction Protocols
RNA samples extracted from HeLa, COS-1, and LCL were run on agarose gel to visualize the differences between various extraction methods. COS-1 RNA samples, extracted with TRIzol reagent and MirVana kit, clearly showed the High Molecular Weight (HMW) 28S and 18S rRNA bands, while LMW RNAs are visualized as faint, smeary bands (Figure 1, lanes 1 and 2, resp.). COS-1 extracted with RNeasy kit displayed only the HMW RNA bands (28S and 18S) (Figure 1, lane 3), while the enriched LMW fraction obtained with MirVana kit is clearly displayed in the lowest part of the gel (Figure 1, lane 4). Similar results (not shown) were also obtained with the other cell lines.
The same samples were also checked with Agilent 2100 Bioanalyzer which is one of the most versatile microfluidics-based platforms for the analysis of DNA, RNA, proteins and cells. In all electropherograms the 28S and 18S RNAs are represented on the right side (Figure 2), and the smaller species (LMW RNAs) are present at a very low concentration and distinguishable on the left side of the profile (see the magnification of the LMW RNAs region in Figure 2). In the electropherogram, all samples present a similar HMW profile region, but the major differences are recognizable in the LMW region. While total RNA recovery is quite good, and similar for the three protocols , the LMW fractions are substantially different (Table 1). In particular, TRIzol reagent allowed the highest LMW RNA recovery (22%–34% of total RNA), while RNEasy Mini Kit the lowest (2.5%–3%). MirVana miRNA Isolation Kit gives good yields for LMW RNA species (16%–19%) even before the enrichment step. While LMW RNA species extracted with TRIzol and MirVana have comparable profiles (Figure 2); RNEasy kit retains only one RNA peak in comparable concentrations to the others (peak in the magnification of Figure 2). Therefore, the lab-on-a-chip analysis is a useful tool to quantify precisely the amount of LMW RNAs of samples extracted with different protocols.
3.2. Low Molecular Weight (LMW) RNAs Characterization
Total RNAs from LCL extracted with TRIzol reagent, and RNEasy kit and enriched with MirVana kit were loaded on a polyacrylamide gel (Figure 3(a)) to analyze the differences between various extraction protocols. All the RNA bands were cut out from the gel, and the RNA was extracted as described. Each species coming from excised bands was run on Bioanalyzer, and the corresponding electropherograms (Figure 3(b)) were superimposed on the large and unresolved LMW RNA profile (black line). This allowed us to localize unambiguously and identify precisely each RNA species previously identified on polyacrylamide gel.
(a)
(b)
(c)
(a) Gel electrophoresis (polyacrylamide 15% stained with ethidium bromide) of LMW RNA samples (LCL) extracted with TRIzol reagent (lane 1), RNEasy kit (lane 2), and small RNA fraction enriched with MirVana kit (lane 3). LMW profile obtained with MirVana kit extraction is similar to that obtained with TRIzol reagent which is not shown for clarity. (b) Agilent 2100 Bioanalyzer electropherogram profile of LMW RNAs (LCL) extracted with TRIzol (black) superimposed on 5.8S (red), 5S (green), and tRNA (blue) bands eluted from polyacrylamide gel. (c) Lymphoblastoid (LCL) LMW RNA profile obtained after plotting the exported raw data from Agilent electropherogram ( ) together with the PeakFit fitted curve (solid line) and the component peak functions. Seven peaks below the LMW RNA profile were fitted by the software . |
The electropherogram profile in the LMW RNAs region is made of several peaks (Figure 3(b)), the greatest part belonging to the 5.8S rRNA (red line), the 5S rRNA (green line), and the tRNAs fractions (blue line). In order to obtain the principal number of peaks under the curve, PeakFit software was used for the deconvolution of the LMW RNA region. Seven peaks underlying the LMW RNA profile were obtained with a good fit . The upper part of Figure 3(c) reports raw data () together with the fitted curve (solid line) and the calculated baseline, while the lower part of Figure 3(c) shows the component peak functions. Calculated peak-peaking (expressed in bp) is also reported (Figure 3(c)). Three of the calculated peaks centered at 145 bp (26.89 seconds), 131 bp (26.33 seconds), and 104 bp (25.24 seconds) are located in the correspondence of 5.8S, 5S, and tRNA peaks, respectively. Another less intense peak is located at 114 bp (25.83 seconds) and may be described as a further tRNA peak. Other peaks found at 171 bp (27.78 seconds), 89 bp (24.84 seconds) and 65 bp (24.31 seconds) might belong to other LMW RNA species, such as smaller tRNAs, pre-miRNAs, and small nuclear or nucleolar RNAs. These small peaks are difficult to extract from polyacrylamide gel owing to their low abundance, and therefore it was not possible to run them on the chip.
Therefore, determination of small RNA species by Bioanalyzer (Table 1) and by polyacrylamide gel electrophoresis (Figure 3(a)) showed that the LMW RNA samples significantly differ from each other. This could reasonably pertain also to the miRNAs fraction.
3.3. MicroRNAs Amount Correlates with the LMW RNA Fraction
In order to quantify miRNAs extracted with different protocols and to study the correlation between LMW RNAs and miRNAs, we carried out a TaqMan quantitative assay for two representative small RNAs: a control non-coding RNA (snoZ30) and a miRNA (hsa-miR-21). A target-specific stem-loop adapter technology was employed to obtain the corresponding cDNA [12].
We started the RT reaction with 10 ng of total RNA from each sample and quantified the absolute expression level of each miRNAs through the analysis of cycle threshold (Ct) values. Ct is the PCR cycle at which the sample reaches the level of detection above the background. LCL RNA samples extracted with different protocols showed different Ct values for snoZ30 (Figure 4(a)) and for miR-21. A similar behavior was also obtained by using HeLa and COS-1 cell lines with both probes (data not shown).
(a)
(b)
Then, we repeated the RT reaction using 10 ng of LMW RNA calculated from the Bioanalyzer electropherogram. As expected, using the same LMW RNA amount similar Ct values both for miR-21 and snoZ30 between the various samples were obtained (Figure 4(b)). This demonstrates that the amount of miR-21 and snoZ30 correlates only to the LMW RNA fraction and not to total RNA amount. In fact, Figure 4(a) clearly shows that using the same amount of total RNA different Ct values may be obtained.
Ct differences between total RNA and LMW RNA Ct values for each extraction protocol and for each cell line are very similar both for hsa-miR-21 and snoZ30 (Table 2). This demonstrates that differences parallel the amount of input RNA that is ultimately linked to the amount extracted by various protocols. Table 2 reports also the calculated LMW RNA concentration (expressed in %) calculated from values applying the formula . MiRNA concentrations calculated from real-time assays are, as expected, in perfect agreement with LMW RNAs concentrations evaluated with the lab-on-a-chip technology.
3.4. Different Tissues Express Different Amounts of Small- and MicroRNAs
Here, we showed that LMW RNAs and miRNAs amount are closely correlated. Then, we asked if different tissues might also express different amounts of LMW and miRNAs and if it might be possible to quantify them individually with lab-on-a-chip technology. For this reason we selected five commercial RNA samples from different tissues extracted with the same protocol (according to manufacturer information). Total RNAs from brain, skeletal muscle, heart, liver, and uterus were run on RNA 6000 Nano kit to quantify the total and the LMW RNA fractions. Total RNA concentrations were quite homogeneous and resulted: brain () heart () uterus () liver () skeletal muscle (). These values are in good agreement with the nominal manufacturer’s concentration of . For all samples the concentration of LMW RNA fraction was quite homogeneous and varied in the order: liver () brain () heart () uterus () skeletal muscle (). Electropherograms of LMW RNAs of different tissues are displayed in Figure 5(a). Calculating the percentage of LMW RNAs with respect to the whole total RNA amount we found the following: liver (4.4%) skeletal muscle (2.6%) brain (2.5%) heart (2.3%) uterus (2.3%). Only some minor differences may be observed among various tissues. Given that the extraction protocol employed by the manufacturer was the same for all samples, we could hypothesize that the only observed difference for liver might depend on the different global expression of small RNAs in this tissue.
(a)
(b)
To evaluate if the difference in LMW RNA amount may also pertain to the miRNA fraction, we ran the total RNA samples on Agilent Small RNA kit, specifically designed for the evaluation of miRNAs. Figure 5(b) shows a magnification of the electropherogram profile obtained by running total RNA samples from different tissues. The displayed region (from 35 to 45 seconds) is specific to the miRNA region (as indicated by the manufacturer). MiRNA concentrations were liver () brain () uterus () heart () skeletal muscle (). The miRNA concentration, expressed in percentage, respect to LMW RNAs resulted brain (1.7%) liver (1.6%) uterus (1.6%) heart (0.9%) skeletal muscle (0.3%). From this lab-on-a-chip quantification it was possible to conclude that the miRNAs amount does not parallel that of LMW RNAs in the same tissue. Moreover, the variation is greater than that displayed by LMW RNAs among different tissues. Again, assuming that the extraction efficiency is the same for all samples, the observed differences are only due to the different global expression of miRNAs in these tissues.
3.5. Evaluation of Endogenous Controls’ Reliability and Their Use for Expression Profile Comparison
The Small RNA kit, specifically designed for the identification and quantification of miRNA species, allowed us to know the exact amount of miRNAs present in the samples used for retrotranscription (RT) reactions. This information is useful, since the Ct values obtained from real-time assays are directly proportional to the amount of miRNAs. In fact, the Ct differences observed for the same miRNA in different tissues may be due essentially (1) to a real difference between different samples and/or (2) to a different sample dilution that generates a shift in Ct absolute values. To eliminate the dilution problem, that occurs when absolute quantification is performed, a relative quantification with respect to an endogenous control is commonly followed. In fact, the endogenous control must have a constant expression in all samples, and hence it may be used to normalize the expression of the other miRNAs. For a methodologically correct comparison, however, one must be confident that the control (i.e., nuclear and/or nucleolar small RNAs) really does have a constant expression in considered samples. The lab-on-a-chip technology that we used is able to give an estimate of this variation and let the researcher choose the right endogenous control (the one that does not significantly vary) from an adequate selection. Therefore, this validation ensures that the miRNA expression comparison among considered samples is methodologically correct.
The expression of three miRNAs (hsa-miR-26a, hsa-miR-26b, and hsa-miR-134) and two small RNAs (U6 and snoZ30) from the same five tissues analyzed before, was assayed with real-time PCR. According to the manufacturer’s suggestion, we started the RT reactions with 5 ng of total RNA. Cycle threshold values for all tissues are reported in Table 3. As expected, different Ct values for endogenous and other miRNAs were obtained. These values reflect the absolute concentrations of these miRNAs in various samples. Since most studies aim to discover differences in expression levels of miRNAs and not absolute levels of expression, the use of an endogenous control is needed. In order to assess if the differences in Ct values of the controls we used (U6 and snoZ30) are because of a different starting concentration or a real differential expression, we corrected the obtained values by taking into account the concentration of miRNAs previously obtained with the lab-on-a-chip technology. Hence, we considered that for double the concentration a correction of one Ct value should be applied. This preliminary correction eliminated the intrinsic variability owing to different sample concentrations and allowed us to estimate the reliability of the selected endogenous controls. We observed slight variation for Ct values of U6 (Average ) while snoZ30 displayed a higher variability (Average ) even after the applied correction. This means that U6 is constitutively expressed, at least in these tissues. Certainly, the lower the Ct difference is, the more reliable are the results. To facilitate this, more lab-on-a-chip quantifications and real time assays are needed.
Therefore, we concluded that U6 is a more reliable endogenous control than snoZ30 for miRNA expression profile comparison in the analyzed tissues. The miRNA expression profile comparison is reported in Figure 6.
We then compared our data with those reported in literature [13]. Figure 7 shows the expression values of hsa-miR-26a, hsa-miR-26b and hsa-miR-134 compared with the expression values of some tissue-specific miRNAs (miR-1 for heart and muscle, miR-122a for liver, miR-124a for brain). Expression values of miRNAs are expressed as copies per ng of RNA.
4. Discussion
One of the most difficult problems with miRNAs experimental studies concerns the efficiency of their quantitative and qualitative recovery after total RNA extraction from cells or tissues. In some cases, traditional total RNA extraction protocols are not efficient methods for extracting both high molecular weight (HMW) and low molecular weight (LMW) RNAs. Moreover, conventional quantification methods (i.e., spectrophotometric measurements) giving an overall quantification of total RNA concentration are therefore inadequate for these small RNA species. Northern blotting is the only current viable technique to visualize and evaluate small RNA molecules, particularly miRNAs. However, their low abundance sometimes prevents the use of such techniques. On the other hand, the lab-on-a-chip technology we employed is an alternative and valuable tool for the precise quantification of small RNAs and miRNAs present in total RNA samples.
In this work, to characterize in detail all the LMW RNA species visible in the Bioanalyzer electropherogram (5.8S, 5S, tRNAs, etc.), we ran RNA samples on polyacrylamide gels to clearly separate each of the predominant species present in the small RNA region (Figure 3(a)). We also optimized an elution protocol from excised polyacrylamide bands after gel electrophoresis. Then, we monitored the eluted fractions by means of a lab-on-a-chip technology that is, at present, the only technique that allows for the comparison of small RNA profiles and the exact quantification of these species (Figure 3(b)). Several other species were present under the peak of LMW RNAs, even if it was not possible to distinguish each of them. For these reasons, we applied a mathematical algorithm (the Fourier transform) for the deconvolution of the profile coming from the Bioanalyzer into different Gaussian peaks. After the fitting procedure, we obtained a mathematical model comprising seven peaks that define the overall profile (Figure 3(c)). These peaks range from a minimum of 65 bp to a maximum of 171 bp, which is well in agreement with the results obtained with Agilent automatic peak-peaking that relies on the comparison of RNA species with an RNA ladder included in the running buffer. Interestingly, after this procedure was performed some other non-assigned bands remained undetermined. These species may belong to smaller RNA species like miRNAs. To clarify this, the Small RNA kit, specifically designed for the analysis of miRNA species, was used to study these unidentified species in more detail.
We initially compared the three commonest RNA extraction protocols, the classic acid phenol/guanidine isothiocyanate solution, and two glass fiber filtration protocols to examine and quantify the recovery of HMW and LMW RNA species. The acid phenol/guanidine isothiocyanate solution (Table 1) maximized not only the recovery of HMW RNA fractions but also the LMW RNAs (from 22% to 36%), as previously observed [3]. Then, we assessed by real-time PCR how significantly affected the recovery of miRNA species was as a function of the total RNA extraction protocol used. The expression levels of a miRNA (hsa-miR-21) and a small nucleolar RNA (snoZ30) were evaluated by real-time PCR specific assays [12]. This technique represents an effective alternative to Northern blotting for miRNA detection and quantification. Real-time quantification results (Table 2) demonstrated that the amount of miRNAs correlates better with the amount of LMW RNAs than with total RNA.
Another aspect that we critically analyzed concerned the correlation between LMW RNA species and miRNAs. The presence of LMW RNAs is not always directly correlated with that of miRNAs. In fact, we proved that even with the same extraction protocol the concentration of miRNAs may differ considerably between various preparations (i.e., extraction from different tissues).
As a potential application of our results, we applied our findings to the comparison of some miRNA expression in different tissues. We assessed the expression of hsa-miR-26a, hsa-miR-26b, and has-miR-134 as model miRNAs and two small RNAs (U6 and snoZ30) as controls. The expression profile comparison is correct only if we assume that the endogenous control (U6 or others) has a constant and equal expression in all the considered tissues. The lab-on-a-chip technology we used allowed precise quantification of input miRNAs, ultimately leading to a fine correction of real-time PCR Ct data for eventual variations. Therefore, we obtained a reliable and correct estimate of the relative quantities of miRNAs present in various tissues. Surprisingly, the expression of miR-26a was the highest of all the miRNAs in all the tissues considered (fold change 10 with respect to U6). Although our data are in good agreement with those reported in literature [13] (Figure 7), the high expression values found for miR-26a will deserve future investigations. However, the expression values of miR-26a reported in literature [13] are also very high (and higher than those of miR-26b) and above the values of those tissue-specific miRNAs considered as highly expressed (i.e., miR-1 for muscle, miR-124a for brain, and miR-122a for liver) (see Figure 7).
In conclusion, all the methods we used in this paper allowed us to study accurately the efficiency of extraction protocols, analyze the small (but significant) differences between various preparations, and suggest a methodologically correct method for the comparison of miRNA expression profiles in various tissues.
Acknowledgments
The authors thank editors and reviewers for their useful suggestions and comments to improve the manuscript. This work was supported by Grants of the Italian Ministry of Health. The first and the second authors equally contributed to the work.