Genetic regulation and heritability of miRNA and mRNA expression link to phosphorus utilization and gut microbiome

Improved utilization of phytates and mineral phosphorus (P) in monogastric animals contributes significantly to preserving the finite resource of mineral P and mitigating environmental pollution. In order to identify pathways and to prioritize candidate genes related to P utilization (PU), the genomic heritability of 77 and 80 trait-dependent expressed miRNAs and mRNAs in 482 Japanese quail were estimated and eQTL (expression quantitative trait loci) were detected. In total, 104 miR-eQTL (microRNA expression quantitative traits loci) were associated with SNP markers (false discovery rate less than 10%) including 41 eQTL of eight miRNAs. Similarly, 944 mRNA-eQTL were identified at the 5% False discovery rate threshold, with 573 being cis-eQTL of 36 mRNAs. High heritabilities of miRNA and mRNA expression coincide with highly significant eQTL. Integration of phenotypic data with transcriptome and microbiome data of the same animals revealed genetic regulated mRNA and miRNA transcripts (SMAD3, CAV1, ENNPP6, ATP2B4, miR-148a-3p, miR-146b-5p, miR-16-5p, miR-194, miR-215-5p, miR-199-3p, miR-1388a-3p) and microbes (Candidatus Arthromitus, Enterococcus) that are associated with PU. The results reveal novel insights into the role of mRNAs and miRNAs in host gut tissue functions, which are involved in PU and other related traits, in terms of the genetic regulation and inheritance of their expression and in association with microbiota components.


Introduction
Phosphorus (P) is an essential mineral and a limited resource. The excessive use of mineral P in agricultural production has an environmental impact. Consequently, vigilant P management is required. The bulk of P in animal feeds originates from plant seeds. However, up to 80% of plant P is in the form of phytates, which cannot be efficiently used by monogastrics due to the lack or scarcity of endogenous phytase in their digestive tract. Microbial phosphatases may contribute to the cleavage of P from phytates in the gut of monogastric animals. Gut properties, the microbiome composition and interactions between the gut tissue and microbiota play a significant role for digestive capacity and need to be understood for improvement of phosphorus utilization (PU). The most promising molecules mediating between gut microbiota and host are microRNAs (miRNAs), which are conserved between species and can regulate gene expression across species. We previously identified differentially abundant miRNAs, mRNAs and microbiota in the gut of Japanese quail representing extremes for PU [1,2].
MiRNAs are small endogenous non-coding molecules ranging from 18 to 24 nt that control target transcripts by inducing transcript cleavage, degradation, destabilization and repression of translation, thereby modulating protein levels. MiRNAs contain a 6-8 nt section known as seed sequence [3]. This sequence is important for targeting mRNA sites via base-pair complementarity, typically in the 3 0 -untranslated region [4,5] but also in the coding sequence [6]. Since a single miRNA can regulate the expression of several genes and a single gene can be influenced by several miRNAs, it was assumed that more than 50% of the genes are controlled by miRNA [7]. Moreover, miRNA-expression profiles have been associated with many complex traits and diseases [8][9][10]. Transcript abundances are endogenous traits that represent an important contribution to the expression of complex traits and that depend on genetic polymorphisms (i.e. they are heritable). Genetic variation associated with variation of expression levels can be investigated through expression-QTL (eQTL) analysis. Many eQTL mapping studies reveal genetically regulated mRNA transcripts in various species and tissues. However, miRNA-eQTL studies are still sparse [11][12][13][14][15][16]. The heritability of transcript abundances of various miRNAs and mRNAs is one criterion to prioritize candidates for further analyses [17].
Our previous study in an experimental F2 population of Japanese quail showed genetic parameters of PU, including a low heritability of 0.14 and a causal structure indicating that a high PU is associated with a low gain per feed and a high body weight gain [18]. In addition, PU was shown to have a strong genetic and phenotypic correlation with bone ash, a trait with a much higher heritability of 0.23-0.34 [19]. New molecular pathways affecting PU were identified by detecting differentially expressed miRNAs and mRNAs in the ileum of PU-divergent Japanese quail, and a biomarker panel associated with PU was provided [1,2]. Further studies looked at microbial gut profiles [2,20,21]. A microbial core community showed shifts in abundances associated with the PU. The list of miRNAs and mRNAs that showed differential expression in ileum tissue of Japanese quail with divergent PU serves as a basis to select transcripts to be used here in a study of a larger number of animals. The identification of SNPs associated with miRNA and mRNA abundance in the ileum has the potential to aid in the understanding of the role of gut cell functions in PU and related traits. We therefore estimated the heritability and examined the genetic regulation of ileum miRNA and mRNA expression and its consequences on mineral utilization including PU, calcium utilization (CaU), bone ash traits (femur ash absolute, F ash abs; femur ash relative, F ash rel; tibial ash absolute, T ash abs; tibial ash relative, T ash rel) and performance traits. At the same time, we considered the association of gut microbiota data obtained from the same animals with transcripts and phenotypes.

Results
After having identified differentially expressed genes, the aim of this study was to identify genetic variants impacting both miRNA and mRNA expression levels and to estimate respective heritabilities. The top differentially expressed transcripts including 77 miRNAs and 80 mRNAs were selected for genomic heritability and genetic regulation study in 482 Japanese quails.
royalsocietypublishing.org/journal/rsob Open Biol. 11: 200182 associated with CTDSP1 or PPP6C, respectively. In other case, SNPs associated with miR-383-5p in cis were also associated with transcript abundances of ETNPPL, KDR and LOC107312819 located on different chromosomes.

Discussion
Extremes for PU traits out of an F2 population of Japanese quail from our previous study reveal molecular interactions occurring in the gut of quail and provided a number of biomarker panels associated with the PU [2]. These biomarkers, which include differentially expressed mRNAs, miRNAs and their targets, are not only linked to PU traits, but also to microbiota enriched in extremes for PU traits. In addition, the mRNA, miRNAs and its targets were enriched in functional pathways involved in energy metabolism, cell proliferation, phosphate or bone metabolism, as well as immune pathways [1,2]. Here we further focused on selected mRNA and miRNA transcripts by estimating heritabilities of expression levels and detecting eQTLs. The molecular traits affect the phenotype at the organismal level including traits related to mineral homeostasis such as PU, CaU and bone ash traits or performance traits such as BWG, FCR and FI.  royalsocietypublishing.org/journal/rsob Open Biol. 11: 200182 In fact, a previous study showed a strong genetic and phenotypic correlation between PU and bone ash traits which are alternative traits reflecting PU [22]. In general, a higher number of eQTL were identified for mRNA than for miRNA. Moreover, the average genomic heritability estimates were slightly higher for mRNA compared to miRNA. In addition, eQTLs of miRNAs or mRNAs in cis dominated compared to trans eQTLs in terms of number and significance, which is consistent with other studies [16,23]. Transcript abundances of a number of miRNAs and mRNAs proved to be heritable molecular features. The majority of these also had significant eQTL which is in line with a previous study [22]. For example, the expression of miR-383-5p was highly heritable and a cis-miR-eQTL existed.
We found a majority of common or shared eQTLs affecting both miRNA and mRNA expression, which was also in line with previously reported results [24]. However, this previous study found that the majority of shared eQTLs affected miRNA and host expression differently [24]. In addition, a number of miRNAs and mRNAs, whose expression correlates with characteristics of the P metabolism, whose transcript abundances are heritable, and for which eQTL have thus been identified, also show functional connections with metabolic pathways of P homoeostasis and/or with the microbiota of relevance for PU (figure 9).
These interrelationships complement results from GWAS and facilitate their interpretation. Remarkably, most of miRNA and mRNA that we previously identified as differentially expressed between Japanese quail divergent for PU were confirmed in this study using a large number of animals representing the whole phenotypic distribution. In addition, the study provides a narrowed list of candidate transcripts mediating interaction between host gut microbiomes and PU characteristics. As shown in our previous study with the same F2 population of Japanese quail, intestinal microbiota are related to the above-mentioned phenotypes [2,20] and specific microbiota have shifted abundance in parallel with PU. This includes SMAD3, CAV1, ENNPP6, ATP2B4, miR-148a-3p, miR-146b-5p, miR-16, miR-194, miR-215, miR-199-3p and miR-1388a-3p. In fact, out of the 29 mRNAs significantly correlated with PU, 13 were identified as enriched with the gene ontology term 'regulation of phosphorus metabolic process' (GO: 0051174), including SMAD3, CAV1, PPP1R1A, ATP2B4 and PKD1. , calcium utilization (CaU)), bone ash traits (femur ash absolute (F ash abs), femur ash relative (F ash rel), tibial ash absolute (T ash abs), tibial ash relative (T ash rel)) and body weight gain (BWG), feed conversion rate (FCR) and feed intake (FI)) after correction of all covariances (method). The colour bar indicates positive correlation (red) and negative correlation (blue).
royalsocietypublishing.org/journal/rsob Open Biol. 11: 200182 ATP2B4 and PKD1 are involved in intestinal calcium uptake and, thus, are closely related to P homeostasis in terms of a constant cellular Ca : P ratio. ATP2B4 is located in the plasma membrane and promotes intestinal Ca transport [25]. For PKD1 knock-out mice, the abundance of this gene was shown to affect the renal expression of several genes involved in mineral homeostasis, including ATP2B4, and consequently leads to altered serum mineral levels [26]. SMAD3 expression was reported as crucial for osteogenesis in studies of stem cell differentiation into cells with an  . Pearson correlation coefficients and p-values between residuals of mRNAs expression (ct value) and phenotypes ( phosphorus utilization (PU), calcium utilization (CaU)), bone ash traits (femur ash absolute (F ash abs), femur ash relative (F ash rel), tibial ash absolute (T ash abs), tibial ash relative (T ash rel)) and body weight gain (BWG), feed conversion rate (FCR) and feed intake (FI)) after correction of all covariances (method). The colour bar indicates positive correlation (red) and negative correlation (blue). royalsocietypublishing.org/journal/rsob Open Biol. 11: 200182 osteoblastic phenotype [27]. The rate of osteogenic differentiation was decreased when SMAD3 signalling was impaired [27]. In addition, SMAD3 cis-acting eQTL operating in primary osteoarthritis and in the aneurysms and osteoarthritis syndrome was reported [28]. In addition, SMAD3 directly mediates TGF-β signalling which is active and regulates the cellular functions of intestinal epithelial cell [29]. Smad3-/-mice showed defective transforming growth factor β-signalling and associated with disease conditions of the gut [30,31]. In our study of Japanese quail, SMAD3 showed a cis-mRNA-eQTL ( p-value = 1.5 × 10 −15 ) located on chromosome 10 (position 171 kb) with the highest association of SNP_ID18249. The expression of this gene was also significantly correlated with PU, CaU, FI, FCR and BWG. The gene was revealed in our previous study and provides an intestinally expressed candidate molecule for divergent PU in quails [2]. Putative mechanisms are the interaction of SMAD3 with the vitamin D receptor and the initiating of effects mediated by vitamin D3 [32]. Moreover, SMAD3 showed a significantly negative correlation with miR-191-5p (r = −0.23, p < 2 × 10 −7 ) followed by miR-215-5p (r = −0.15, p < 7 × 10 −4 ).
ENPP6 and miR-16-5p were significantly correlated with PU, CaU, FI, FCR and BWG. A cis-eQTL of miR-16-5p mapped to chromosome 1 ( position 151.6 kb) with two SNPs (id10120 and id07989) (chromosome 1 position 152 kb) being associated at p = 2.0 × 10 −4 . A previous study reported that insertion mutation disrupts genesis of miR-16 and causes increased body weight in domesticated chicken [33]. ENPP6 has also a cis-eQTL located on chromosome 4 ( position 35 kb) with SNP id17811. Both miR-16-5p and ENPP6 were reported in our previous study to link to PU [2]. ENPP6 activities mediate bone mineralization which is a key process in the formation of bone [34]. In addition, ENPP6 was found positively correlated with Candidatus Arthromitus while miR-16-5p positively correlated to Enterobacter.
The link between the host and gut microbiota can be mediated by miRNA. Due to their strong evolutionary conservation, faecal miRNA-mediated inter-species gene regulation facilitates host control of the gut microbiota [35]. In other case, microbiota-derived metabolites may influence the host's cellular functions. Recently, a study reported that microbiota-derived nitric oxide directly alters the host's Argonaute family protein activity, and consequently altered miRNAs and gene expression of the host [36]. CAV1 belong to nitric oxide (NO) signalling and negatively regulates NO  Figure 8. Pearson correlation coefficients and p-value between residuals of microbiota and phenotype ( phosphorus utilization (PU), calcium utilization (CaU)), bone ash traits (femur ash absolute (F ash abs), femur ash relative (F ash rel), tibial ash absolute (T ash abs), tibial ash relative (T ash rel)) and body weight gain (BWG), feed conversion rate (FCR) and feed intake (FI)) after correction of all covariances (method). The colour bar indicates positive correlation (red) and negative correlation (blue).
royalsocietypublishing.org/journal/rsob Open Biol. 11: 200182 synthase activity [37,38]. The reduction in NO concentration in the gut of mice after treatment with a polysaccharide from Spirulina platensis was accompanied by an increase in the abundance of beneficial bacteria, including Candidatus Arthromitus [39]. Candidatus Arthromitus belongs to segmented filamentous bacteria (SFB) and selectively attach to the ileal epithelium of the host [40]. SFB genomes contain genes predicted for two catalases, a peroxidase and an arginase, which might limit nitric oxide production through catabolism of arginine [41,42]. The previous study showed the role of SFBs (in particular Candidatus Arthromitus) in the ileum of turkeys that were related to weight gain which is also in line with our study in Japanese quail [43]. The previous study of the same animals demonstrated the impact of the microbiota composition on quantitative trait variation in particular PU, CAU, FI, FCR and BWG [21]. The integration of analyses of microbiota abundances with the genetic variation of the host, transcript levels and complex traits reveals new relationships. For instance, the microbiome-wide association analyses to study the relationship of abundances of microbiota with other traits revealed association of some genera including Candidatus Arthromitus with multiple traits [21]. In our study, CAV1 ct value was negatively correlated to miR-1388a-5p and miR-1388-3p expression and positively correlated with Candidatus Arthromitus abundance. We found CAV1, miR-1388a-5p, miR-1388-3p and Candidatus Arthromitus being linked to the phenotype of BWG, FI, FCR and PU. We also found that the miR-1388a-5p and miR-1388-3p were highly heritable.
The correlations shown here provide further evidence for the complex interactions that exist between microbiota, mRNA and miRNA expression, and phenotypes, which in line with detailed previously described connections. In particular, we found evidence that the interaction of miRNAs and mRNAs and the gut microbiota play an important role in PU and other related phenotypes. In addition, this study demonstrated genetic regulation by eQTL analysis and the heritability of expression of both miRNAs and mRNAs correlated with mineral utilization and other related traits and associated with variable microbiota abundances. The study also revealed that approaches to improve a complex trait such as PU, which allow for resource-efficient and environmentally friendly animal husbandry, need to consider multiple levels of the genotype-phenotype map.

Experimental design and samples selection
Japanese quails (Coturnix japonica) used in this study were originated from a previous study [18,19]. The experiment  Figure 9. Network of top 20 links between miRNA, mRNA, microbiota and phenotypes including phosphorus utilization (PU), calcium utilization (CaU), bone ash traits (femur ash absolute (F ash abs), femur ash relative (F ash rel), tibial ash absolute (T ash abs), tibial ash relative (T ash rel)) and body weight gain (BWG) and feed conversion rate (FCR) and feed intake (FI) traits based on Pearson correlation coefficients.
royalsocietypublishing.org/journal/rsob Open Biol. 11: 200182 was conducted in accordance with the German Animal Welfare Legislation approved by the Animal Welfare Commissioner of the University of Hohenheim. An F2 design using two Japanese quail lines divergently selected on social reinstatement behaviour was established with 12 males from line A and 12 females from line B of the F0 generation. A total of 17 roosters and 34 hens were randomly selected from the F1 birds to generate F2 animals [18]. In order to let the birds express their full genetic potential of PU, the F2 animals were fed with a low-P diet (4.0 g P kg −1 dry matter) without a mineral P supplement or phytase. In total, 482 birds of the F2 quail population were used. All the traits used in this study including PU, CaU, were descripted in detail by Beck et al. [18]. , respectively, and FC was the individual feed intake measured over a 5day collection period [g] [18]. The corresponding calculation applies for CaU.

RNA extraction and quantitative PCR
After feeding the low-P diet for a period of 10 days and when the animals were 15 days old, the ileum samples were dissected, cut open and rinsed with a sterilized saline buffer to remove digesta residue. The intestinal epithelium was scraped using sterilized razor blades and were immediately submerged in RNAlater solution (Sigma) and stored at −80°C until RNA extraction. Total RNA was extracted from approximately 50 mg sample using TRIzol Reagent (Invitrogen) and the RNeasy Mini kit (Qiagen) with DNaseI treatment according to manufacturer's recommendations for mRNA. A 2100 Bioanalyzer Instrument (Agilent, Santa Clara, CA) was used for RNA quality assessment. The cDNA synthesis of miRNA was performed according to a previous study [44]. In brief, 100 ng of total RNA were poly(A) tailed and reverse transcribed using 1 unit of poly(A) polymerase 1 µM (BioLab), RT-primers (CAGGTCCAGTTTTTTTTTTTTTTTVN where V is A, C and G and N is A, C, G and T), 0.1 mM of NTPs, 100 units of MuLV reverse transcriptase (Invitrogen). The reaction was incubated at 42°C for 1 h followed by 95°C to inactivate the enzyme. For cDNA synthesis from mRNAs 200 ng mRNA was mixed with 1 µl Reverse Transcription Master Mix (Fluidigm PN 100-6297) in 5 µl volume. The reaction was incubated at 25°C for 5 min, 42°C for 30 min followed by 85°C for 5 min. The cDNA was used for further quantitative PCR (qPCR). In total, 77 miRNA and 80 mRNA from 482 F2-quail samples were used for qPCR with the Fluidigm BioMark HD System. Specific target amplification (STA) was done per manufacturer's recommendations. Pre-amplification sample mixtures were prepared using PreAmp Master Mix (Fluidigm PN 1005581) containing 1.25 µl of cDNA, 1 µl PreAmp Master Mix, and 0.5 µl Pooled Delta Gene Assay Mix (500 nM) containing DNA suspension buffer and primers mixes (electronic supplementary material, file S1 for miRNA and electronic supplementary material, file S2 for mRNA) in 5 µl total volume. The pre-amplification reaction was incubated at 95°C for 2 min, followed by 10 cycles at 95°C for 15 s and 60°C for 4 min. The pre-amplification reaction was cleaned up using exonuclease I, followed by 10 × dilution of STA with DNA suspension buffer (TEKnova, PN T0221). Fluidigm quantitative measurement runs were carried out with 96.96 dynamic arrays (Fluidigm Corporation, CA, USA) according to manufactures instructions. In brief, 2.5 µl of 2 × SsoFast Evagreen Supermix with Low ROX, 0.25 µl 20 × sample-loading reagent, and 2.25 µl of treated samples were prepared. Separately, an assay mixture was prepared for each primer pair and this included 2.25 µl of DNA Suspension buffer, 0.25 µl of 100 µM forward and reverse primer and 2.5 µl of 2 × assay-loading reagent. The dynamic arrays were first primed with control line fluid and then loaded with the sample and assay mixtures via the appropriate inlets using an IFC controller. The array chips were placed in the BioMark Instrument for PCR at 95°C for 10 min, followed by 30 cycles at 95°C for 15 s and 60°C for 1 min. The data were analysed with real-time PCR analysis software in the BioMark HD instrument (Fluidigm Corporation, San Francisco, CA). The internal controls of cel-miR-39-3p, 18S, SNORD21 were used for miRNA and six housekeeping gene (RAB35, RPS11, ACTB, AGFG1, UB2A and UBE2D1) were used for mRNA.

eQTL detection
Genotyping data were generated in a parallel genome-wide association and linkage study of PU and related traits [45]. In brief, genotyping was performed using costumer's Illumina iSelect Chips for 5388 SNPs genotyped. After normalization and filtering, 3469 SNPs of the 482 individuals were used. Based on the estimated linkage map of these SNPs markers, this map was updated and compared with the new reference genome of the Japanese quail (Coturnix japonica 2.0) [45]. The 77 miRNA and 80 mRNA Ct values and SNPs of 482 individuals were subjected to a mixed-model analysis of variance using JMP Genomics (Proc Mixed; SAS Institute, Rockville, MD). An identity by descent (IBD) matrix giving the relationships between individuals was used as a random effect. Additionally, genotype, gender and the batch of RT-PCR were used as fixed effects. Ct values of 18S, SNORD21 and cel-miR-39-3p and the housekeeping genes were used as covariates for analyses of miRNAs and mRNAs, respectively. To correct for multiple testing, a FDR at 10% was used for miRNA and 5% for mRNA. We defined an eQTL as 'cis' if an associated SNP was located within an area less than 10 Mb from the corresponding gene. All other eQTL were considered as 'trans'.

Genomic heritability estimation
Genomic heritability represents the proportion of genetic variance explained by SNPs to the phenotypic variant (i.e. miRNA or mRNA expression), and was calculated using the JMP Genomics. The estimated genetic influence on traits was based on the SNP data matrix instead of using a formal experimental design [46]. The measured SNP-level variation is used to estimate the genetic similarity between individuals, and this estimated genetic similarity is compared to the phenotypic similarity to produce a heritability estimate. To estimate the heritability of the expression for a specific miRNA or mRNA, sex and technical covariates were included as a fixed effect, the genetic similarity matrix between individuals was first computed as identityby-descent of each pair for the k-matrix (SNPs) used as the random effect.
royalsocietypublishing.org/journal/rsob Open Biol. 11: 200182 4.5. Trait-correlated ileum transcripts expression and gut microbiota Operational taxonomic units (OTUs) deduced from 16S rRNA sequencing of the same animals was obtained from a recent study [20]. Initially, OTUs were assigned to taxa at the genus level and OTU counts belonging to the same genera were summarized. Moreover, the dataset was filtered so that only taxa with more than one observation in at least half of the samples were considered. In total 91 gut microbiota genera from 411 individuals were used for the correlation study. After adjustment for the effect of either miRNA-or mRNA expression, gut microbiota abundances or phenotypic traits as described above, the residuals were retained for further analysis. The correlation of miRNA or mRNA Ct values as well as gut microbiota with PU and other related traits were estimated as Pearson coefficients and corrected for multiple comparisons by calculating the FDR [47]. Networks of mRNA, miRNA, gut microbiota and the phenotype were visualized with cytoscpe (http:// cytoscape.org).
Ethics. The experiment was conducted in accordance with the German Animal Welfare Legislation approved by the Animal Welfare Commissioner of the University of Hohenheim (approval no. S371/13TE).
Data accessibility. All data are available in the electronic supplementary material.
Authors' contribution. S.P. and K.W. contributed to the conception and