mRNA N⁶-methyladenosine methylation of postnatal liver development in pig

<div><p>N⁶-methyladenosine (m⁶A) is a ubiquitous reversible epigenetic RNA modification that plays an important role in the regulation of post-transcriptional protein coding gene expression. Liver is a vital organ and plays a major role in metabolism with numerous functions. Information concerning the dynamic patterns of mRNA m⁶A methylation during postnatal development of liver has been long overdue and elucidation of this information will benefit for further deciphering a multitude of functional outcomes of mRNA m⁶A methylation. Here, we profile transcriptome-wide m⁶A in porcine liver at three developmental stages: newborn (0 day), suckling (21 days) and adult (2 years). About 33% of transcribed genes were modified by m⁶A, with 1.33 to 1.42 m⁶A peaks per modified gene. m⁶A was distributed predominantly around stop codons. The consensus motif sequence RRm⁶ACH was observed in 78.90% of m⁶A peaks. A negative correlation (average Pearson’s <i>r</i> = -0.45, <i>P</i> < 10⁻¹⁶) was found between levels of m⁶A methylation and gene expression. Functional enrichment analysis of genes consistently modified by m⁶A methylation at all three stages showed genes relevant to important functions, including regulation of growth and development, regulation of metabolic processes and protein catabolic processes. Genes with higher m⁶A methylation and lower expression levels at any particular stage were associated with the biological processes required for or unique to that stage. We suggest that differential m⁶A methylation may be important for the regulation of nutrient metabolism in porcine liver.</p></div

Time-series modules and co-expression network of lncRNAs and protein-coding genes.

<p>(A) Time-series modules of protein-coding genes and lncRNAs. The top panel shows protein-coding genes and the second panel shows lncRNAs. Numbers in the top left corner indicate module number. Numbers in lower left corners indicate numbers of protein-coding genes or lncRNAs in each module. The same color was used to represent each cluster. Functional categories of genes in green (B) and red modules (C). Benjamini adjusted <i>P</i> values were transformed by ‒log₁₀. (D) Heat map showing the largest two co-expression networks of protein-coding genes. Values represent log₂(FPKM+1) of each gene in each sample minus the mean value of each gene across all samples.</p

Expression profile and PCA of protein-coding genes.

<p>(A) Heat map showing the expression profile of protein-coding genes. The top panel is the tree constructed by Pearson correlation. (B) Two-way PCA plot of protein-coding genes based on expression profile.</p

Temporal expression profiles of protein-coding genes and lncRNAs.

<p>(A) Dynamic changes in expression profiles of protein-coding genes and lncRNAs. The top panel shows protein-coding genes and the bottom panel shows lncRNAs. Values represent the pairwise Pearson correlation. Correlation between every two samples was calculated by log₂-transformed (FPKM+1) gene expression values. Three main expression patterns can be distinguished. (B) Distributions of Shannon entropy-based temporal specificity scores were calculated for distinct classes of lncRNAs and protein-coding genes.</p

Differentially expressed protein-coding genes and lncRNAs, and PCA of PSI values.

<p>Venn diagram of common differentially expressed protein-coding genes (A) and lncRNAs (B) in five developmental stages. (C) Dynamic expression profiles of <i>CP</i> and TU78568. (D) Two-way PCA plot of protein-coding genes based on PSI values.</p

Functions of genes with higher m⁶A peak enrichment and lower expression.

<p>Functions of genes with higher m⁶A peak enrichment and lower expression.</p

Overview of m⁶A methylation in porcine liver.

<p>(A) Venn diagram showing the overlap of m⁶A peaks in newborn (8,855), suckling (7,350) and adult (7,961). There are 5,848 common peaks among the three stages, which with ≥ 50% length overlap between stages. (B) Venn diagram showing the overlap of m⁶A modified genes. Respectively, 4,676 genes in newborn, 4,103 in suckling and 4,339 in adult were m⁶A methylated. For all three stages, 3,481 genes were consistently modified. (C) Proportion of genes containing variant numbers of m⁶A peaks. Majority of modified genes (74.60%) contain one or two m⁶A peaks, while the rest contains more. (D) Sequence logo representing the most common consensus motif (RRm⁶ACH) in the m⁶A peaks. The consensus sequence was detected by DREME (version: 4.10.2), using the 101 nucleotides centered on the summits of called original narrow peaks.</p

Relationship between m⁶A methylation and expression of modified genes.

<p>(A) Fraction of genes with m⁶A peaks in each of the segments as a function of expression level. Most of the modified genes were expressed at moderate levels. Genes expressed at the two extremes were less methylated. (B) Plot of m⁶A peak enrichment and mRNA abundance in the three stages. Obvious negative correlation between m⁶A peak enrichment and modified mRNA abundance was found (Pearson’s <i>r</i> = -0.47 to -0.42, <i>P</i> < 10⁻¹⁶). Lines represent the linear trend for the obtained values.</p

Number of genes showing differential transcript levels and differential m⁶A methylation.

<p>Number of genes showing differential transcript levels and differential m⁶A methylation.</p

m⁶A enrichment and gene expression profile of <i>GATM</i> in three stages.

<p>Opposite trends of the m⁶A methylation level (left panel) and gene expression level (right panel) of <i>GATM</i> are shown. Gene expression level is presented by the accumulation of input reads.</p