- 1Key Laboratory of Swine Genetics and Breeding of Ministry of Agriculture & Key Laboratory of Agriculture Animal Genetics, Breeding and Reproduction of Ministry of Education, College of Animal Science, Huazhong Agricultural University, Wuhan, China
- 2Key Lab of Swine Healthy Breeding of Ministry of Agriculture and Rural Affairs, Guangxi Yangxiang Co., Ltd., Guigang, China
- 3Hubei Key Laboratory of Agricultural Bioinformatics, College of Informatics, Huazhong Agricultural University, Wuhan, China
Somatic cloning has had a significant impact on the life sciences and is important in a variety of processes, including medical research and animal production. However, the application of somatic cloning has been limited due to its low success rate. Therefore, potential epigenetic variations between cloned and donor animals are still unclear. DNA methylation, one of the factors which is responsible for phenotypic differences in animals, is a commonly researched topic in epigenetic studies of mammals. To investigate the epigenetic variations between cloned and donor animals, we selected blood and ear fibroblasts of a donor pig and a cloned pig to perform whole-genome bisulfite sequencing (WGBS). A total of 215 and 707 differential methylation genes (DMGs) were identified in blood and ear fibroblasts, respectively. Functional annotation revealed that DMGs are enriched in many pathways, including T/B or natural killer (NK) cell differentiation, oocyte maturation, embryonic development, and reproductive hormone secretion. Moreover, 22 DMGs in the blood and 75 in the ear were associated with immune responses (e.g., CD244, CDK6, CD5, CD2, CD83, and CDC7). We also found that 18 DMGs in blood and 53 in ear fibroblasts were involved in reproduction. Understanding the expression patterns of DMGs, especially in relation to immune responses and reproduction, will reveal insights that will aid the advancement of future somatic cloning techniques in swine.
Introduction
Somatic cell nuclear transplantation (SCNT) technology has been widely studied in basic research in medical science and livestock production. SCNT has been used to successfully clone different mammals (Li et al., 2006; Ogura et al., 2013). SCNT can produce genetically identical animals, and the production performance of the cloned animal can be similar to that of the donor. Consequently, SCNT may be an effective method for the production of indigenous target animals for use in epigenetic research (Adachi et al., 2014; Keefer, 2015; Lee et al., 2018). However, somatic cell cloning has a low efficiency success rate, which may be affected by somatic cell type and oocyte activation status (Akagi et al., 2014). Although several methods have been used to improve cloning efficiency, abnormal phenotypes that occur due to cloning continue to limit the utilization of this technology in industry. Somatic cell cloning inefficiency and phenotypic abnormalities are probably related to nuclear reprogramming, which may lead to a shortened life span for the cloned animals as well as diseases including, but not limited to severe pneumonia and thymic hypertrophy (Rideout et al., 2001; Ogura et al., 2002; Shimozawa et al., 2006). All of the aforementioned phenomena suggest that cloning may lead to low immunity in the cloned animals. However, several studies have reported that cloned animals had similar characteristics to the donor animals. For example, it has been reported that cloned Landrace and Jinhua sows grew normally and had similar carcass characteristics and reproductive performance compared to their offering and/or control group (Shibata et al., 2006; Adachi et al., 2014). Therefore, cloning holds great promise as a tool to produce elite breeding livestock. However, there is a lack of comparative studies on the reproductive performance and methylation of cloned and donor pigs which means that phenotypic and epigenetic differences between cloned and donor pigs are as yet unknown.
DNA methylation, the most common epigenetic modification in eukaryotes, alters the activity of non-coding elements such as introns, resulting in changes in the mechanism of action of DNA-proteins, and subsequent alterations in mammalian phenotypes (Jones and Takai, 2001). Mammalian DNA methylation mainly occurs in 5-methylcytosine and is associated with gene silencing (Li and Zhang, 2014), which is essential for normal mammalian development and normal biological phenotypes (Auclair and Weber, 2012; Bogdanovic and Gomez-Skarmeta, 2014). Similarly, epigenetic differences in monozygotic twins are an important cause of autoimmune diseases such as rheumatoid arthritis and type I diabetes in one of the twins (Hannon et al., 2018). A previous study has shown that, sows with different levels of DNA methylation have a difference litter size and different methylation gene mRNA expression levels (Hwang et al., 2017b). Sows cloned using ear fibroblasts display DNA methylation variation during the cloning process, and this variation leads to a short lifespan and organ dysplasia of cloned sows from the same pig (Shen et al., 2012). To investigate whether cloning will lead to a variation in DNA methylation, and subsequent changes in immune responses and reproductive performance in pigs. This work aims to describe the characterization of these epigenetic changes in the examined tissues of a cloned sow, in order to provide some insights for genetic studies of cloning.
Materials and Methods
Animals and Sample Collection
We cloned a normal female piglet from a fourth-parity Danish Yorkshire sow with an average litter size ≥18. Next, we collected the piglet's ear fibroblasts on the day of delivery for SCNT. Fibroblasts were separated by trypsin and collagenase for cell culture. The surrogate sow was normal, four-litter, gentle, and manageable. This sow also had a well-developed reproductive system, with two consecutive estrous periods. Its previous two litters consisted of 14 live piglets. A cloned female piglet was successfully cloned from the cultured piglet somatic cell population according to the nuclear transfer procedure of the cloned pig (Polejaeva et al., 2000). The cloned female pig and the donor female piglet were referred to as CP and DP, respectively. The CP and DP were both raised under the same environmental conditions and management. Blood is an important tissue which is associated with disease, immunity, and hormones. The methylation patterns of blood from the CP and DP were analyzed and found to be associated with litter size (Zhang et al., 2017a; Hwang et al., 2017b). Both blood and ear fibroblast samples from the CP and DP were collected for subsequent WGBS studies. Hereafter, blood and ear samples from the DP will be referred to as DP-blood and DP-ear, respectively. The same nomenclature will also be used for the CP. Both animals were raised in a fully enclosed barn with negative pressure ventilation, water curtain cooling, and a fully slatted floor. Both the DP and CP were fed alongside normal pregnant sows with a restricted feed regime and sufficient water during their pregnancy. All animal experiments were approved by the Scientific Ethic Committee of Huazhong Agricultural University (HZAUCA-2016-008).
DNA Extraction, Bisulfite Treatment, and cDNA Preparation
Ear vein blood (3 ml) and ear fibroblast (2 g) samples were collected from the 594-day-old DP and CP for DNA extraction and subsequent WGBS studies. Genomic DNA was extracted using the TIANamp genomic DNA kit (TIANamp Technologies, Beijing, China) and bisulfite conversion was carried out using the EZ DNA Methylation Direct Kit (Zymo Research Corporation, Orange County, CA, USA). DNA quality and quantity were determined by NanoDrop (NanoDrop Technologies, Wilmington, DE, USA).
WGBS Library Preparation and Data Quality Control
Four samples (DP-blood, DP-ear, CP-blood, and CP-ear) from the two pigs were used for WGBS sequencing. After a quality control test for the library preparation, these DNA samples were fragmented by Scientz18-A interrupter (SCIENTZ, Ningbo, China). The fragments were then end-repaired, plus A at the 3′ end, and ligated with adapters. Fragments of 400–500 base pairs (bp) were selected by agarose gel electrophoresis, then treated with bisulfite and subjected to PCR amplification to form the sequencing library. The sequencing library which passed a quality control test was sequenced using Illumina HiSeq Xten (Biomarker Technologies, Beijing, China). The sequencing reads were converted into FASTQ files by Illumina Casava1.8. To ensure the quality of the sequencing reads, we used FastQC software with the following criteria: (1) the removal of duplicate reads due to PCR amplification, (2) reads have adapters, (3) reads with more than 10% N content or more than 50% low quality bases. The raw sequencing data were stored in the National Centre for Biotechnology Information (NCBI) public database of Sequence Read Archive (SRA) (data number SUB5307644).
Mapping Reads to the Reference Genome
After bisulfite treatment and PCR amplification, the unmethylated cytosines (C) from the genome were converted to thymine (T), and the methylated C remained unchanged. All of the C T conversions in the coding strand and guanine adenine (G A) conversions in the noncoding strand of the genome were compared with the reference genome (Sus_scrofa11.1) using Bowtie2, and the unique matches were subsequently analyzed. The filtered clean reads that aligned with the same region of the genome were used to summarize the sequencing depth and coverage (Krueger and Andrews, 2011). The percentage of methylated clean reads compared to the total number of clean reads in the lambda genome was calculated as the conversion rate of bisulfite using the Bismark software. The Bismark_methylation_extractor was used to identify and calculate the distribution of different types (CG, CHG, and CHH) of methylated C reads (Krueger and Andrews, 2011). Considering the sample size, we mapped the clean reads to the reference genome four times and retained only the reads that mapped to the same position of the reference genome four times for subsequent bioinformatic analysis. A binomial distribution test was performed to identify 5mC for each C site, and then the P value of binomial distribution test was subjected to false discovery rate (FDR) correction using the Benjamini & Hochberg method in the p.adjust package of the R software. The C site with coverage > 4X and an FDR < 0.05 was considered to be a true methylation site (Zhang et al., 2017b; Furci et al., 2019).
Estimating 5mC Methylation Levels and Identifying Differential Methylation Regions (DMRs)
The methylation level of the C site with methylation frequencies are higher than the expected background in the bisulfite conversion reaction and sequencing errors in the binomial test detection were estimated. The methylation level was defined as the ratio between the number of reads of a methylated C to the total number of reads of that site (Schultz et al., 2012). We used the Kolmogorov–Smirnov test to determine the significance of the methylation differences between the CP and the DP, p < 0.05 was considered to be statistically significant. We also calculated the correlation of methylation levels by using the Pearson correlation coefficient of the samples. The regions with different methylation sites greater than 3, and in which the difference in methylation levels was no less than 0.2 (0.3 for CG type) with p value from Fisher's exact test of less than 0.05.) were used for DMR analysis with MOABS (Sun et al., 2014).
Gene Ontology and Pathway Analysis
In order to identify genes associated with DMRs and to understand their functions, we used a homemade perl script (Enrich_analysis.pl) to perform a Kyoto Encyclopedia of Genes and Genomes (KEGG) pathway enrichment analysis. We also used topGO to perform gene function enrichment analysis via Gene Ontology (GO). GO terms and KEGG pathways with corrected p-value <0.05 were defined as enriched clusters. The genes that were located within the DMRs or closest to the DMRs of the intergenic region were detected and defined as differential methylation genes (DMGs) by ChIPseeker on the R package. All the DMGs underwent functional analysis using GO and KEGG.
Results
Reproductive Phenotype of the CP and DP
Under the same feeding conditions, the litter size of first three parities of the CP was 8, 12, and 21. The litter size of the first three parities of the DP was 12, 20, and 14. The total litter size of the CP was 4 piglets less than that of the DP in the first two parities and one piglet less in the first three parities (Table 1).
DNA Methylation Mapping and Patterns
Genome bisulfite sequencing results are shown in Table 2. A total of 366.43 GB of clean sequence data was obtained for the four samples. Appropriately 91.61 GB of clean data was obtained for each sample, with an average sequencing depth of equal to or greater than 25X. Methylation conversion with more than 87.82% of reads have >10X coverage and the quality control analysis showed that Q30 of the sequence data were more than 90.72%. In the mapping step, 74.93% of reads were aligned to the reference genome. These results suggested that our sequencing data was of high quality, which is beneficial for subsequent bioinformatic analyses.
Three methylation patterns of CG, CHG, and CHH (H = A, T, or C) were obtained, and the methylated contexts and proportions of the four samples were similar. Methylated C in the blood sample of the DP was composed of 97.10% CG, 0.73% CHG, and 2.17% CHH compared to 96.68% CG, 0.83% CHG, and 2.49% CHH in the CP. Methylated C from ear tissue samples was composed of 96.37% CG, 0.86% CHG, and 2.77% CHH in the DP, and 95.92% CHG, 0.98% CHG, and 3.1% CHH in the CP (Figure 1). The majority of the methylation was CG, indicating that this type of methylation is dominant in this breed of pigs.
Figure 1 The average ratio of DNA methylation types in the genome of donor and cloned pigs. Panels (A–D) are the methylation patterns of DP-blood, CP-blood, DP-ear and CP-ear, respectively. The color of the pie chart indicates the type of methylation, and the number represents the proportion of the methylation type.
Methylation of Different C Base Types
The C bases (CG, CHG, and CHH) of different distribution types had different levels of methylation between different animals and tissues. We counted the distribution of C-base methylation levels for each type using a violin plot. CG methylations were the most abundant, and most of them were hypermethylated with a methylation level greater than 0.7 (Figure 2A). The CHG and CHH methylation levels were mostly below 0.2 (Figures 2B, C). The chromosomal methylation profile also showed that CG methylation occurred more frequently, and that mCHG and mCHH methylation levels were considerably lower than mCG methylation level, which is consistent with the results of the chromosome methylation map (Supplementary Figure 1). Therefore, the methylation status of blood or ear tissue samples and the number of mCs between DP and CP showed existential discrepancy. Next, we analyzed the relationship between sequence context and methylation preference, mapping the 9 bp sequence characteristics around the mCs (Figure 3). We observed that the mCG was completely hypermethylated, and that CAG and CAT were the most common motifs of CHG and CHH in the two pigs.
Figure 2 Violin plot for the overall distribution of methylation levels for different methylation types. Panels (A–C) are the context of CG, CHG, and CHH. H = A, C, or T. DP-blood, CP-blood, DP-ear, and CP-ear were the sample name. The abscissa represents the different samples, the ordinate represents the level of methylation of the samples; the points represented the methylation levels and cross-sectional area indicated the density of the methylation level; the boxplot shows the methylation levels in each violin.
Figure 3 Methylation features of 9 bp spanning 5mC methylcytosine sites. CG, CHG, CHH is the methylation type, H is A, C, or T. Panels (A–D) are the 4 bp sequences of the upstream and downstream of the methylated cytosine of DP-blood, CP-blood, DP-ear, and CP-ear, respectively. The abscissa denotes the base number of the methylation site, the total height of each position is the sequence conservation of the base, which represents the relative frequency of the base at that position.
DNA Methylation of CpG Islands in Functional Regions
We divided all mCs into functional 5′ untranslated region (UTR), first exon, inner exon, first intron, inner intron, last exon, and 3′ UTR. We also measured the methylation of CpG islands in each functional region of the four samples (Figure 4). The CpG islands of CG-type methylation in the functional regions of tissues from both pigs were similar. Furthermore, the levels of methylation in the intron regions, exon regions (except the first exon), and downstream were significantly higher than in other regions. CHH-type CpG islands were hypomethylated and stable in each functional region, the level of methylation in the CP was higher than that of the DP, and CHG type CpG islands in the four samples were almost unmethylated.
Figure 4 CpG islands distribution in different functional elements. Panels (A–D) are CpG islands distribution in the functional elements of DP-blood, CP-blood, DP-ear, and CP-ear, respectively. The methylation region was divided into upstream, general gene body, downstream, TSS and TTS regions. Each region was divided into 20 bins to count weight methylation levels. The abscissa represents the gene functional elements, line color represent methylation type. The left ordinate represents the methylation levels of CG and CHG, and the right ordinate represents the methylation levels of CHH.
Methylation Annotation of the Hypermethylated CpG Islands
We took the 3,000 bp upstream of a gene as promoter region, made an overlap annotation on CpG islands with methylation levels >0.7 and mC coverage >5X, but not including the hypermethylated CpG island with C-degree confidence less than 0.1 (Figure 5). In the blood, hypermethylated CpG islands were found in the distal intergenic regions and other exons (except the first exon), but the frequency of methylation were different, while other hypermethylated CpG islands of the two pigs were found in different gene functional regions. In the DP, 73.08% of hypermethylated CpG islands were found at distal intergenic regions, 7.69% of hypermethylated CpG islands were found at other exon (exon region except first exon), 3.85% of hypermethylated CpG islands occurred in the 5′UTR, and 15.38% of hypermethylated CpG islands were found at other introns (intron region except first intron). In contrast, 76.92% of CpG islands were found in distal intergenic regions, 15.38% of hypermethylated CpG islands were found at other exon, 3.85% of hypermethylated CpG islands occurred in a promoter region, and 3.85% of hypermethylated CpG islands were found at first exon regions in the CP. In the ear tissue, hypermethylated CpG islands of the two pigs occurred in the 5′UTR, distal intergenic, intron, and other exon, but their rate of hypermethylation was different. In the DP, 11.11% of hypermethylated CpG islands were found in other exons, 22.22% of hypermethylated CpG islands were found in other introns, 5.56% of hypermethylated CpG islands were found in the first intron, and 61.11% of hypermethylated CpG islands were found in the distal intergenic region. In the CP, these figures were 13.64%, 9.09%, 4.55%, and 68.18%, respectively. Hypermethylated CpG islands were found at distal intergenic region in the DP compared to the CP. There were also some differences in the functional regions and proportions of hypermethylated CpG islands in the two pigs, in which hypermethylated CpG islands in the blood were simultaneously present in the intergenic region and other exons, while the functional regions of hypermethylated CpG islands in the ear tissue were found to be the same just at different proportions.
Figure 5 The ratio of hypermethylation CpG islands distribution in different functional elements. Panels (A, B) are the CpG islands annotation in functional areas of DP-blood and CP-blood, respectively. Panels (C, D) are the CpG islands annotation in functional areas of DP-ear and CP-ear, respectively. The color of the pie chart indicates the functional area and number in parentheses represents the percentage of methylation in the functional area.
Correlation Analysis of DP and CP
Correlation analysis between samples with different types of methylation (Figure 6) showed that the correlation coefficient of the CG type of the four samples ranged from 0.63 to 0.80 while CHH type ranged from 0.59 to 0.70, and CHG type ranged from 0.43 to 0.55. Principal component analysis also showed that the type of methylation in each functional region of the four samples was different. Furthermore, CG methylation in the blood samples of the two pigs showed some differences, while CHG methylation between DP-blood sample and CP-ear tissue samples also exhibited some differences. The same was true for CHH methylation of each functional element in DP-blood and CP-ear samples (Figure 4).
Figure 6 The correlation analysis of samples. The histogram of methylation level in the figure, the abscissa is the methylation level, and the ordinate is proportion of the corresponding methylation level; the value in the figure is the correlation coefficient of two samples. Panels (A–C) represented the ratio of DMRs in different gene functional regions of CG, CHG, and CHH of the blood, respectively.
Analysis of DMRs
We annotated the DMGs with gene functional elements (with 3,000 bp 5′ to the TSS and 3′ to transcription termination site as the gene upstream and downstream functional regions, respectively) that associate with DMRs according to the DMR position and genome annotation information (Zhang et al., 2017a). DMRs was annotated according to different methylation types (Figure 7). Finally, a total of 210 CG-type, 1 CHG-type, and 7 CHH-type DMRs were identified in all blood samples. Additionally, 764 CG-type, 2 CHG-type, and 9 CHH-type DMRs were identified in the ear tissue samples. Most of these DMRs were located in the gene coding regions and only a few were found in the intergenic regions, in first exon, and in the UTR. The majority of CG type DMRs were found in intergenic regions, second only to intron regions. CHG-type DMRs were mainly located in exons, while CHH-type DMRs were found in the distal intergenic and exon regions. In summary, a total of 215 DMGs were identified in the blood (Supplementary File 1) and 707 DMGs in the ear (Supplementary File 2).
Figure 7 The ratio of DMRs in gene functional regions. Panels (A–C) represented the ratio of DMRs in different gene functional regions of CG, CHG, and CHH of the blood, respectively. Panels (D–F) represented the ratio of DMRs in different gene functional regions of CG, CHG, and CHH of the ear, respectively.
The genes that were located within or closest to the DMRs of the intergenic region were detected and defined as DMGs by ChIPseeker in R. All of the DMGs are enriched in biological processes, cell components, and molecular functions. In the blood, DMGs were enriched in all biological processes except those involved in the rhythmic process, hormone secretion, and cell killing. Cellular components except for the virion, virion part and collagen trimer, molecular functions of antioxidants, receptor regular, and translation regular activity, chemoattractant, and protein tag, in which many DMGs were enriched in immune-related and female reproduction pathways (Figure 8C, D). In the ear, DMGs were enriched in all the GO terms except cell killing, molecular function of chemointermediary activity, and translation regular activity. Moreover, many DMGs are enriched on immune-related and reproduction terms (Figure 8A, B). KEGG analysis showed that DMGs are mainly enriched on ontogeny, metabolism, and tumor disease, in which many DMGs are enriched in immunity and reproduction associated pathways.
Figure 8 Top GO and top KEGG pathway analysis of CG type DMGs. Panels (A, C) were top GO analysis and top KEGG analysis in blood. Panels (B, D) were top GO analysis and top KEGG analysis in ear, respectively.
DMGs Enriched in GO and KEGG Pathways Related to Immunity
To determine whether the cloning process will lead to immune resistance variations between the two pigs, we performed functional enrichment analysis and annotation via GO and KEGG for DMGs. Some DMGs are enriched in T/B cell differentiation, thymus development, and other immune-related terms/pathways. In the blood, 16 DMGs are enriched in immunologically related GO terms (Supplementary File 3), in which CARD11, CD244, RFTN1, and THBS1 are enriched in three or more immune-related GO terms. Interestingly, a novel gene (gene ID: ENSSSCG00000016737) is enriched in antibiotic response processes. KEGG analysis identified eight DMGs which are enriched in immune-related pathways. PTPRC, GRB2, and CARD11 are enriched in three or more immune-related pathways (Supplementary File 4), in which a new gene (gene ID: ENSSSCG00000035088) associated with B and T cell receptor activity and cytotoxicity of natural killer cells was identified following KEGG analysis. Moreover, CARD11 and CD244 are enriched in immune-related terms/pathways. In ear tissue samples, 51 annotated genes and 7 new genes are enriched for the terms “immunity” and “differential methylated” (Supplementary File 5). Among these, IRAK3, IL17RA, FGR, PDE4B, CD244, TGFB2, MX1, F2RL1, MEF2C, GLI3, CD83, FOXP1, CDK6, LEP, ENSSSCG00000037390, and ENSSSCG00000037643 were enriched in three or more immune-related terms. Furthermore, 22 DMGs were enriched in the immunologically related pathways (Supplementary File 6) with ENSSSCG00000039951 and PTPRC being enriched in three or more immune-related pathways. NCK2, PLCG1, TANK, TGFB2, and ENSSSCG00000006851 are simultaneously enriched on immune-related terms and pathways.
DMGs Enriched in GO and KEGG Pathways Related to Reproductive Performance
A total of 55 and 20 DMGs are enriched on oocyte maturation, embryo development, gonad development, estrogen synthesis, and estrogen activity. Additionally, 11 DMGs in blood (Supplementary File 7) and 40 DMGs in ear samples (Supplementary File 8) are enriched in reproductive related biological processes, whereas seven DMGs in blood (Supplementary File 9) and 13 DMGs in ear (Supplementary File 10) are enriched in reproductive related pathways. GO analysis has shown that most DMGs are enriched in processed involving embryonic development, in which NR2F1, TUBGCP5, RUNX2, and ASH2L are related to steroid hormone signaling, steroid receptor activity, oocyte meiosis cycle, female gonad development, or estrogen responses. In the KEGG analysis, most of the DMGs are mainly enriched in the signaling pathways of oxytocin, estrogen, GnRH, prolactin, and oocyte meiosis; among which, GANAS is located in the reproductive-hormone-related pathways and is associated with follicular development. Furthermore, in blood GNAS is enriched in post-embryonic body morphogenesis and embryonic hindlimb morphogenesis terms, as well as in ovarian steroidogenesis, oxytocin, estrogen, and GnRH signaling pathways (Table 3). EGFR, another DMG in ear samples is enriched to embryonic placental development, estrogen, GnRH, and Oxytocin signaling pathways (Table 3).
Discussion
Phenotype and Methylation Were Different Between the CP and DP
Under the same feeding regime and environmental conditions, the first calving age of the CP was 51 days older than the DP, and the litter size of the CP was four piglets less and two live piglets less than the DP in the first two parities and in the first three parities, respectively (Table 3). DNA methylation was a common epigenetic modification and affected oocyte meiotic maturation and embryonic development, which in turn affecting the litter size (Gao et al., 2016; Glanzner et al., 2017; Hwang et al., 2017a; Wang et al., 2018). In this study, we found that the level of methylation and the methylated functional regions of the same tissue between the two pigs were different. Furthermore, in both pigs the majority of methylation in the blood occurs in the intergenic region and other exon regions, while other methylations occurred in different regions. In the ear, hypermethylated CpG islands in both animals were found in the same functional region, but the hypermethylation rate of each functional region was different. 3.85% of hypermethylated CpG islands occurred in the promoter regions in the blood of cloned pig, and the proportion of hypermethylated CpG islands occurred in exons was significantly higher than that of donor pigs in blood and ear (with 19.23% vs 7.69% in blood, 13.64% vs 11.11% in ear). This suggested that DNA methylation status in the promoter and exon regions of genes altered in cloned pig, indicating a potential risk of abnormal gene expression thus leading to changes in the cloned pig's phenotype. This was consistent with the results of the DNA methylation of promoter regions were changed in some genes and the DNA methylation status with a trend towards slightly increased methylation levels in cloned pigs (Gao et al., 2011). DNA methylation levels are negatively and monotonically correlated with gene expression levels around the TSS of genes (Zou et al., 2016), this change may lead a poor reproduction performance in the cloned pig. In our study, the DMRs (CP vs DP) suffered more hypomethylation (391/775 DMRs) than hypermethylation (384/775 DMRs) in the ear, suggested the DNA methylation status in the DMRs altered and phenotype changed in cloned pig, this was consistent with the DNA methylation context of skeletal muscle in abnormal and normal cloned pig (Li et al., 2014). Therefore, the methylation context of the ear and blood were different, which was consistent with the results of the DNA methylation in the ovary and intestine of pigs (Yuan et al., 2016). Furthermore, the average sequencing depth was more than 25X, the methylation conversion of reads with a coverage of >10X were more than 87.82% and >74.93% of bases were aligned to the reference genome. In our study, the majority of methylation sites were the CG type and most of them were hypermethylated, this was consistent with the maximum number of CG-type DMRs and DMGs. This shows that we have high quality data allowing us to be confident that our bioinformatic analyses are reliable.
DMGs Involved in Immunity and Reproduction
This study has identified a number of immune-related regulatory factors such as CD244, CD2, CD83, CD5, and CDK6 with differences in methylation between the CP and the DP. CD244 is known as an important regulator of various autoimmune defects, where low expression of CD244 on monocytes can cause systemic lupus erythematosus and the overexpression of CD244 could accelerate the cancer pathogenesis (Zhang et al., 2017c; Mak et al., 2018). Inhibition of CD244 expression may promote the expression of IFN-γ/TNF-α in CD8+ T cells and may reduce the incidence and mortality of bovine tuberculosis (Yang et al., 2013; Wang et al., 2015). Binding of CD244 to CD2 can reduce the production of antigen-specific interleukin-2 and inhibit the immune response of antiviral CD8+ T cells (Clarkson and Brown, 2009; Pacheco et al., 2013; McArdel et al., 2016). CD83, a specific marker of dendritic cells (DCs) for antigen presentation and lymphocyte activation, can regulate immune responses to prevent colitis by regulating DC surface proteins. CD83 can also promote the expression of some immune-related genes or the secretion of immune factors (Tze et al., 2011). Found in T lymphocytes, CD5 is a type of cluster differentiation antigen. High CD5 expression in peripheral T cells allows them to better respond to foreign antigens and enrich memory populations to produce stronger antigenic responses when exposed to the same antigen again (Freitas et al., 2017). Differences in CD5 methylation may lead to a change in the expression of CD5 in the peripheral blood of the CP and DP in this study, which may have had an effect on the immune responses of the two pigs during immunization programs. CD5 antigen-positive T lymphocytes play a role in signal transduction. These cells work with some immune-related genes and immune factors to maintain immune homeostasis and normal immune function of the thymus (Freitas et al., 2018). Cyclin-dependent kinase, CDK6, is a driving factor for cell proliferation and can bind D cell cycle protein 3 to promote tumor cell apoptosis (Wang et al., 2017). A lack of CDK6 usually leads to anemia, while overexpression leads to hematologic malignancies (Kollmann et al., 2016; Uras et al., 2017). Most of the aforementioned DMGs were enriched in cancer, T/B, NK cell differentiation, and other immune pathways, which are important for the prevention and development of disease (Kaur et al., 2013; Kuwabara et al., 2018; Mazzone et al., 2019). The differences in methylation of these genes may have an effect on immune responses in the CP. However, as the sampled pig or its offspring showed no significant difference in disease characteristics, and the sample size in our study was limited, further investigation is required to confirm this hypothesis.
To explain the difference between the first birth age and total number born (TNB) between the DP and the CP, we annotated the DMGs involved in reproduction. ITPR2, CALML5, TAPT1, GNAS, and EGFR were enriched in three or more reproduction related pathways (Table 3). ITPR2 is an important regulator of calcium-channel activity, which is the basis of animal fertilization and embryonic development (Kuroda et al., 1999; Miao and Williams, 2012). TAPT1 defects can lead to fatal chondrodysplasia of the fetus (Symoens et al., 2015) and CALML5 regulates epidermal differentiation by binding to stratifin (Sun et al., 2015). Differences in the expression of these genes may explain the delayed first calving age and the fewer offspring of the CP. GNAS has been shown to have a direct impact on follicular development. Moreover, inhibition of GNAS expression improved the success rate of embryonic implantation (Meng et al., 2014). Abnormal methylation of GNAS has also been shown to lead to neural tube defects (NTDs) during pregnancy, and even to embryonic development failure (Wang et al., 2017). Differences in GNAS expression may explain the differences in the number of litters between the DP and the CP in our study. EGFR is involved in the signaling pathways of mature oocytes and embryonic development, and expression of EGFR can inhibit oocyte development and maturation by stimulating PI3K/AKT signaling (Park et al., 2017). Inhibition of EGFR activates extracellular signal-regulated kinases 1 and 2 (ERK1/2) and accelerates oocyte meiosis (Lei et al., 2018). Conversely, overexpression of EGFR promotes progesterone and estrogen secretion thereby obstructing embryonic implantation and development (Wang et al., 2018). Differences in the methylation of EGFR may be an important factor in the delayed first calving age and the smaller litter size in the CP. In addition to ITPR2, CALML5, TAPT1, and GNAS, many DMGs were enriched in functions related to reproduction or pathways involving oxytocin, prolactin, GnRH, estrogen, ovarian hormones, follicular development, luteinizing hormon (LH) production, gonadal development, or reproductive activities of the females. Consequently, these DMGs may have an impact on the reproductive performance of sows (Napso et al., 2018).
Conclusion
In this study, we investigated the whole-genome methylation patterns in blood and ear tissue of a CP and a DP. Using this approach, we were able to identify DMGs in both sample types. A total of 60 DMGs were enriched in immune pathways, of which 19 were enriched in more than three immunity pathways. A total of 68 DMGs were enriched in reproductive related pathways, of which 13 participated in more than three pathways related to reproduction. Among which ITPR2, TAPT1, GNAS, and EGFR were related to reproduction. Our research provides basic genomic reference for clonal epigenetic research.
Data Availability Statement
The raw data for this study can be found in the NCBI public database at this URL link: https://www.ncbi.nlm.nih.gov/Traces/study/?acc=PRJNA530044.
Ethics Statement
This study was carried out in accordance with the recommendations of ' Scientific Ethic Committee of Huazhong Agricultural University (HZAUCA-2016-008), Wuhan, China '. The protocol was approved by the ' Scientific Ethic Committee of Huazhong Agricultural University(HZAUCA-2016-008)'.
Author Contributions
XL, BZ, YM, and SZ designed the experiments. XL, BZ, YM, GM, SF, and MW carried out the experiments. MW, XD, JR, and HW conducted the statistical analysis and discussion. XL, XD, SZ, JR, and MW organized and wrote the manuscript. All the authors have read the paper and have agreed to be co-authors.
Funding
This work was supported by the Fundamental Research Funds for the Central Universities of China (2662016PY006), National Natural Science Foundation of China (31572375), the Fundamental Research Funds for the Central Universities of China (2662018JC033), the Fundamental Research Funds for the Central Universities (2662017JC027), the National Swine Industry Technology System (No.CARS-35), and Da Bei Nong Group Promoted Project for Young Scholar of HZAU (2017DBN019).
Conflict of Interest
Authors MW, GM, XL and XD were employed by company Guangxi Yangxiang Co., Ltd, China.
The remaining authors declare that the research was conducted in the absence of any commercial or financial relationships that could be construed as a potential conflict of interest.
Acknowledgement
We are grateful to Guangxi Yangxiang Co., Ltd. For providing the samples of the high-yield lean-type sows for the DP. We would like to thank Yu Zhou and Xiangdong Xu for their help during the sampling process. We are also grateful to Dr. Zhiquan Wang and Dr. Hafiz Ishfaq Ahmad for their suggestions in the manuscript writing. Therefore, we would like to thank Editage (www.editage.cn) for English language editing.
Supplementary Material
The Supplementary Material for this article can be found online at: https://www.frontiersin.org/articles/10.3389/fgene.2020.00023/full#supplementary-material
Abbreviations
CP, cloned pig; DMGs, differential methylation genes; DMRs, differential methylation regions; DP, donor pig; SCNT, somatic cell nuclear transplantation.
References
Adachi, N., Yamaguchi, D., Watanabe, A., Miura, N., Sunaga, S., Oishi, H., et al. (2014). Growth, reproductive performance, carcass characteristics and meat quality in F1 and F2 progenies of somatic cell-cloned pigs. J. Reprod. Dev. 60 (2), 100–105. doi: 10.1262/jrd.2012-167
Akagi, S., Matsukawa, K., Takahashi, S. (2014). Factors affecting the development of somatic cell nuclear transfer embryos in cattle. J. Reprod. Dev. 60 (5), 329–335. doi: 10.1262/jrd.2014-057
Auclair, G., Weber, M. (2012). Mechanisms of DNA methylation and demethylation in mammals. Biochimie 94 (11), 2202–2211. doi: 10.1016/j.biochi.2012.05.016
Bogdanovic, O., Gomez-Skarmeta, J. L. (2014). Embryonic DNA methylation: insights from the genomics era. Brief Funct. Genomics 13 (2), 121–130. doi: 10.1093/bfgp/elt039
Clarkson, N. G., Brown, M. H. (2009). Inhibition and activation by CD244 depends on CD2 and phospholipase C-gamma1. J. Biol. Chem. 284 (37), 24725–24734. doi: 10.1074/jbc.M109.028209
Freitas, C. M. T., Hamblin, G. J., Raymond, C. M., Weber, K. S. (2017). Naïve helper T cells with high CD5 expression have increased calcium signaling. PloS One 12 (5), e0178799–e0178799. doi: 10.1371/journal.pone.0178799
Freitas, C. M. T., Johnson, D. K., Weber, K. S. (2018). T cell calcium signaling regulation by the co-receptor CD5. Int. J. Mol. Sci. 19 (5), 1295. doi: 10.3390/ijms19051295
Furci, L., Jain, R., Stassen, J., Berkowitz, O., Whelan, J., Roquis, D., et al. (2019). Identification and characterisation of hypomethylated DNA loci controlling quantitative resistance in Arabidopsis. eLife 8, e40655. doi: 10.7554/eLife.40655
Gao, F., Luo, Y., Li, S., Li, J., Lin, L., Nielsen, A. L., et al. (2011). Comparison of gene expression and genome-wide DNA methylation profiling between phenotypically normal cloned pigs and conventionally bred controls. PloS One 6 (10), e25901. doi: 10.1371/journal.pone.0025901
Gao, Y. Y., Chen, L., Wang, T., Nie, Z. W., Zhang, X., Miao, Y. L. (2016). Oocyte aging-induced Neuronatin (NNAT) hypermethylation affects oocyte quality by impairing glucose transport in porcine. Sci. Rep. 6, 36008. doi: 10.1038/srep36008
Glanzner, W. G., Wachter, A., Coutinho, A. R., Albornoz, M. S., Duggavathi, R., GonÇAlves, P. B., et al. (2017). Altered expression of BRG1 and histone demethylases, and aberrant H3K4 methylation in less developmentally competent embryos at the time of embryonic genome activation. Mol. Reprod. Dev. 84 (1), 19–29. doi: 10.1002/mrd.22762
Hannon, E., Knox, O., Sugden, K., Burrage, J., Wong, C. C. Y., Belsky, D. W., et al. (2018). Characterizing genetic and environmental influences on variable DNA methylation using monozygotic and dizygotic twins. PloS Genet. 14 (8), e1007544–e1007544. doi: 10.1371/journal.pgen.1007544
Hwang, J. H., An, S. M., Kwon, S., Park, D. H., Kim, T. W., Kang, D. G., et al. (2017). DNA methylation patterns and gene expression associated with litter size in Berkshire pig placenta. PloS One 12 (9), e0184539–e0184539. doi: 10.1371/journal.pone.0184539
Jones, P. A., Takai, D. (2001). The role of DNA methylation in mammalian epigenetics. Science 293 (5532), 1068. doi: 10.1126/science.1063852
Kaur, G., Trowsdale, J., Fugger, L. (2013). Natural killer cells and their receptors in multiple sclerosis. Brain: A J. Neurol. 136 (Pt 9), 2657–2676. doi: 10.1093/brain/aws159
Keefer, C. L. (2015). Artificial cloning of domestic animals. Proc. Natl. Acad. Sci. U. States America 112 (29), 8874–8878. doi: 10.1073/pnas.1501718112
Kollmann, K., Heller, G., Schneckenleithner, C., Warsch, W., Scheicher, R., Ott, R. G., et al. (2016). A kinase-independent function of CDK6 links the cell cycle to tumor angiogenesis. Cancer Cell 30 (2), 359–360. doi: 10.1016/j.ccell.2016.07.003
Krueger, F., Andrews, S. R. (2011). Bismark: a flexible aligner and methylation caller for bisulfite-seq applications. Bioinf. (Oxf. Engl.) 27 (11), 1571–1572. doi: 10.1093/bioinformatics/btr167
Kuroda, Y., Kaneko, S., Yoshimura, Y., Nozawa, S., Mikoshiba, K. (1999). Are there inositol 1,4,5-triphosphate (IP3) receptors in human sperm? Life Sci. 65 (2), 135–143. doi: 10.1016/S0024-3205(99)00230-1
Kuwabara, T., Matsui, Y., Ishikawa, F., Kondo, M. (2018). Regulation of t-cell signaling by post-translational modifications in autoimmune disease. Int. J. Mol. Sci. 19 (3), 819. doi: 10.3390/ijms19030819
Lee, S. H., Oh, H. J., Kim, M. J., Kim, G., Setyawan, E., Ra, K., et al. (2018). Dog cloning—no longer science fiction. Reprod. Dev. doi: 10.1111/rda.13358
Lei, R., Bai, X., Chang, Y., Li, J., Qin, Y., Chen, K., et al. (2018). Effects of Fullerenol nanoparticles on rat oocyte meiosis resumption. Int. J. Mol. Sci. 19 (3), 699–699. doi: 10.3390/ijms19030699
Li, E., Zhang, Y. (2014). DNA methylation in mammals. Cold Spring Harbor Perspect. in Biol. 6 (5), a019133–a019133. doi: 10.1101/cshperspect.a019133
Li, Z., Sun, X., Chen, J., Liu, X., Wisely, S. M., Zhou, Q., et al. (2006). Cloned ferrets produced by somatic cell nuclear transfer. Dev. Biol. 293 (2), 439–448. doi: 10.1016/j.ydbio.2006.02.016
Li, G., Jia, Q., Zhao, J., Li, X., Yu, M., Samuel, M. S., et al. (2014). Dysregulation of genome-wide gene expression and DNA methylation in abnormal cloned piglets. BMC Genomics 15, 811. doi: 10.1186/1471-2164-15-811
Mak, A., Thornhill, S. I., Lee, H. Y., Lee, B., Poidinger, M., Connolly, J. E., et al. (2018). Brief report: decreased expression of CD244 (SLAMF4) on monocytes and platelets in patients with systemic lupus erythematosus. Clin. Rheumatol. 37, 3, 811–816. doi: 10.1007/s10067-017-3698-2
Mazzone, R., Zwergel, C., Artico, M., Taurone, S., Ralli, M., Greco, A., et al. (2019). The emerging role of epigenetics in human autoimmune disorders. Clin. Epigenet. 11 (1), 34–34. doi: 10.1186/s13148-019-0632-2
McArdel, S. L., Terhorst, C., Sharpe, A. H. (2016). Roles of CD48 in regulating immunity and tolerance. Clin. Immunol. (Orlando Fla.) 164, 10–20. doi: 10.1016/j.clim.2016.01.008
Meng, Y., Guo, Y., Qian, Y., Guo, X., Gao, L., Sha, J., et al. (2014). Effects of GnRH antagonist on endometrial protein profiles in the window of implantation. Proteomics 14 (20), 2350–2359. doi: 10.1002/pmic.201400145
Miao, Y. L., Williams, C. J. (2012). Calcium signaling in mammalian egg activation and embryo development: the influence of subcellular localization. Mol. Reprod. Dev. 79 (11), 742–756. doi: 10.1002/mrd.22078
Napso, T., Yong, H. E. J., Lopez-Tello, J., Sferruzzi-Perri, A. N. (2018). The role of placental hormones in mediating maternal adaptations to support pregnancy and lactation. Front. in Physiol. 9, 1091–1091. doi: 10.3389/fphys.2018.01091
Ogura, A., Inoue, K., Ogonuki, N., Lee, J., Kohda, T., Ishino, F. (2002). Phenotypic effects of somatic cell cloning in the mouse. Cloning Stem Cells 4 (4), 397–405. doi: 10.1089/153623002321025078
Ogura, A., Inoue, K., Wakayama, T. (2013). Recent advancements in cloning by somatic cell nuclear transfer. Philos. Trans. R. Soc. Lond. Ser. B Biol. Sci. 368 (1609), 20110329–20110329. doi: 10.1098/rstb.2011.0329
Pacheco, Y., McLean, A. P., Rohrbach, J., Porichis, F., Kaufmann, D. E., Kavanagh, D. G. (2013). Simultaneous TCR and CD244 signals induce dynamic downmodulation of CD244 on human antiviral T cells. J. Immunol. 191 (5), 2072–2081. doi: 10.4049/jimmunol.1300435
Park, H. J., Chae, S. K., Kim, J. W., Yang, S. G., Jung, J. M., Kim, M. J., et al. (2017). Ganglioside GM3 induces cumulus cell apoptosis through inhibition of epidermal growth factor receptor-mediated PI3K/AKT signaling pathways during in vitro maturation of pig oocytes. Mol. Reprod. Dev. 84 (8), 702–711. doi: 10.1002/mrd.22848
Polejaeva, I. A., Chen, S.-H., Vaught, T. D., Page, R. L., Mullins, J., Ball, S., et al. (2000). Cloned pigs produced by nuclear transfer from adult somatic cells. Nature 407 (6800), 86–90. doi: 10.1038/35024082
Rideout, W. M., Eggan, K., Jaenisch, R. (2001). Nuclear cloning and epigenetic reprogramming of the genome. Science 293 (5532), 1093. doi: 10.1126/science.1063206
Schultz, M. D., Schmitz, R. J., Ecker, J. R. (2012). ‘Leveling' the playing field for analyses of single-base resolution DNA methylomes. Trends In Genet.: TIG 28 (12), 583–585. doi: 10.1016/j.tig.2012.10.012
Shen, C.-J., Cheng, W. T. K., Wu, S.-C., Chen, H.-L., Tsai, T.-C., Yang, S.-H., et al. (2012). Differential differences in methylation status of putative imprinted genes among cloned swine genomes. PloS One 7 (2), e32812–e32812. doi: 10.1371/journal.pone.0032812
Shibata, M., Otake, M., Tsuchiya, S., Chikyu, M., Horiuchi, A., Kawarasaki, T. (2006). Reproductive and growth performance in Jin Hua pigs cloned from somatic cell nuclei and the meat quality of their offspring. J. Reprod. Dev. 52 (5), 583–590. doi: 10.1262/jrd.18004
Shimozawa, N., Sotomaru, Y., Eguchi, N., Suzuki, S., Hioki, K., Usui, T., et al. (2006). Phenotypic abnormalities observed in aged cloned mice from embryonic stem cells after long-term maintenance. Reproduction 132 (3), 435–441. doi: 10.1530/rep.1.00745
Sun, D., Xi, Y., Rodriguez, B., Park, H. J., Tong, P., Meong, M., et al. (2014). MOABS: model based analysis of bisulfite sequencing data. Genome Biol. 15 (2), R38–R38. doi: 10.1186/gb-2014-15-2-r38
Sun, B. K., Boxer, L. D., Ransohoff, J. D., Siprashvili, Z., Qu, K., Lopez-Pajares, V., et al. (2015). CALML5 is a ZNF750- and TINCR-induced protein that binds stratifin to regulate epidermal differentiation. Genes Dev. 29 (21), 2225–2230. doi: 10.1101/gad.267708.115
Symoens, S., Barnes, A. M., Gistelinck, C., Malfait, F., Guillemyn, B., Steyaert, W., et al. (2015). Genetic defects in TAPT1 disrupt ciliogenesis and cause a complex lethal osteochondrodysplasia. Am. J. Hum. Genet. 97 (4), 521–534. doi: 10.1016/j.ajhg.2015.08.009
Tze, L. E., Horikawa, K., Domaschenz, H., Howard, D. R., Roots, C. M., Rigby, R. J., et al. (2011). CD83 increases MHC II and CD86 on dendritic cells by opposing IL-10-driven MARCH1-mediated ubiquitination and degradation. J. Exp. Med. 208 (1), 149–165. doi: 10.1084/jem.20092203
Uras, I. Z., Scheicher, R. M., Kollmann, K., Glosmann, M., Prchal-Murphy, M., Tigan, A. S., et al. (2017). Cdk6 contributes to cytoskeletal stability in erythroid cells. Haematologica 102 (6), 995–1005. doi: 10.3324/haematol.2016.159947
Wang, Y., Zhong, H., Xie, X., Chen, C. Y., Huang, D., Shen, L., et al. (2015). Long noncoding RNA derived from CD244 signaling epigenetically controls CD8+ T-cell immune responses in tuberculosis infection. Proc. Natl. Acad. Sci. U. S. A. 112 (29), E3883–E3892. doi: 10.1073/pnas.1501662112
Wang, H., Nicolay, B. N., Chick, J. M., Gao, X., Geng, Y., Ren, H., et al. (2017). The metabolic function of cyclin D3-CDK6 kinase in cancer cell survival. Nature 546 (7658), 426–430. doi: 10.1038/nature22797
Wang, L., Chang, S., Wang, Z., Wang, S., Huo, J., Ding, G., et al. (2017). Altered GNAS imprinting due to folic acid deficiency contributes to poor embryo development and may lead to neural tube defects. Oncotarget 8 (67), 110797–110810. doi: 10.18632/oncotarget.22731
Wang, W., Bai, G., Zhang, Y., Zhang, T., Li, C., Yuan, Y., et al. (2018). HBxAg suppresses cell apoptosis and promotes the secretion of placental hormones in human placental trophoblasts via activation of the EGFR/Akt pathway. Cell Biol. Int. 42 (2), 237–247. doi: 10.1002/cbin.10891
Wang, Y. K., Yu, X. X., Liu, Y. H., Li, X., Liu, X. M., Wang, P. C., et al. (2018). Reduced nucleic acid methylation impairs meiotic maturation and developmental potency of pig oocytes. Theriogenology 121, 160–167. doi: 10.1016/j.theriogenology.2018.08.009
Yang, B., Wang, X., Jiang, J., Cheng, X. (2013). Involvement of CD244 in regulating CD4+ T cell immunity in patients with active tuberculosis. PloS One 8 (4), e63261–e63261. doi: 10.1371/journal.pone.0063261
Yuan, X. L., Gao, N., Xing, Y., Zhang, H. B., Zhang, A. L., Liu, J., et al. (2016). Profiling the genome-wide DNA methylation pattern of porcine ovaries using reduced representation bisulfite sequencing. Sci. Rep. 6, 22138. doi: 10.1038/srep22138
Zhang, F., Liu, X., Chen, C., Zhu, J., Yu, Z., Xie, J., et al. (2017a). CD244 maintains the proliferation ability of leukemia initiating cells through SHP-2/p27(kip1) signaling. Haematologica 102 (4), 707–718. doi: 10.3324/haematol.2016.151555
Zhang, Y., Li, F., Feng, X., Yang, H., Zhu, A., Pang, J., et al. (2017b). Genome-wide analysis of DNA methylation profiles on sheep ovaries associated with prolificacy using whole-genome Bisulfite sequencing. BMC Genomics 18 (1), 759. doi: 10.1186/s12864-017-4068-9
Zhang, Y., Wilson, R., Heiss, J., Breitling, L. P., Saum, K.-U., Schöttker, B., et al. (2017c). DNA methylation signatures in peripheral blood strongly predict all-cause mortality. Nat. Commun. 8, 14617–14617. doi: 10.1038/ncomms14617
Keywords: DNA methylation, gene, immunity, methylation differential, pig, reproduction, somatic cloning
Citation: Wang M, Feng S, Ma G, Miao Y, Zuo B, Ruan J, Zhao S, Wang H, Du X and Liu X (2020) Whole-Genome Methylation Analysis Reveals Epigenetic Variation in Cloned and Donor Pigs. Front. Genet. 11:23. doi: 10.3389/fgene.2020.00023
Received: 05 March 2019; Accepted: 08 January 2020;
Published: 20 February 2020.
Edited by:
Kazuhiko Nakabayashi, National Center for Child Health and Development (NCCHD), JapanReviewed by:
Zhe Zhang, South China Agricultural University, ChinaEveline M. Ibeagha-Awemu, Agriculture and Agri-Food Canada (AAFC), Canada
Copyright © 2020 Wang, Feng, Ma, Miao, Zuo, Ruan, Zhao, Wang, Du and Liu. This is an open-access article distributed under the terms of the Creative Commons Attribution License (CC BY). The use, distribution or reproduction in other forums is permitted, provided the original author(s) and the copyright owner(s) are credited and that the original publication in this journal is cited, in accordance with accepted academic practice. No use, distribution or reproduction is permitted which does not comply with these terms.
*Correspondence: Xiaoyong Du, duxiaoyong@mail.hzau.edu.cn; Xiangdong Liu, liuxiangdong@mail.hzau.edu.cn