Identification of the candidate reference genes from RNA-Seq data
As previously reported by Gerdol et al. [21] Rpl32, Rpl14 and Rpl34 were identified as the first, third and fourth most stable core genes across 51 RNA-seq datasets from multiple adult tissues and experimental conditions in M. galloprovincialis. Rpl32, with a stability score equal to 0.24, achieved a mean expression level of 3061.45 TPMs; Rpl14, with a stability score equal to 0.25, achieved a mean expression level of 2747.49 TPMs; Rpl34, with a stability score equal to 0.27, achieved a mean expression level of 2288.28 TPMs (Figure S1). These three candidate genes were confirmed to lack paralogous gene copies in the M. galloprovincialis genome, and in silico primer check ruled out the possibility of non-specific amplification of other genomic regions or transcribed RNAs.
Expression stability of candidate reference genes
In this work, we analyzed ten genes, considering three novel candidates selected with a transcriptomic approach (Rpl14, Rpl32, Rpl34) and seven genes commonly used as references (Act, Cyp-A, Ef1α, Gapdh, 18S, 28S, Rps4), verifying their expression stability in various mussel tissues. The range of cycle threshold (Ct) values was quite similar among samples, with 18S and 28S exhibiting the highest range of variation in Ct in all tissues analyzed (Fig. 1, S2).
An overview of the results obtained from the analysis of the selected reference genes using three different algorithms is presented in Fig. 2. In the gills, Rps4 was identified as the most stable gene by geNorm evaluation, whereas Act and Cyp-A were selected by the NormFinder and BestKeeper, respectively. However, all these three genes were generally considered stable by the three algorithms but with a different ranking (Table S1). Interestingly, geNorm suggested the use of Gapdh/Rpl34 as the best combination of reference genes, as both showed higher M values than the other analyzed genes (Fig. 2A, Table S1).
In the digestive gland, Gapdh and Rps4 showed the lowest (and identical) M value according to geNorm; nevertheless, the algorithm suggested Rps4/Rpl14 as the best possible combination since the latter gene had an M value slightly higher than Gapdh (Table S2). Gapdh has ranked first also in the NormFinder analysis (Table S2). In addition, both Rps4 and Rpl14 displayed good stability values. Conversely, Gapdh, Rps4 and Rpl14 were found to be beyond the recommended value (SD < 1) in the BestKeeper analysis, while Cyp-A ranked first (Table S2). It is important to note that Cyp-A was also identified as a stable gene by the geNorm and NormFinder analysis (Table S2, Fig. 2B).
The comparison among the different algorithms produced less consistent results in the gonads, with Gapdh resulting as the most stable gene by geNorm, Rps4 by NormFinder and Rpl34 by BestKeeper (Table S4). In addition, the best combination, according to geNorm was Cyp-A/Gapdh. All the named genes were generally considered stable by the three algorithms, except for Rps4 that showed an SD value much higher than recommendations (SD = 1.25) by BestKeeper (Table S4, Fig. 2C).
In the mantle, Cyp-A showed the highest stability according to both geNorm and BestKeeper, in contrast with NormFinder, which indicated Act as the most stable gene (Table S3). Nevertheless, Cyp-A and Act were considered stable by the three analyses, even though they occupied different positions in the rankings. Furthermore, geNorm suggested Cyp-A/Rpl14 as the best combination of reference genes (Table S3, Fig. 2D).
In the foot, Gapdh, Act and Cyp-A had the best stability values in geNorm, NormFinder and BestKeeper ranking, respectively (Table S5). The best pair of reference genes suggested by geNorm was Rpl14/Rpl32. As previously discussed for the other tissues, all these genes were considered relatively stable from the different algorithms, with the exception of act, which had an SD value higher than 1 (SD = 1.29) (Table S5, Fig. 2E).
Finally, we reanalyzed the data by considering all the tissues. As previously mentioned for the analysis of single tissues, 18S and 28S were the most variable candidate reference genes across all the tissues (Fig. 3). NormFinder defined Act and Cyp-A as the best gene pair combination.
Expression level of target genes
We validated the bonafide selection of the previously described reference genes for the interpretation expression data of a set of GOI. To this aim, we analyzed the expression of four genes differentially expressed in adult tissues, selected according to their strong specific expression inferred from RNA-seq data in the mantle and digestive gland, respectively. In detail, MGAL10A005692, specifically expressed at high level (averaging ~ 3900 TPM) in the mantle, encodes a low-complexity protein, which bears distant homology with SCO-spondin in the C-terminal region. The second selected mantle-specific GOI, MGAL10A087091, was expressed at similar levels (averaging 3600 TPM) and was orthologous with M. coruscus MUSP-3, previously reported to be involved in shell mineralization [28].
On the other hand, MGAL10A040115, a meprin A-like metalloprotease, was one of the most distinctive digestive gland-specific genes in mussel, reaching expression values as high as 100,000 TPM in some samples. The second selected digestive gland GOI, which also displayed very high mean expression values (32,000 TPM) was MGAL10A054097, encoding a protein with multiple ShKT domains.
As reported in Fig. 4, the expression profile of MGAL10A005692 and MLGAL10A087091 confirmed the high mantle-specificity indicated by transcriptome data. In particular, the Ct values were higher than 30 in all the tissues other than mantle and the expression of MGAL10A087091 was almost undetectable in gonad and gills (data not shown). The expression profiles of MGAL10A040115 and MGAL10A054097 are reported in Fig. 5. The specificity of these genes for the digestive gland was strongly confirmed, as in all the tissues other than digestive gland the Ct values were higher than 30 (data not shown).
We tested a few other specific GOIs in the digestive gland (MGAL10A065261, MGAL10A066134), in the mantle (MGAL10A006232, MGAL10A059842, MGAL10A079192) and in the foot (MGAL10A033477, MGAL10A064866, MGAL10A070506). However, possibly owing to the low sequence complexity of these gene targets (which often encoded highly repetitive proteins), the efficiency of amplification observed by RT-qPCR data in these cases was very low (data not shown).
For all analyzed GOIs, no significant differences could be evidenced between the use of Act as a single reference gene and the use of the Act-Cyp-A pair (Figs. 4 and 5). As reported in Figure S3, the use of the last ranking genes, such as 18S and the couple 18S and 28S still allow to detect the tissue specific of the GOIs (i.e MLGAL10A087091 and MGAL10A040115, mantle and digestive gland, respectively) but with an higher standard error.