ORIGINAL RESEARCH article
Transcriptome Analysis Identifies Candidate Genes and Signaling Pathways Associated With Feed Efficiency in Xiayan Chicken
- College of Animal Science and Technology, Guangxi University, Nanning, China
Feed efficiency is an important economic factor in poultry production, and the rate of feed efficiency is generally evaluated using residual feed intake (RFI). The molecular regulatory mechanisms of RFI remain unknown. Therefore, the objective of this study was to identify candidate genes and signaling pathways related to RFI using RNA-sequencing for low RFI (LRFI) and high RFI (HRFI) in the Xiayan chicken, a native chicken of the Guangxi province. Chickens were divided into four groups based on FE and sex: LRFI and HRFI for males and females, respectively. We identified a total of 1,015 and 742 differentially expressed genes associated with RFI in males and females, respectively. The 32 and 7 Gene Ontology (GO) enrichment terms, respectively, identified in males and females chiefly involved carbohydrate, amino acid, and energy metabolism. Additionally, Kyoto Encyclopedia of Genes and Genomes (KEGG) analysis identified 11 and 5 significantly enriched signaling pathways, including those for nutrient metabolism, insulin signaling, and MAPK signaling, respectively. Protein–protein interaction (PPI) network analysis showed that the pathways involving CAT, ACSL1, ECI2, ABCD2, ACOX1, PCK1, HSPA2, and HSP90AA1 may have an effect on feed efficiency, and these genes are mainly involved in the biological processes of fat metabolism and heat stress. Gene set enrichment analysis indicated that the increased expression of genes in LRFI chickens was related to intestinal microvilli structure and function, and to the fat metabolism process in males. In females, the highly expressed set of genes in the LRFI group was primarily associated with nervous system and cell development. Our findings provide further insight into RFI regulation mechanisms in chickens.
Poultry is one of the healthiest meat sources due to its low fat and high protein content, and its health factors have led to an increase in its consumption in recent years (Liu et al., 2020). This increase in demand requires an increase in chicken feed on farms, leading to increased production costs. The cost of feed already makes up most of the total cost of poultry production (Aggrey et al., 2010). Therefore, improving feed efficiency plays an essential role in the production of poultry products.
Feed efficiency is generally estimated using residual feed intake (RFI), which was proposed in 1963 and is considered the most functional parameter for the evaluation of feed efficiency (Koch et al., 1963). At present, RFI has been applied in the artificial selection of feed efficiency in mammals and poultry (Mebratie et al., 2019; Banerjee et al., 2020; Liu and VandeHaar, 2020). There is a general agreement that RFI is a moderately inherited characteristic, making it easy to improve the feed efficiency of commercial breeding companies (Sell-Kubiak et al., 2017). The current RFI breeding process is mainly based on genomic selection and is expensive (VanRaden, 2020). However, the genes and biological processes that account for feed efficiency remain largely unknown. Therefore, there is an urgent need to develop effective biomarkers to facilitate RFI selection.
With the development of next-generation sequencing, it is possible to investigate the mechanisms underlying RFI and accelerate the breeding process of broiler chicken feed efficiency using molecular bioinformatics (Mardis, 2008; Shendure and Ji, 2008). High-throughput sequencing techniques have become powerful and effective tools for gaining a deeper understanding of the basic molecular mechanisms of complicated systems (Wang et al., 2009; Ozsolak and Milos, 2011). RNA sequencing, a high-throughput technique, has been generally used in livestock to find the expression patterns of functional genes (Cui et al., 2017; Ni et al., 2019; Peng et al., 2019).
In chickens, the duodenum is considered the main feed absorption organ. A previous duodenal transcriptome study indicated that RFI differences may be associated with digestion, metabolism, biosynthetic processes, and energy homeostasis (Yi et al., 2015). In addition, many studies have demonstrated that feed efficiency has a particular influence on mitochondrial function, intestinal cellularity and absorption, appetite regulation, and fat metabolism (Montanholi et al., 2013; Lancaster et al., 2014; Perkins et al., 2014; Reyer et al., 2017).
To date, abundant sequencing analyses have been carried out on commercial broiler breeds, but only a few studies have focused on indigenous chickens (Sun et al., 2013; Wolc et al., 2020). The physiological characteristics and genetic backgrounds of indigenous and commercial chicken breeds are quite distinctive. Thus, gaining insight into genetic resources is necessary for genetic studies of indigenous breeds. The Xiayan chicken is an indigenous chicken of the Guangxi province, in southern China. Due to its excellent meat quality, the Xiayan chicken is likely to become one of the preferred poultry breeds for consumers in the Guangxi and Guangdong provinces. Hence, this study was designed to identify the candidate genes and signaling pathways related to feed efficiency through the transcriptional sequencing of the duodenum in male and female Xiayan chickens. This will contribute to uncovering the molecular mechanisms of feed efficiency of indigenous chickens in Guangxi.
Materials and Methods
All animal experiments and methods in this study have been evaluated and approved by the Animal Ethics Committee of Guangxi University (GXU2018-058).
Chicken Breed and RFI Calculation
According to the standard breed program, 340 indigenous Xiayan chicken (male = 173, female = 167) were bred in the Guangxi University experiment farm. The normal experiment started 70 to 90 days after the 10-day pre-experiment. The chickens were raised in 3-layer metal cages; the average stocking density was 15 birds per square meter. The basal diet was a corn-soybean broiler diet (13 MJ metabolizable energy/kg of diet, 220 g/kg crude protein) formulated without antibiotics or coccidiostats. Ad libitum access to fresh water was provided. Chickens were ranked by RFI, where the 10 most extreme duodenum samples from high (male = 3, female = 3) and low (male = 2, female = 2) RFI chickens were selected for RNA extraction. The feed intake (FI) and body weight (BW) were measured at 70–90 days of age. Feed conversion ratio (FCR) was calculated by FI and Body weight gain (BWG). The metabolic body weight (MBW0.75), BWG, and average daily body weight gain (ADG) were calculated according to the bodyweight of 70 and 90 days. ADFI was the average daily feed intake. The RFI value was used to measure the feed efficiency of Xiayan chickens, and we estimated it using the model as follows (Izadnia et al., 2019): RFI = ADFI-(b0 + b1MBW0.75 + b2ADG) where b0, b1, and b2 were regression intercept, the partial regression coefficient of ADFI on MBW0.75, and the partial regression coefficient of ADFI on ADG, respectively. SAS procedures t-test (SAS Version 9.4) was used to analyze the feed efficiency difference between HRFI and LRFI groups. The probability value was P < 0.05, indicating statistical significance. Detailed methods were described in our previous studies (Du et al., 2020).
RNA Extraction and RNA Sequencing
Total RNA was extracted using the TRIzolTM reagent (Invitrogen, Carlsbad, CA, United States) according to the manufacturer’s instructions. RNA purity and concentration were measured using the NanoPhotometer® spectrophotometer (IMPLEN, Westlake Village, CA, United States) and Qubit® RNA Assay Kit in Qubit® 2.0 Fluorometer (Life Technologies, Carlsbad, CA, United States). RNA integrity was assessed using the RNA Nano 6000 Assay Kit of the Bioanalyzer 2100 system (Agilent Technologies, Santa Clara, CA, United States). Subsequently, sequencing libraries were generated using the rRNA-depleted RNA by NEBNext® UltraTM Directional RNA Library Prep Kit for Illumina® (NEB, United States) following manufacturer’s recommendations. Fragmentation was carried out using divalent cations under elevated temperature in NEBNext First Strand Synthesis Reaction Buffer (5X). First-strand cDNA was synthesized using random hexamer primer and M-MuLV Reverse Transcriptase (RNaseH-). Second-strand cDNA synthesis was subsequently performed using DNA Polymerase I and RNase H. In the reaction buffer, dNTPs with dTTP were replaced by dUTP. Remaining overhangs were converted into blunt ends via exonuclease/polymerase activities. After adenylation of 3′ ends of DNA fragments, NEBNext Adaptor with hairpin loop structure was ligated to prepare for hybridization. In order to select cDNA fragments of preferentially 150∼200 bp in length, the library fragments were purified with the AMPure XP system (Beckman Coulter, Beverly, United States). Then 3 μl USER Enzyme (NEB, United States) was used with size-selected, adaptor-ligated cDNA at 37°C for 15 min followed by 5 min at 95°C before PCR. Then PCR was performed with Phusion High-Fidelity DNA polymerase, Universal PCR primers, and Index (X) Primer. Then, products were purified (AMPure XP system) and library quality was assessed on the Agilent Bioanalyzer 2100 system. At last the libraries were sequenced on an Illumina Hiseq 4000 platform, and 150 bp paired-end reads were generated.
RNA-Seq Data Analysis
Before read alignment, raw data (raw reads) in fastq format were first processed through Trimmomatic (Bolger et al., 2014). In this step, read bases with a phred base quality score less than 20, sequencing adapters, and reads containing poly-N are filtered to generate reliable clean data. Furthermore, Fastqc was used further to control the overall quality level of clean data to ensure the reliability of subsequent bioinformatics analysis (Schmieder and Edwards, 2011). The reference genome sequence files and annotation files used in this study were downloaded from the Ensembl genome browser1. Hisat2v2.1.0 was used to align clean data to a reference genome (Kim et al., 2015; Pertea et al., 2016). Afterward, stringtie (version 2.1.1) was used to assemble the transcriptome of each sample to generate a comprehensive transcript set (Pertea et al., 2015). With the fragments per kilobase of exon per million reads (FPKM) value, the gene expression levels can be quantified to a certain extent (Mortazavi et al., 2008). Furthermore, the assembly transcripts for all samples were integral to the enhancement of the overall quality of assembly by combining novel and mapped transcripts with a single one and removing the manual structures. The read count assignment was performed with the HTSeq-count tool of the HTSeq software (v0.6.1p1) (Anders et al., 2015). Prior to DEseq2 (version 2.2.1) (Love et al., 2014) read count standard normalization and expression analysis, genes with counts <1 were removed. Differently expressed genes (DEGs) were identified with | fold change| > 1 and p-value < 0.05. DEGs were picked out and analyzed in Kyoto Encyclopedia of Genes and Genomes (KEGG) and Gene Ontology (GO) by KOBAS software and Goseq in the R package. Benjamini-Hochberg (B-H) p-value < 0.05 was regarded as statistically significant.
Protein–Protein Interaction (PPI) Network Construction and Modules Selection
Differently expressed genes were submitted to the online search tool STRING database2 to search for interacting genes to obtain gene interaction relationships, with a confidence score >0.9 defined as significant (Szklarczyk et al., 2019). The open-source software Cytoscape (Shannon et al., 2003) was used to visualize the PPI network of DEGs, and to visualize complex networks that provide a wide range of applications to analyze the interactive network further. The application of Molecular Complex Detection (MCODE) (Bader and Hogue, 2003) in Cytoscape was utilized to screen the modules of the PPI network. The standard settings of MCODE rest with: degree cutoff = 2, node score cutoff = 0.2, k-core = 2, maximum depth = 100. Additionally, the GO and KEGG enrichment analysis for the function and pathway of genes were carried out in the modules.
Gene Set Enrichment Analysis (GSEA)
All expressed genes in both males and females could be applied in GSEA analysis by GSEA software, which ranked all expressed genes based on the fold-change (HRFI/LRFI) between the HRFI and LRFI groups (Bertocchi et al., 2019). The enrichment score of each gene set was then calculated with the full ranking, thus shedding light on the distribution of each gene set in the list. Furthermore, a normalized enriched score (NES) was determined for each gene set. The significant enrichment of the gene set followed the absolute values of NES > 1 and false discovery rate (FDR) ≤ 0.05 (Liu et al., 2018).
Validation of RNA-Seq
Six DEGs were randomly selected. The primers used for qPCR were designed using Oligo 6.0 software. Five micrograms of RNA were reverse-transcribed into cDNA using RT Reagent Kit (TaKaRa, Dalian, China). The volume of the reaction mixture was 20 μl, with 2 μl of cDNA, 0.5 μl for each primer, 10 μl of SYBR (TaKaRa, Dalian, China), and 7 μl of RNA-free water. The following RT-PCR reaction was performed as follows: 95°C for 3 min; 95°C for 10 s, annealing temperature for the 30 s for 35 cycles; finally, melting curve collection at 65 to 95°C. The expression levels were calculated according to the 2−ΔΔCt method normalized with β-actin. Primers used were synthesized by Sangon Biotech (Shanghai, China) and analyzed using Oligo 7.0 software. Primers for the samples KMT2E, PPA2, IQGAP2, NPLOC4, SSH1, and TRAFD1 genes are listed in Supplementary Table 1.
Performance and Feed Efficiency
The disparity in ADFI, MBW0.75, ADG, RFI, and FCR is illustrated in Table 1. As expected, the FCR and RFI of the LRFI group were considerably lower than those of the HRFI group (P < 0.05), and the ADFI of the LRFI group accounted for approximately 80% of the HRFI group. In females, the ADFI and ADG of HRFI birds were significantly higher than those of LRFI birds (P < 0.05). Whereas in males, the RFI value of LRFI birds was −13.03 ± 0.92 compared with 13.56 ± 1.47 for the HRFI birds for 20 experimental days (day 70–90). In females, the LRFI birds had an RFI value of −11.79 ± 2.55 compared with 11.89 ± 1.62 for the HRFI birds (day 70–90). Moreover, there was no significant difference in metabolic body weight (MBW0.75) between the two groups (p > 0.05).
Table 1. Characterization of performance and feed efficiency traits of male and female (Least square means and SEM).
All duodenal samples of males and females (total n = 10) were gathered for RNA-seq. The amount of raw reads, clean reads, total mapped reads (%), uniquely mapped reads (%), Q20(%), Q30(%), and GC content (%) for each sample is shown in Supplementary Table 2. After the filter, the number of clean reads per sample ranged from 50,275,148 to 102,421,430. The general Q30 percentage of clean data was above 91%. Comparing the sequencing reads with the chicken reference genome, we found that the total map rate is between 81.92 and 91.79%, and the unique map rate is between 74.43 and 81.42%. The GC content of 10 samples ranged from 48.45% to 50.34%.
Identification of DEGs
In this study, differential expression analysis was performed to detect gene expression differences between the HRFI and LRFI groups. A total of 1015 and 742 genes were identified as being DEGs in males and females, respectively. Of the 1015 DEGs, 381 DEGs were upregulated in the HRFI groups, while 634 were downregulated than the LRFI groups (Figure 1A). In the females, 481 DEGs were upregulated in the HRFI groups, while 261 were downregulated compared with the LRFI groups (Figure 1B). 57 DEGs are shared between females and males (Figure 1C).
Figure 1. Differently expressed genes (DEGs) of transcriptome analysis. (A,B) Represent the volcano map of DEGs in the male and female, respectively. (C) Venn diagram showing the intersection of the DEGs between males and females.
GO and KEGG Analysis
To further explore the functions of DEGs, we conducted a functional enrichment analysis. GO enrichment analysis showed that 8 and 11 GO items related to biological processes were significantly enriched in males and females, respectively. Notable among these were the metabolic process, cellular metabolic process, cellular metabolic process, and protein metabolic process. Other significant GO entries related to cellular components and molecular functions were oxidoreductase activity, mitochondrial part, enzyme binding cytoplasmic part, and membrane-bounded organelle (Figure 2).
Figure 2. Function enrichment analysis of differently expressed genes. (a) GO enrichment analysis of DEGs in the Male. (b) KEGG enrichment analysis of DEGs in the male. (A) GO enrichment analysis of DEGs in the female. (B) KEGG enrichment analysis of DEGs in the female.
Kyoto Encyclopedia of Genes and Genomes pathways enrichment analysis was performed to reveal the biological functions of DEGs further. In the male, 11 signaling pathways were significantly enriched, which mainly involved nutrient metabolism, energy metabolism, and insulin signaling pathways. Genes related to nutrient metabolism and energy metabolism were upregulated in the LRFI group, including UQCRFS1, PCK1, ALDOB, EHHADH, and GCH1. Simultaneously, genes involved in the insulin signaling pathway and calcium signaling pathway were upregulated in the HRFI group, including VDAC1, PHKA2, PHKG1, PTGFR, and CACNA1C. In the female, 5 pathways were significantly enriched, which were involved in endoplasmic reticulum protein processing and actin cells. The genes related to skeleton regulation and MAPK signaling pathway were upregulated in LRFI groups, while the genes involved in small molecule metabolism and synthesis were upregulated in the HRFI group.
Identification of Hub Genes and Pathways Through PPI Network Analysis
Protein–protein interaction network analysis was performed on DEGs; we could get more insight into the interaction relationship among them. The interaction relationship of DEGs was shown in Figure 3. We built the top three critical modules in DEGs’ PPI network through the MCODE application (Figure 4). And then, GO and KEGG enrichment analysis were performed on these genes in the modules. In the male, the top three modules’ genes were significantly enriched in the peptide biosynthetic process, peptide metabolic process, fatty acid degradation, peroxisome, fatty acid metabolism, and PPAR signaling pathway. In the female, these genes’ functions may be explained by the cellular response to heat, aminoacyl-tRNA biosynthesis, and ribosome. A complete result of the enrichment analysis of genes in each module was shown in Supplementary Table 3.
Figure 3. Protein–protein interaction (PPI) network analysis. (A) PPI network for DEGs in males. (B) PPI network for DEGs in the female. Red circles represent upregulated genes, and green circles represent downregulated genes. The size of the circle indicates the fold change of each gene (HRFI/LRFI).
Figure 4. The top three protein–protein interaction (PPI) hub network modules in the male (A–C) and female (a–c). Red circles represent upregulated genes, and green circles represent downregulated genes. The size of the circle indicates the fold change of each gene (HRFI/LRFI).
We further investigated the difference in gene expression levels between HRFI and LRFI groups by GSEA. The results of the GSEA analysis were presented in Tables 2–5. As for the GO-based list, in the male, higher expression gene sets in the LRFI group were mainly associated with intestinal digestion and absorption, such as brush border, brush border membrane, and intestinal absorption. KEGG-base gene set enrichment analysis enriched primarily for xenobiotics’ metabolism by cytochrome p450, fatty acid metabolism, drug metabolism cytochrome p450, peroxisome, and PPAR signaling pathway. From the GO-based list, in the female, higher expression gene sets in the LRFI group were mainly connected to neurodevelopment, such as neuron fate commitment, central nervous system neuron differentiation, and neuron fate specification. The KEGG-base Gene set enriched primarily for basal cell carcinoma, hedgehog signaling pathway, neuroactive ligand-receptor interaction, and melanogenesis. The most enriched GO and KEGG items in the male and female were shown in Figure 5.
Figure 5. Gene set enrichment analysis (GO-base and KEGG-base). GSEA algorithm scored the enrichment of genes in the pathway in the ranked gene list. Enrichment score (ES) > 0 indicates that the distribution of the gene set is biased upstream of the ranking list, and ES < 0 shows the gene set distribution is biased downstream of the ranking list. (A,B), (a,b) represent the GO entries and KEGG signaling pathways with the lowest FDR values. (A) Brush border and (B) Metabolism of xenobiotics by cytochrome p450 are enriched in the male LRFI groups. (a) Cell fate commitment and (b) Basal cell carcinoma are enriched in LRFI groups in the female.
Validation of RNA-Seq
Six genes were selected randomly from DEGs for qPCR validation. The results showed that KMT2E, PPA2, and IQGAP2 were upregulated in the LRFI group, while TRAFD1, NPLOC4, and SSH1 were upregulated in the HRFI group. The expression patterns of six genes in qPCR were consistent with those in RNA-seq (Figure 6).
Figure 6. Validation of the accuracy of RNA-seq. (A) Correlations of the expression level of 6 random DEGs between high and low abdominal fat using RNA-Seq and qPCR. The x- and y-axis correspond to the genes and log2 (ratio of H/L) measured by RNA-Seq and qPCR. (B) Gene expression abundance of HRFI and LRFI in Male and female; results are expressed as means ± standard deviation (n = 6); a, b means p < 0.05.
Feed efficiency plays an important role in improving profits and the environmental footprint in broiler production. In this study, the duodenum transcriptome data came from four groups of Xiayan chickens with extreme opposing RFI and different sex using RNA-seq. The gene expression profile was further explored by differential expression analysis, GO and KEGG enrichment analysis, PPI network analysis, and GSEA. All bioinformatics analyses were conducted to study gene expression differences, associations, and enrichment to further gain more widespread biological insight into indigenous chickens’ feed efficiency.
Among males and females, the values of RFI, FI, and FCR of high-RFI chickens were significantly higher than those of low-RFI chickens. However, there was no significant difference between ADG, BW, and MBW, which was consistent with the findings of previous studies (Zhang X. et al., 2017; Zhang et al., 2019). RFI uses feed intake to measure feed efficiency; using RFI to select feed efficiency in breeding work can avoid affecting other growth traits. RFI could be used as an ideal breeding index for chicken feed efficiency trait.
Traditional differential expression analysis of RNA sequencing data would produce many DEGs, and further analysis was needed to understand the function of the DEGs between different comparison groups (Hekman et al., 2015). Therefore, we further performed functional enrichment analysis on these DEGs. The analysis results of GO and KEGG showed that these genes were mainly involved in metabolism progress, including material metabolism processes, protein metabolism, carbohydrate metabolism, amino acid metabolism, energy metabolism, and lipid metabolism processes. Many studies have shown that metabolism is an essential factor affecting residual feed intake (Jing et al., 2015; Yi et al., 2015; Kong et al., 2016).
In males, some genes related to the metabolism process were upregulated in the LRFI group, such as ALDOB (Aldolase, Fructose-Bisphosphate B), UQCRFS1 (ubiquinone-cytochrome C reductase iron-sulfur protein subunit Base 1), EHHADH (3-hydroxyacyl-CoA dehydrogenase), and GCH1 (guanosine triphosphate cyclase I). A notable gene is ALDOB. The aldolase family is a biological enzyme necessary for bioenergy metabolism. It catalyzes the formation of dihydroxyacetone phosphate and glyceraldehyde 3-phosphate by fructose 1, 6-bisphosphate, which plays a vital role in glycolysis (Almon et al., 2007). ALDOB (aldolase B), as a member of the aldolase family, is mainly distributed in liver tissues, participates in liver metabolism, and also participates in glycolysis and gluconeogenesis processes (Kallio et al., 1998). More recent studies already point to ALDOB as a candidate gene for feed efficiency in chickens. A previous study in chicken found the ALDOB gene that participates in glycolysis and gluconeogenesis was highly expressed in LRFI chickens (Shah et al., 2019). Gene expression in breast muscle associated with feed efficiency in a single male broiler line using a chicken 44K oligo microarray also identified that ALDOB was more highly expressed in high feed efficiency chickens (Kong et al., 2011). We speculate that ALDOB could have an essential role in enhancing energy metabolism in the high-FE broiler phenotype. Some genes enriched in the insulin signaling pathway were upregulated in the HRFI group. In the chicken, the insulin-signaling pathway has anabolic effects in glucose transport and utilization, glycogen synthesis, control of liver lipogenic enzymes, amino acid transport, and protein synthesis (Fujita et al., 2019). Previous genome-wide association analysis identified an SNP site significantly related to feed efficiency that may regulate feed efficiency through the insulin signaling pathway (Yuan et al., 2017). A study conducted on Duroc pigs with differences in RFI has found that the insulin signaling pathway may be a biological process that is significantly associated with RFI (Banerjee et al., 2020). Our findings were consistent with the results above, indicating the body’s digestion and metabolism are essential to RFI, and high feed efficiency chickens may have a stronger ability to metabolize nutrients. In the female, genes related to the MAPK signaling pathway and actin cytoskeleton regulation were upregulated in the LRFI group. Many studies have proven that MAPK plays a key role in the regulation of energy balance (Kikuchi et al., 2020; Okubo et al., 2020). Previous studies have shown that the MAPK signaling pathway is widely involved in the regulation of growth and development (Xu et al., 2019; Zhao et al., 2019); the genes contained in the MAPK signaling pathway may regulate feed efficiency through energy distribution and homeostasis in the digestion, absorption, and metabolism of feed.
Protein–protein interaction network analysis can capture the interaction information of DEGs, which helps to understand the molecular mechanism of fat deposition from the perspective of biological systems (Athanasios et al., 2017). In this research, the PPI network analysis was constructed with DEGs. To identify the essential part of the PPI network, we used the “module” program to extract the top three modules and performed GO and KEGG enrichment analysis. Some GO items and KEGG signaling pathways were worth noting, including fatty acid degradation, fatty acid metabolism, PPAR signaling pathway, and cellular response to heat. Some biological pathways, like lipid metabolism and cholesterol biosynthesis, were identified to be associated with RFI (Nafikov and Beitz, 2007; Karisa et al., 2014). Our research found that some key hub-genes are significantly enriched in fat metabolism-related pathways. These genes include ECH1, EHHADH, CAT, ACSL1, ECI2, ABCD2, ACOX1, and PCK1. Interestingly, these genes are both highly expressed in the high feed efficiency group in male/females. ACSL1 (acyl-CoA synthetase long-chain family member 1) plays an important role in the transportation and activation of fatty acids. A previous study in chickens showed that high expression of the ACSL1 gene could promote fat synthesis, while chickens with high feed efficiency tend to have more fat deposits (Neijat et al., 2017). The PCK1 gene is associated with obesity, insulin resistance, type II diabetes in mammals, and abdominal fat content (Beale et al., 2007; Rees et al., 2009; Millward et al., 2010). Previous research showed that fat deposition was positively correlated with PCK1 expression in birds.
Similarly, a prior study in chicken performance transcriptome sequencing also found PCK1 higher express in the LRFI group than HRFI; it was also speculated that LRFI chickens have more fatty deposits (Duan et al., 2013), which is consistent with the results in this study. A previous genome-wide association analysis related to pig feed efficiency found a significant correlation between the ECI2 (enoyl-CoA delta isomerase 2) gene and feed efficiency. Earlier studies in pigs showed that the expression level of the ABCD2 (ATP binding cassette subfamily D member 2) gene was downregulated in the high feed efficiency group, and this result was also verified in our study. In a meat duck liver transcriptome study, ACOX1 (acyl-CoA oxidase 1) expression was significantly negatively correlated with RFI, which was also consistent with our results. As for ECH1, EHHADH, and CAT genes, no relevant studies have shown that they have a relationship with RFI. It is challenging to predict chickens’ capability with low or high RFI to cause specific fat deposition changes. But our research confirmed that the regulatory relationship between feed efficiency and fat deposition needs to be further explored.
Interestingly, two key hub genes encode for Heat Shock Proteins (HSPs): HSPA2/HSP90AA1. Chickens are susceptible to heat stress due to their overall feather coverage and lack of sweat glands (Zhang C. et al., 2017). Numerous research has reported on the negative influence of heat stress on poultry production, such as decreased body antioxidant capacity and intestinal immunity and impaired intestinal morphology (Sahin et al., 2017; Song et al., 2018). Under HS conditions, chickens may spend more energy on maintenance and acclimation, which reduces the energy for growth and leads to a decrease in BWG (Mujahid et al., 2007). The activation of heat stress is energetically costly, and long-term stimulation can negatively impact feed efficiency. These genes involved in heat stress can be used as candidate markers related to feed efficiency. HSPs protein is highly conserved in all organisms and plays a crucial role in cellular stress response (Stetler et al., 2010). Many studies have demonstrated the link between heat stress and feed affection. Heat stress can lead to a decrease in animal feed intake and damage the intestine’s integrity and barrier function (Pearce et al., 2015). Studies have reported that when pigs undergo heat stress, the expression of HSP mRNA will increase. The abundance of various metabolic enzymes in the ileum will decrease, suggesting that the body’s metabolic process has changed (Pearce et al., 2015). Similar to our findings, a study explored the host transcriptome and microbiome interaction modulates of full-sibs broilers with divergent feed conversion ratio and identified that HSP90AA1 is a significantly differently expressed gene. In poultry production, strictly controlling the environment and avoiding heat stress will effectively improve feed efficiency.
We used the GSEA approach (Wang and Cairns, 2013) to compare biological differences in all gene expression patterns between the high-RFI and low-RFI groups. Based on the GO-based list, all higher expressed gene sets in the LRFI group were mainly divided into three categories. lipid metabolism, structure and function of small intestinal villi, and development of the nervous system. Some genes are enriched to the gene set involved in border brush, which is considered to be related to the intestine’s digestion and absorption function, including MYO1D, MYO1E, MYO1A, USH1C, and EZR. These genes in this pathway are highly expressed in the high feed efficiency group, which shows that some genes related to the intestinal structure may play an important role in regulating high feed efficiency chicken intestines with better absorption capacity. Biologists have long appreciated the intimate connection between morphology and function. Here, morphological adaptations at both the tissue and cellular level, allow the intestinal epithelium to make close and prolonged contact with luminal contents, promoting efficient uptake of available nutrients that may lead to the change of feed efficiency. The intestinal brush border is home to several class I myosins, with myosin-1a (myo1a) being the most abundant (McConnell et al., 2011). The expression of the myo1A gene has strong tissue specificity. Expression of myo1A is limited to the intestinal tract, where it localizes almost exclusively to the brush border (Skowron et al., 1998; Skowron and Mooseker, 1999). MYO1A, MYO1D, and MYO1E are considered to be key genes in the rapid differentiation of neonatal epithelial cells into mature intestinal epithelial cells (Tyska et al., 2005; Beale et al., 2007; Benesh et al., 2010). Previous studies have shown that a large deletion mutation in the human USH1C gene can cause severe gastrointestinal dysfunction (Bitner-Glindzicz et al., 2000; Hussain et al., 2004). Consistent with this, USH1C knockout mice, which were developed to model Type 1 Usher syndrome, display significant perturbations in intestinal brush border morphology (Crawley et al., 2014). Our results show the potential connection between chicken intestinal structure and function and feed efficiency. It is worth noting that in our GSEA enrichment analysis, some genes were enriched and related to nervous system development. A previous study integrated genome-wide co-association and gene expression found that some candidate genes directly or indirectly affect FE-related traits were mainly associated with immunity, nervous system, behavior, and energy metabolism (Ramayo-Caldas et al., 2019). The influence mechanism of nervous system development on feed efficiency needs further study.
To summarize, this work is the first to use RNA-seq in males and females to identify and annotate DEGs in the duodenum tissues of high- and low-FE native chickens. In male and females, we identified 1015 and 742 DEGs associated with RFI, respectively. Those genes were mainly enriched in the pathways of metabolism, oxidation reaction, and energy homeostasis. Moreover, we found that genes involved in fat metabolism and heat stress may affect feed efficiency through PPI network analysis. Through GSEA analysis, we found that many genes related to small intestinal villi structure and absorption function can be used as candidate genes to RFI, including MYO1D, MYO1E, MYO1A, USH1C, and EZR. Finally, the research shows that differences in gene expression patterns related to nervous system development may cause different feed efficiencies.
Data Availability Statement
The sequencing data have been deposited in the NCBI SRA (http://www.ncbi.nlm.nih.gov/sra), BioProject accession PRJNA672913.
The animal study was reviewed and approved by the Animal Ethics Committee of Guangxi University (GXU2018-058).
CX, JD, and ZY: data curation. CX, TS, ZY, and XY: formal analysis. XY: funding acquisition, project administration, and writing – review and editing. CX, JD, and LZ: investigation. CX and XY: methodology. JD and XY: resources. ZY and JD: supervision. TS and LZ: validation. CX and JD: writing – original draft. All authors contributed to the article and approved the submitted version.
This work was supported in part by the Science and Technology Major Project of Guangxi (GK AA17204027).
Conflict of Interest
The 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.
The Supplementary Material for this article can be found online at: https://www.frontiersin.org/articles/10.3389/fgene.2021.607719/full#supplementary-material
Almon, R. R., DuBois, D. C., and Jusko, W. J. (2007). A microarray analysis of the temporal response of liver to methylprednisolone: a comparative analysis of two dosing regimens. Endocrinology 148, 2209–2225. doi: 10.1210/en.2006-0790
Athanasios, A., Charalampos, V., Vasileios, T., and Ashraf, G. M. (2017). Protein-protein interaction (PPI) network: recent advances in drug discovery. Curr. Drug Metab. 18, 5–10. doi: 10.2174/138920021801170119204832
Banerjee, P., Carmelo, V. A. O., and Kadarmideen, H. N. (2020). Genome-wide epistatic interaction networks affecting feed efficiency in duroc and landrace pigs. Front. Genet. 11:121. doi: 10.3389/fgene.2020.00121
Benesh, A. E., Nambiar, R., McConnell, R. E., Mao, S., Tabb, D. L., and Tyska, M. J. (2010). Differential localization and dynamics of class I myosins in the enterocyte microvillus. Mol. Biol. Cell 21, 970–978. doi: 10.1091/mbc.E09-07-0638
Bertocchi, M., Sirri, F., Palumbo, O., Luise, D., Maiorano, G., Bosi, P., et al. (2019). Exploring Differential Transcriptome between Jejunal and Cecal Tissue of Broiler Chickens. Animals 9:221. doi: 10.3390/ani9050221
Bitner-Glindzicz, M., Lindley, K. J., Rutland, P., Blaydon, D., Smith, V. V., Milla, P. J., et al. (2000). A recessive contiguous gene deletion causing infantile hyperinsulinism, enteropathy and deafness identifies the Usher type 1C gene. Nat. Genet. 26, 56–60. doi: 10.1038/79178
Crawley, S. W., Shifrin, D. A., Grega-Larson, N. E., McConnell, R. E., Benesh, A. E., Mao, S. L., et al. (2014). Intestinal brush border assembly driven by protocadherin-based intermicrovillar adhesion. Cell 157, 433–446. doi: 10.1016/j.cell.2014.01.067
Cui, X., Marshall, B., Shi, N., Chen, S. Y., Rekaya, R., and Liu, H. X. (2017). RNA-Seq analysis on chicken taste sensory organs: an ideal system to study organogenesis. Sci. Rep. 7:9131. doi: 10.1038/s41598-017-09299-7
Du, W., Deng, J., Yang, Z., Zeng, L., and Yang, X. J. P. S. (2020). Metagenomic analysis reveals linkages between cecal microbiota and feed efficiency in Xiayan chickens. Poult. Sci. 99, 7066–7075. doi: 10.1016/j.psj.2020.09.076
Duan, J., Shao, F., Shao, Y., Li, J., Ling, Y., Teng, K., et al. (2013). Androgen inhibits abdominal fat accumulation and negatively regulates the PCK1 gene in male chickens. PLoS One 8:e59636. doi: 10.1371/journal.pone.0059636
Fujita, S., Honda, K., Yamaguchi, M., Fukuzo, S., Saneyasu, T., and Kamisoyama, H. (2019). Role of insulin-like growth factor-1 in the central regulation of feeding behavior in chicks. J. Poult. Sci. 56, 270–276. doi: 10.2141/jpsa.0180127
Hekman, J. P., Johnson, J. L., and Kukekova, A. V. (2015). Transcriptome analysis in domesticated species: challenges and strategies. Bioinform. Biol. Insights 9(Suppl. 4), 21–31. doi: 10.4137/BBI.S29334
Hussain, K., Bitner-Glindzicz, M., Blaydon, D., Lindley, K. J., Thompson, D. A., Kriss, T., et al. (2004). Infantile hyperinsulinism associated with enteropathy, deafness and renal tubulopathy: clinical manifestations of a syndrome caused by a contiguous gene deletion located on chromosome 11p. J. Pediatr. Endocrinol. Metab. 17, 1613–1621. doi: 10.1515/jpem.2004.17.12.1613
Izadnia, H. R., Tahmoorespur, M., Bakhtiarizadeh, M. R., Nassiri, M., and Esmaeilkhanien, S. (2019). Gene expression profile analysis of residual feed intake for Isfahan native chickens using RNA-SEQ data. Ital. J. Anim. Sci. 18, 246–260. doi: 10.1080/1828051x.2018.1507625
Jing, L., Hou, Y., Wu, H., Miao, Y., Li, X., Cao, J., et al. (2015). Transcriptome analysis of mRNA and miRNA in skeletal muscle indicates an important network for differential Residual Feed Intake in pigs. Sci. Rep. 5:11953. doi: 10.1038/srep11953
Kallio, P. J., Okamoto, K., O’Brien, S., Carrero, P., Makino, Y., Tanaka, H., et al. (1998). Signal transduction in hypoxic cells: inducible nuclear translocation and recruitment of the CBP/p300 coactivator by the hypoxia-inducible factor-1alpha. EMBO J. 17, 6573–6586. doi: 10.1093/emboj/17.22.6573
Karisa, B., Moore, S., and Plastow, G. (2014). Analysis of biological networks and biological pathways associated with residual feed intake in beef cattle. Anim. Sci. J. 85, 374–387. doi: 10.1111/asj.12159
Kikuchi, S., Piraino, G., O’Connor, M., Wolfe, V., Ridings, K., Lahni, P., et al. (2020). Hepatocyte-specific deletion of AMPKalpha1 results in worse outcomes in mice subjected to sepsis in a sex-specific manner. Front. Immunol. 11:210. doi: 10.3389/fimmu.2020.00210
Kong, B. W., Song, J. J., Lee, J. Y., Hargis, B. M., Wing, T., Lassiter, K., et al. (2011). Gene expression in breast muscle associated with feed efficiency in a single male broiler line using a chicken 44K oligo microarray. I. Top differentially expressed genes. Poultry Sci. 90, 2535–2547. doi: 10.3382/ps.2011-01435
Kong, R. S., Liang, G., Chen, Y., Stothard, P., and Guan le, L. (2016). Transcriptome profiling of the rumen epithelium of beef cattle differing in residual feed intake. BMC Genomics 17:592. doi: 10.1186/s12864-016-2935-4
Lancaster, P. A., Carstens, G. E., Michal, J. J., Brennan, K. M., Johnson, K. A., and Davis, M. E. (2014). Relationships between residual feed intake and hepatic mitochondrial function in growing beef cattle. J. Anim. Sci. 92, 3134–3141. doi: 10.2527/jas.2013-7409
Liu, E., and VandeHaar, M. J. (2020). Relationship of residual feed intake and protein efficiency in lactating cows fed high- or low-protein diets. J. Dairy Sci. 103, 3177–3190. doi: 10.3168/jds.2019-17567
Liu, J., Puolanne, E., Schwartzkopf, M., and Arner, A. (2020). Altered sarcomeric structure and function in woody breast myopathy of avian pectoralis major muscle. Front. Physiol. 11:287. doi: 10.3389/fphys.2020.00287
Liu, Z., Meng, J., Li, X., Zhu, F., Liu, T., Wu, G., et al. (2018). Identification of Hub Genes and Key Pathways Associated with Two Subtypes of Diffuse Large B-Cell Lymphoma Based on Gene Expression Profiling via Integrated Bioinformatics. Biomed. Res. Int. 2018:3574534. doi: 10.1155/2018/3574534
McConnell, R. E., Benesh, A. E., Mao, S. L., Tabb, D. L., and Tyska, M. J. (2011). Proteomic analysis of the enterocyte brush border. Am. J. Physiol. Gastrointest. Liver Physiol. 300, G914–G926. doi: 10.1152/ajpgi.00005.2011
Mebratie, W., Madsen, P., Hawken, R., Rome, H., Marois, D., Henshall, J., et al. (2019). Genetic parameters for body weight and different definitions of residual feed intake in broiler chickens. Genet. Sel. Evol. 51:53. doi: 10.1186/s12711-019-0494-2
Millward, C. A., DeSantis, D., Hsieh, C. W., Heaney, J. D., Pisano, S., Olswang, Y., et al. (2010). Phosphoenolpyruvate carboxykinase (Pck1) helps regulate the triglyceride/fatty acid cycle and development of insulin resistance in mice. J. Lipid Res. 51, 1452–1463. doi: 10.1194/jlr.M005363
Montanholi, Y., Fontoura, A., Swanson, K., Coomber, B., Yamashiro, S., and Miller, S. (2013). Small intestine histomorphometry of beef cattle with divergent feed efficiency. Acta Vet. Scand. 55:9. doi: 10.1186/1751-0147-55-9
Mujahid, A., Akiba, Y., and Toyomizu, M. (2007). Acute heat stress induces oxidative stress and decreases adaptation in young white leghorn cockerels by downregulation of avian uncoupling protein. Poult. Sci. 86, 364–371. doi: 10.1093/ps/86.2.364
Neijat, M., Eck, P., and House, J. D. (2017). Impact of dietary precursor ALA versus preformed DHA on fatty acid profiles of eggs, liver and adipose tissue and expression of genes associated with hepatic lipid metabolism in laying hens. Prostaglandins Leukot Essent Fatty Acids 119, 1–17. doi: 10.1016/j.plefa.2017.01.010
Ni, L., Song, C., Wu, X., Zhao, X., Wang, X., Li, B., et al. (2019). RNA-seq transcriptome profiling of porcine lung from two pig breeds in response to Mycoplasma hyopneumoniae infection. PeerJ 7:e7900. doi: 10.7717/peerj.7900
Pearce, S. C., Lonergan, S. M., Huff-Lonergan, E., Baumgard, L. H., and Gabler, N. K. (2015). Acute heat stress and reduced nutrient intake alter intestinal proteomic profile and gene expression in pigs. PLoS One 10:e0143099. doi: 10.1371/journal.pone.0143099
Peng, L. Y., Cui, Z. Q., Wu, Z. M., Fu, B. D., Yi, P. F., and Shen, H. Q. (2019). RNA-seq profiles of chicken type II pneumocyte in response to Escherichia coli infection. PLoS One 14:e0217438. doi: 10.1371/journal.pone.0217438
Perkins, S. D., Key, C. N., Garrett, C. F., Foradori, C. D., Bratcher, C. L., Kriese-Anderson, L. A., et al. (2014). Residual feed intake studies in Angus-sired cattle reveal a potential role for hypothalamic gene expression in regulating feed efficiency. J. Anim. Sci. 92, 549–560. doi: 10.2527/jas.2013-7019
Pertea, M., Kim, D., Pertea, G. M., Leek, J. T., and Salzberg, S. L. (2016). Transcript-level expression analysis of RNA-seq experiments with HISAT, StringTie and Ballgown. Nat. Protoc. 11, 1650–1667. doi: 10.1038/nprot.2016.095
Pertea, M., Pertea, G. M., Antonescu, C. M., Chang, T. C., Mendell, J. T., and Salzberg, S. L. (2015). StringTie enables improved reconstruction of a transcriptome from RNA-seq reads. Nat. Biotechnol. 33, 290–295. doi: 10.1038/nbt.3122
Ramayo-Caldas, Y., Marmol-Sanchez, E., Ballester, M., Sanchez, J. P., Gonzalez-Prendes, R., Amills, M., et al. (2019). Integrating genome-wide co-association and gene expression to identify putative regulators and predictors of feed efficiency in pigs. Genet. Select. Evol. 51:48. doi: 10.1186/s12711-019-0490-6
Rees, S. D., Britten, A. C., Bellary, S., O’Hare, J. P., Kumar, S., Barnett, A. H., et al. (2009). The promoter polymorphism-232C/G of the PCK1 gene is associated with type 2 diabetes in a UK-resident South Asian population. BMC Med. Genet. 10:83. doi: 10.1186/1471-2350-10-83
Reyer, H., Shirali, M., Ponsuksili, S., Murani, E., Varley, P. F., Jensen, J., et al. (2017). Exploring the genetics of feed efficiency and feeding behaviour traits in a pig line highly selected for performance characteristics. Mol. Genet. Genomics 292, 1001–1011. doi: 10.1007/s00438-017-1325-1
Sahin, N., Hayirli, A., Orhan, C., Tuzcu, M., Akdemir, F., Komorowski, J. R., et al. (2017). Effects of the supplemental chromium form on performance and oxidative stress in broilers exposed to heat stress. Poult. Sci. 96, 4317–4324. doi: 10.3382/ps/pex249
Sell-Kubiak, E., Wimmers, K., Reyer, H., and Szwaczkowski, T. (2017). Genetic aspects of feed efficiency and reduction of environmental footprint in broilers: a review. J. Appl. Genet. 58, 487–498. doi: 10.1007/s13353-017-0392-7
Shah, T. M., Patel, J. G., Gohil, T. P., Blake, D. P., and Joshi, C. G. (2019). Host transcriptome and microbiome interaction modulates physiology of full-sibs broilers with divergent feed conversion ratio. NPJ Biofilms Microb. 5:24. doi: 10.1038/s41522-019-0096-3
Shannon, P., Markiel, A., Ozier, O., Baliga, N. S., Wang, J. T., Ramage, D., et al. (2003). Cytoscape: a software environment for integrated models of biomolecular interaction networks. Genome Res. 13, 2498–2504. doi: 10.1101/gr.1239303
Song, Z. H., Cheng, K., Zheng, X. C., Ahmad, H., Zhang, L. L., and Wang, T. (2018). Effects of dietary supplementation with enzymatically treated Artemisia annua on growth performance, intestinal morphology, digestive enzyme activities, immunity, and antioxidant capacity of heat-stressed broilers. Poult. Sci. 97, 430–437. doi: 10.3382/ps/pex312
Stetler, R. A., Gan, Y., Zhang, W., Liou, A. K., Gao, Y., Cao, G., et al. (2010). Heat shock proteins: cellular and molecular mechanisms in the central nervous system. Prog. Neurobiol. 92, 184–211. doi: 10.1016/j.pneurobio.2010.05.002
Sun, Y. F., Zhao, G. P., Liu, R. R., Zheng, M. Q., Hu, Y. D., Wu, D., et al. (2013). The identification of 14 new genes for meat quality traits in chicken using a genome-wide association study. BMC Genomics 14:458. doi: 10.1186/1471-2164-14-458
Szklarczyk, D., Gable, A. L., Lyon, D., Junge, A., Wyder, S., Huerta-Cepas, J., et al. (2019). STRING v11: protein-protein association networks with increased coverage, supporting functional discovery in genome-wide experimental datasets. Nucleic Acids Res. 47, D607–D613. doi: 10.1093/nar/gky1131
Tyska, M. J., Mackey, A. T., Huang, J. D., Copeland, N. G., Jenkins, N. A., and Mooseker, M. S. (2005). Myosin-1a is critical for normal brush border structure and composition. Mol. Biol. Cell 16, 2443–2457. doi: 10.1091/mbc.e04-12-1116
Wang, X., and Cairns, M. J. (2013). Gene set enrichment analysis of RNA-Seq data: integrating differential expression and splicing. BMC Bioinformatics 14(Suppl. 5):S16. doi: 10.1186/1471-2105-14-S5-S16
Wolc, A., Drobik-Czwarno, W., Jankowski, T., Arango, J., Settar, P., Fulton, J. E., et al. (2020). Accuracy of genomic prediction of shell quality in a White Leghorn line. Poult. Sci. 99, 2833–2840. doi: 10.1016/j.psj.2020.01.019
Xu, Y., An, J. J., Tabys, D., Xie, Y. D., Zhao, T. Y., Ren, H. W., et al. (2019). Effect of lactoferrin on the expression profiles of long non-coding RNA during osteogenic differentiation of bone marrow mesenchymal stem cells. Int. J. Mol. Sci. 20:4834. doi: 10.3390/ijms20194834
Yi, G., Yuan, J., Bi, H., Yan, W., Yang, N., and Qu, L. (2015). In-depth duodenal transcriptome survey in chickens with divergent feed efficiency using RNA-Seq. PLoS One 10:e0136765. doi: 10.1371/journal.pone.0136765
Yuan, J., Chen, S., Shi, F., Wu, G., Liu, A., Yang, N., et al. (2017). Genome-wide association study reveals putative role of gga-miR-15a in controlling feed conversion ratio in layer chickens. BMC Genomics 18:699. doi: 10.1186/s12864-017-4092-9
Zhang, C., Zhao, X. H., Yang, L., Chen, X. Y., Jiang, R. S., Jin, S. H., et al. (2017). Resveratrol alleviates heat stress-induced impairment of intestinal morphology, microflora, and barrier integrity in broilers. Poult. Sci. 96, 4325–4332. doi: 10.3382/ps/pex266
Zhang, D., Zhang, X., Li, F., Li, C., La, Y., Mo, F., et al. (2019). Transcriptome analysis identifies candidate genes and pathways associated with feed efficiency in hu sheep. Front. Genet 10:1183. doi: 10.3389/fgene.2019.01183
Zhang, X., Wang, W., Mo, F., La, Y., Li, C., and Li, F. (2017). Association of residual feed intake with growth and slaughtering performance, blood metabolism, and body composition in growing lambs. Sci. Rep. 7:12681. doi: 10.1038/s41598-017-13042-7
Keywords: RNA-seq, feed efficiency, candidate genes, protein–protein interaction, gene set enrichment analysis
Citation: Xiao C, Deng J, Zeng L, Sun T, Yang Z and Yang X (2021) Transcriptome Analysis Identifies Candidate Genes and Signaling Pathways Associated With Feed Efficiency in Xiayan Chicken. Front. Genet. 12:607719. doi: 10.3389/fgene.2021.607719
Received: 18 September 2020; Accepted: 25 February 2021;
Published: 17 March 2021.
Edited by:Aline Silva Mello Cesar, University of São Paulo, Brazil
Reviewed by:Romdhane Rekaya, University of Georgia, United States
Gabriel Costa Monteiro Moreira, University of São Paulo, Brazil
Copyright © 2021 Xiao, Deng, Zeng, Sun, Yang and Yang. 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: Xiurong Yang, firstname.lastname@example.org
†These authors have contributed equally to this work