Yuan Chai1,2
Zaixia Liu2,3
Shaoyin Fu3,4Bin Liu5
Lili Guo2
Lingli Dai2
Yanyong Sun6Wenguang Zhang2,7*Chun Li8*Taodi Liu9*- 1College of Agronomy, Animal Husbandry and Bioengineering, Xing’an Vocational and Technical College, Ulanhot, China
- 2College of Animal Science, Inner Mongolia Agricultural University, Hohhot, China
- 3Inner Mongolia Academy of Agricultural and Animal Husbandry Science, Hohhot, China
- 4Institute of ATCG, Nei Mongol Bio-Information, Hohhot, China
- 5Nei Mongol BioNew Technology Co., Ltd., Hohhot, China
- 6College of Animal Science and Veterinary Medicine, Tianjin Agricultural University, Tianjin, China
- 7Kunming Institute of Zoology, Chinese Academy of Sciences, Kunming, China
- 8College of Animal Science and Technology, Inner Mongolia University for Nationalities, Tongliao, China
- 9School of Basic Medical Sciences, Inner Mongolia Medical University, Hohhot, China
The interplay between melatonin and immune system is well recognized in humans. The true integration of research on cashmere goat is still far from clear, especially for cashmere goat maintained in wool and cashmere growth. In this study, we applied various approaches to identify the complex regulated network between the immune-related genes and transcription factors (TFs) and to explore the relationship between melatonin and gene expression in cashmere goats. In total, 1,599 and 1756 immune-related genes were found in the blood and skin of cashmere goats, respectively, and 24 differentially expressed immune-related GO terms were highly expressed in blood after melatonin implantation. We studied the melatonin-dependent networks between the TFs and immune-related genes in cashmere goat. The 3 major regulatory networks were interconnected through TFs. The TFs, such as PHF5A, REXO4, STRAP, JUNB, GATAD2A, ZNF710, and VDR, were also expressed in the blood and skin tissue of cashmere goat. In addition, most genes in these networks, such as VDR, JUNB, and Trib3, were involved in WNT pathway, which is related to cashmere wool growth regulation. On the network basis, we developed a knockout mouse model to identify the network interaction. We observed that 8 high-sulfur protein genes, 12 keratin (KRT) genes, and 19 keratin associated protein (KRTAP) genes related to the growth of cashmere wool were almost not expressed in Trib3−/− rat skin. Our results suggested that the expression of genes related to wool and cashmere growth may be regulated by the interaction network between genes affected by melatonin and immune-related genes. In summary, we outlined some particularly promising ways for future research on immune-related genes of cashmere goats and the role of melatonin in wool and cashmere growth.
Introduction
The wool and cashmere growth is regulated by several mechanisms. Overall, the hair follicle (HF) is one of the most important factors of wool and cashmere growth (Jiao et al., 2019; Pazzaglia et al., 2019; He et al., 2020; Nocelli et al., 2020). Biological implications of the molecular regulatory mechanisms underlying the periodic development of HFs are important for the future studies on cashmere wool growth. Decades ago, it was reported that melatonin and immune are involved in regulating hair growth and the periodicity of HFs (Fischer et al., 2010; Wang and Higgins, 2020).
The correlation between the immune and HF has been confirmed in previous studies (Paus et al., 1999). The testis, brain, adrenal cortex, and fetal placenta are some well-defined immune “privileged” tissue compartments in mammals, and HF is one of them. Paus et al. studied the immune privilege function of HFs during hair bulb growth and speculated that a constitutive active protection system can protect HFs from immune injury. System studies of the regulatory T-cells in skin revealed that they play a major role in HF biology by promoting the function of the stem cells of HF (Ali et al., 2017). Interestingly, immune cell types change and number involved in hair growth and dependent on remodeling all skin layers (Castellana et al., 2014; Joost et al., 2020). With the continuous progress in immune and HF research, the immune cells in various stages of HF cycle and the immune cells in diseases caused by normal HF interference are becoming more and more complex. For example, recruitment of immune cells to the developing HF contributes to the embryogenesis and the skin immune system changes with the hair cycle (Paus et al., 1998). Dendritic epidermal T-cells are distributed in the epithelium of HFs in mice and humans (Hashizume et al., 1994; Paus et al., 1994; Nagao et al., 2012), which act closely with basal keratinocyte stem cells and promote wound healing via fibroblast growth factor (Jameson et al., 2004).
A study reported that melatonin, a major neuroendocrine regulator, is crucial for the function of fur phenotype and photoperiod dependent environment. The study further addressed various aspects, including melatonin changing the yield of wool and cashmere, fur development and cycle frequency, seasonal molting, and fur color (Paterson et al., 1994; Kobayashi et al., 2005). Guerrero reported that the production of melatonin is related to the circadian rhythm and seasonal changes of the immune system, the impact of melatonin injection on the immune system, and the existence of melatonin receptors in the immune system (Guerrero and Reiter, 2002). However, the role of melatonin and immunity in HF biology is still not completely clear. Because of the complexity of molecular interaction and metabolism, the exact function of melatonin and immunity in wool and cashmere biology is still unclear.
In this study, we proposed a hypothesis that the regulation of villus growth is regulated by the interaction between immune-related genes and melatonin. We systematically analyzed the immune-related genes of cashmere goats to identify the important genes and potential factors in the regulatory network of villus growth. The present study revealed a complex immune landscape consisting of multiple regulation by melatonin and immune-related genes in the blood of cashmere goat.
Materials and methods
Ethics approval and consent to participate
Consent for publication
All samples were collected in accordance with the International Guiding Principles for Biomedical Research Involving Animals and approved by the Special Committee on Scientific Research and Academic Ethics of Inner Mongolia Agricultural University, responsible for the approval of biomedical research ethics of Inner Mongolia Agricultural University (Approval No. [2020] 056). We confirm that all methods were performed in accordance with relevant guidelines and regulations and the experimental procedure complied with the Animal Research Reporting of In Vivo Experiments (ARRIVE) guidelines.
Experimental animals and sample collection
The data were obtained from the Genome Sequence Archive in BIG Data Center, Beijing Institute of Genomics (BIG), Chinese Academy of Sciences, including the RNA-seq data of 72 skin samples and 72 blood samples from 6 16-month-old female cashmere goats over 12 months. The accession number is CRA004598, which can be accessed from https://ngdc.cncb.ac.cn/gsa. Six cashmere goats were divided into two groups, including melatonin-implanted group (treatment group) and melatonin-nonimplanted group (control group), with 3 cashmere goats in each group (Chai et al., 2021). During the experiment, 2 mg/kg body weight melatonin was implanted subcutaneously behind the ears of the treatment group goats, whereas no melatonin was implanted in the control group. One month after implantation, scapula skin samples (size, 1 cm2) and venous blood samples were collected. Samples were quickly frozen in liquid nitrogen and were stored in a −80°C freezer (Chai et al., 2021).
The Trib3 knockout (KO) mice were established by Dr. Taodi Liu. The adult male homozygous (−/−), heterozygous (+/−), and wild-type (+/+) Wistar rats were selected. Detailed genetic information of Trib3 KO mice can be found in a previous article (Wang, 2020). The skin (size, 1 cm2), brain, ductus deferens, epididymics, kidney, liver, testis, and spleen tissues were collected, all 24 samples were quickly frozen in liquid nitrogen, and kept in a -80°C freezer.
RNA-seq library construction and sequencing
For the RNA-seq library, total RNA was extracted from all rat tissue samples using TRIzol Reagent as per the manufacturer’s instructions. After the RNA samples were quantified, mRNAs were enriched with oligo (DT) magnetic beads. Further, the fragment buffer was added to break the mRNA into shorter fragments. Using mRNA as the template, single-stranded cDNA was synthesized using a randomized hexamer, followed by the synthesis of double-stranded cDNA by adding buffer, dNTPs, DNA polymerase I, and RNase H. Further, the double-stranded cDNA was purified using AMPure XP beads, and the second strand of cDNA containing U was degraded using USER enzyme. Purified double-stranded cDNA was repaired, followed by the addition of a tail and connection to the sequencing connector. Next, AMPure XP beads were used to select the fragment size. Finally, the final library was obtained via PCR amplification and purification. After RNA-seq library construction, Qubit 2.0 was used for preliminary quantification, followed by library insert size detection using Agilent 2,100. After the inserts met the expectations, qPCR was used to accurately quantify the effective concentration of the cDNA library and ensure its quality. In total, 24 libraries were constructed from the samples as the final sequencing library. High-throughput sequencing was performed using the illumina novaseq 6,000 platform with a paired-end sequencing.
Quantification of the transcript levels and differentially expressed gene analysis
High quality; clean data of rats obtained by filtering data according to strict standards can be used for further research and publication (Supplementary Table S1). When the N content in any sequencing read is >10% of the base number of the read or when the low quality base number (Q ≤ 5) base present in any sequencing read was >50% of the base number of the read, the paired reads were removed. Further, we downloaded the rat genome (Rnor_6.0) and annotation files from the ensemble database. Hisat2 was used to align clean reads for each sample against the rat reference genome. StringTie (Pertea et al., 2015) program was used to quantify the transcript levels in rats.
We used Hisat2 (Kim et al., 2015) to match the read segments to the reference genome, which was obtained from the goat reference genome and its corresponding genome annotation file in the NCBI database (ARS1) (Dong et al., 2013). The differentially expressed genes (DEGs) were detected using limma in the R package (Ritchie et al., 2015) with default parameters to screen DEGs according to FPKM. The genes with |log2fc| > 1 and p < 0.05 were identified as DEGs.
Construction of consensus weighted gene coexpression network and GCNs
The TFs with FPKM>1 that were expressed in at least two samples were kept, 686 TFs common in the control and treatment groups were retained. The expression data from 36 control and 36 treatment samples underwent consensus module analysis (Langfelder and Horvath, 2008) using the weight gene coexpression network analysis (WGCNA) package for R. The first step was to choose the soft thresholding power to which coexpression similarity was raised to calculate adjacency. Soft thresholding power was set to 12 and minimum module size to 30. Modules underwent hierarchical clustering, and similar modules were merged at a cut height of 0.25. The expressed TFs were identified from 4 modules that were highly conserved between the two datasets. Module eigengene was calculated for each module as the first principal component of gene expression for that module. Correlation analysis was performed to relate module eigengene to external traits including 71 immune-related DEGs.
The regulatory network between the 12,444 expressed genes (FPKM > 1 genes in at least 2 of the 72 samples) and 686 TFs was subjected to the gene coexpression network (GCN) (Chang et al., 2019) analysis for the identification and classification of gene coexpression network. FPKMs of all genes including TFs were combined and used for the calculation of the Pearson’s correlation coefficient (PCC) to predict potential gene regulatory networks. First, the PCC of the FPKM values was calculated for all pairs of 686 TF genes and 12,444 genes under the control and treatment groups, respectively. We obtained the positive and negative cutoff values with p < 0.05 for each condition in the Cutoff program (Supplementary Table S2). In the second step, we construct eight coexpression types of GCN under two conditions (C1 and C2): C1+C2+, C1+C20, C1+C2–, C10C2+ C1–C2+, C1–C2–, C1–C20, and C10C2–, where +, -, 0 represents the positive, negative, and no coexpression, respectively. Such as GCN in January, we defined PCC 0.85 and p < 0.05 as the criteria for the coexpression for TF and non-TF genes in control group. Therefore, the genes were defined as positively coexpressed (denoted as C1+ according to the dataset) if PCC was ≥0.85. Furthermore, the genes were identified asnoncoexpressed (denoted as C10) if −0.85 < PCC <0.85. If PCC was ≤-0.85, the genes were identified as negatively coexpressed (denoted as C1−). Considering the coexpression status under C1 and C2 together, if both genes were positively coexpressed under C1 and C2, they belonged to the set of C1+ C2+. Similarly, if they were positively coexpressed only under C1 (or C2), they belonged to C1 + C20 (or C10C2 +) and were not positively coexpressed under C2 (or C1).
Weight gene coexpression network and expression-based genome-wide association studies analyses
In the blood control group, 5,497 genes (FPKM > 5 in all samples) were highly expressed. Therefore, a coexpression network was constructed using the WGCNA. The parameter β was given a default value from 1 to 30, as determined by the function SFT $powerestimate, with the minimum number of modules chosen as 30. The expression pattern of the module genes in each sample was represented by the module eigenvalue (ME). The MEs were calculated for the same expression pattern in the treatment group.
Further, the MEs of the control and treatment groups as the trait input data and eGWAS was conducted. The BAM file was used to detect single nucleotide polymorphisms (SNPs) using the samtools. SNPs were filtered using the following parameters: BIg -q 20 -d 1,000. The Binary Call Format (BCF) output file was redirected to the bcftools program to convert it to Variant Call Format (VCF) format. A total of 72 blood samples were genotyped for 7,919 (DP ≥ 10 in all samples) SNPs. To detect expression-associated eSNPs, the MEs of the control and treatment groups as the trait and expression-based genome-wide association studies (eGWAS) was conducted with 7,919 genome-wide SNPs.
Enrichment analysis and biomarker genes selection using random forest
For the TFs and expressed genes at each month of a GCN, a corresponding set of coexpressed genes can be identified. The Gene Set Enrichment Analysis (GSEA) (Subramanian et al., 2005) was conducted with the background set of all genes in every month GCN. Only gene-sets that passed conservative significance thresholds (FDR < 25%) were selected for display in the Enrichment Map (Daniele et al., 2010). We presented a random forest-based classification model to predict a biomarker gene from immune-related genes (Statistics and Breiman, 2001). A total of 645 immune-related genes (n = 72) were used as input. For the random forest-based modelling, the discovery cohort was randomly split into a training (75%) and test (25%) sets. Our mRNA-based classifier demonstrated excellent performance in discriminating the control and treatment groups in the test. The area under the receiver operating characteristic curve (AUC) for the sensitivity analysis was 89% (Supplementary Figure S1).
Results
Transcriptome-wide identification, characterization, and expression profile of immune-related genes/pathways and transcription factors in cashmere goats
To comprehensively examine the role of expression of immune-related genes during the HF development in the skin and blood tissues in the control and treatment groups, we collected the skin and blood samples at 12 continuous time points (from January to December). In total, 23,108 and 23,050 genes were identified in the skin of the control and treatment groups, respectively. Moreover, 21,378 and 21,304 genes were identified in the blood of the control and treatment groups, respectively. We defined a gene as expressed if FPKM was > 1 gene in at least 2 of the 72 samples. This resulted in 848 and 686 TFs being expressed in the skin and blood transcriptomes, respectively. Next, we searched all immune-related genes in cashmere goat genome from all annotated immune metabolic pathways in GO database (http://geneontology.org/). We found 3,308 immune related genes in 157 species covering 16,149 immune metabolic pathways, including 1,599 immune-related genes in cashmere goat blood. Furthermore, 1756 immune-related genes were found in cashmere goat skin samples. We observed that the immune-related genes expressed variably in the skin and blood over 1 year. The expression of immune-related genes in the skin demonstrated upward trends twice, from January to April and August to December, and a fluctuating trend was observed from May to July (Figure 1A). The trends of change in expression were not consistent in blood samples, with only January to July exhibiting a similar decreasing trend and August to December exhibiting an increasing trend (Figure 1B). The changes in gene expression in the skin and blood tissues associated with the treatment and control groups were analyzed. Overall, 680 DEGs were obtained in the skin (|log2fc| > 1, p < 0.05), with 366 and 314 genes upregulated in the control and treatment groups (Figure 2A). The number of upregulated DEGs in the blood was 160 and 247 in the control and treatment groups, respectively (Figure 2B). To determine the immune-related DEGs in the blood, 1,599 immune-related DEGs were analyzed. After screening, we obtained 71 immune-related DEGs in the blood of the control and treatment groups. A total of 18 genes were identified to be related to immune-related DEGs in the skin of the control and treatment groups. These cumulative changes in DEGs in the control and treatment groups between the blood and skin tissues may be minimal at the level of all genes. However, the concerned changes became evident as differential expression of two groups of immune-related genes in blood and skin. Our previous study reported that tissue-specific module genes involved in 687 pathways, such as adaptive immune response, innate immune response, immune system development, and regulation of leukocyte migration, were upregulated in the blood compared with skin (Chai et al., 2021). The functions involved immune-related GO terms were grouped for analysis in the blood. A total of 272 sets of immune-related GO terms were observed to be associated with 1,599 immune-related genes in the blood. Assessment of the expression of immune-related GO terms revealed 57 differentially expressed GO terms between the control and treatment groups in the blood of cashmere goats. In addition, 28 of the 47 differentially expressed immune-related GO terms were highly expressed in the treatment group (Figure 3).
FIGURE 1. Abundance of immune-related genes in each month. (A) FPKM of immune-related genes in skin. (B) FPKM of simmune-related genes in blood. Abscissa represents month, and ordinate represents ratio to max. D represents the control group.
FIGURE 2. Identification of differentially expressed genes. (A) Volcano plot significantly differentially expressed genes in skin. (B) Volcano plot shows the distribution of differentially expressed genes in blood.
FIGURE 3. Identification of differentially expressed immune-related pathway. The red node indicates that the pathway is up-regulated in the treatment group of the blood, The green node indicates that the pathway is up-regulated in the control group of the blood.
Identification of coexpression relationships and melatonin-dependent networks between the transcription factors and immune-related genes in cashmere goat
Consensus module provides us with a method to assess whether gene modules co expressed under one condition are also strongly coexpressed under another condition. (Langfelder and Horvath, 2007; Langfelder et al., 2011). We performed a comprehensive analysis of coexpression relationships between the TFs and immune-related genes in the developmental stages across 1 year after melatonin implantation in blood. The TFs with FPKM > 1 that were expressed in at least two samples were kept, 686 TFs of common between the control and treatment groups were retained for the WGCNA. Coexpression networks were constructed on the same expression trends across the control and treat groups. The Consensus module is defined as a highly interconnected TFs cluster, and TFs in the same cluster have a high correlation coefficient. In this analysis, we identified 4 distinct consensus modules (labeled with different colors) shown in the dendrogram (Supplementary Figure S2). The 3 consensus modules (excluding grey module; the grey module contained unassigned TFs) correlated with 71 immune-related DEGs in the control and treatment groups (Figure 4). The brown module identified 118 TFs specific to the 10 immune-related genes in the control group and 6 immune-related genes in the treatment group. Gene coexpression turquoise modules were significantly associated with 11 immune-related genes in the control group and 8 immune-related genes in the treatment group (
FIGURE 4. Module–differentially expressed immune-related genes association. Each row corresponds to a module, labeled with a color. Each column corresponds to differentially expressed immune-related genes. The color of each cell at the row-column intersection indicates the correlation coeffcient between the module and the tissue type.
Table 1). This result prompted us to further identify whether the alteration in gene regulation networks is driven by melatonin and alteration of network is affected by the change in the regulatory relationship between TFs and genes. We used a comparative, GCN method to analyze the relationship between 686 TFs and 12,444 expressed genes (FPKM > 1 gene in at least 2 of the 72 samples) in blood. A total of 8 GCNs between the control and treatment groups were obtained. It is worth noting that for the C10C2- set, the coexpression relationships were dependent on melatonin (Supplementary Table S3) and independent of the months of May and December. The coexpression network between TFs and expressed genes was based on a PCC cutoff of 0.85 and p < 0.05. By identifying the functional categories among the coexpression genes and TFs at each month of the melatonin-dependent GCN, we reported 12 immune-related function sets containing 645 immune-related genes (Figure 5) using GSEA and Enrichment Map analysis. A total of 12 immune-related function sets were upregulated throughout the year except in May and December in the treatment group. A random forest machine learning approach was used for the analysis. In total, 30 immune-related genes were identified as the significant indicators for the control and treatment group (Figure 6). Based on random forest machine learning approach, the 30 immune-related genes assemblages that were indicative of high and low expressions were present throughout the sampling period. The 4 genes (immune-related genes such as RAB7B, CFD, SLC11A1, and EMP2) were differentially expressed in the treatment and control cashmere goats (|log2fc| > 1, p < 0.05). We explored and characterized the coexpression relationships and melatonin-dependent networks between the 249 TFs and 4 immune-related genes from the C10C2- set (Figure 7C).
TABLE 1. Module−trait relationships in control and treat.
FIGURE 5. Global network of enrichment results of C10C2+ genes. The red indicate that the pathway is upregulated in the whole year except for May and December.
FIGURE 6. Expression pattern of 30 immune-related genes in blood. D is the control group, M is the implantation group, and the number is the month.
FIGURE 7. Co-expression regulatory network of high-sulfur protein gene, transcription factor, melatonin-dependent network and TIRB3-gene network. The red diamond node represents transcription factor, the pink node represents transcription factor cofactor. Gray node represents high sulfur protein gene in network (A). The gray triangular node represents biomarker gene analysis by random forest machine learning approach in network (B). The red arrow node represents biomarker gene of immune-related gene.
Integration of gene coexpression networks and expression-based genome-wide association studies analyses highlighted the candidate regulators implicated in melatonin-dependent immune gene interaction networks
To identify the regulatory mechanisms underlying the induction of the melatonin-dependent networks by immune-related genes, we performed eGWAS and WGCNA in blood samples. First, we used WGCNA for describing the correlation patterns among 5,497 genes (FPKM > 5 in all samples) in the control group. Using WGCNA, we identified 6 gene coexpression modules in the control group (Supplementary Figure S3). The number of genes in each module ranged from 6 (module grey) to 1,485 (module turquoise). Assessment of MEs for the 6 gene modules in the treatment group revealed that a great variation existed in the MEs compared with the control group. Further, the MEs of the control and treatment groups as the WGCNA trait input data and eGWAS was conducted with 7,919 genome-wide SNPs. Using eGWAS, the identified 25 eSNPs were mapped in intergenic regions andassociated with 20 genes in the control and treatment groups (Supplementary Table S4). Remarkably, we observed a significant genetic variation in the treatment group, and brown module associated genes are the most. Exploratory analyses were conducted with brown module genes and 9 biomarker genes. The assessment of regulatory network between brown module genes and biomarker genes was conducted based on the PCC in the treatment group. Regulatory gene networks (cor ≥ 0.8, p ≤ 0.05) were identified involving 991 genes and 8 biomarker genes, such as LOC102189514, LOC106501899, LOC102185917, MYADM, LOC108633278, Trib3, TIRAP, and SUDS3 (Figure 7B). The genes in the interaction network between brown module genes and biomarker genes seem to be closely related to RNA transport, RNA degradation, spliceosome, protein export, protein processing in the endoplasmic reticulum KEGG pathway, and gene enrichment in the regulatory network (Supplementary Figure S4). Because the expression of many genes in the genome is regulated by TFs, we observed that the network included 24 TFs for target gene sets in the treatment group, and these 24 TFs simultaneously regulated the melatonin-dependent networks between the TFs and immune-related genes in cashmere goats. We previously reported that high-sulfur gene expression is important for wool growth regulation (Chai et al., 2021). A total of 47 high-sulfur protein genes exhibited an interaction with TFs and cofactors with ASE in the skin of cashmere goat. These TFs and cofactors were inhibited after melatonin implantation. Our data are particularly alarming because overlapping TFs, such as PHF5A, REXO4, STRAP, JUNB, GATAD2A, ZNF710, and VDR, were already observed in 3 networks (Figure 7C), and these TFs were also expressed in blood and skin tissues of cashmere goat. This indicated that many biological functions cannot be explained by the function of a single gene or protein. Instead, they are the result of a network of interactions between two genes or between proteins and other molecules.
Role of Trib3 in the regulation of the expression of high-sulfur protein genes in the rat skin
The identification and understanding of the gene regulatory network of wool growth regulation process after melatonin implantation is one of the main challenges. Another challenge is to identify the main driving factors and major regulatory genes that control the wool growth regulation. In this study, we observed that 3 major regulatory networks were interconnected through TFs, and some TFs, such as PHF5A, REXO4, STRAP, JUNB, GATAD2A, ZNF710, and VDR, were also expressed in the blood and skin tissue of cashmere goat. In addition, most genes in the network, such as VDR, JUNB, and Trib3, were involved in WNT pathway, which is related to wool growth regulation. The close cooperation between 3 networks gives rise to an interesting question: how does regulatory network collaborate with genes for controlling gene expression in skin and blood. To address this question, the network of Trib3 gene in blood was selected as candidate genes because of the involvement of this gene in the WNT pathway. In the network of Trib3 genes, the TFs, namely, PHF5A, REXO4, STRAP, and JUNB interact with the regulatory network of Trib3 in blood and with high-sulfur protein genes in skin. Based on this, we guessed that the regulatory network of high-sulfur protein gene and Trib3 gene also interacted in the skin. However, this is just a conjecture. Further, we investigated the impact of Trib3 gene knockout on the high-sulfur protein genes expression in skin, brain, ductus deferens, epididymics, kidney, liver, testis, and spleen tissues in rat. Using homologous gene sequences of high-sulfur protein genes in cashmere goat as our query terms, we searched for 8 candidate genes from the rat using BLAST. We used the “absolute” standard for filtering, which simply removes BLAST hits that do not meet a top1 criterion; various possible criteria include bit score, E-value, and alignment length. The RNA-seq analysis revealed specific expression of 8 homologous genes (Supplementary Table S5) of high-sulfur protein genes in rat skin tissue (Supplementary Table S6). We investigated the expression of high-sulfur protein genes in more detail using IGV (Figure 8). The 8 high-sulfur protein genes, 12 KRT genes, and 19 KRTAP genes were not expressed in adult male homozygous (−/−) skin tissue. The expression of high-sulfur protein genes in the wild-type (+/+) Wistar rat skin was significantly increased. However, we observed that the decreased expression became conspicuous in adult male heterozygous (+/−) rat (Supplementary Table S7).
FIGURE 8. The high sulfur protein genes expression in more detail by IGV. The high sulfur protein gene were not expressed in adult male homozygous (−/−) skin tissue, the expression of high sulfur protein genes in the wild-type (+/+) Wistar rat skin was marked blue. The green node represent the expression of high sulfur protein genes in heterozygous (+/−) rat.
Discussion
It is increasingly clear that the immune and melatonin might contribute to HF development. However, comprehensive analysis of the differentially expressed profiles of immune-related genes after implantation of melatonin and potential interactions between TFs and immune-related genes in blood and skin through integrating RNA-seq and construction of a rat gene knockout model are not available.
Researchers have confirmed that TFs, such as Lhx2 and csdc2, regulate the secondary HF cycle of cashmere goats, which are potentially important TFs in the hair growth cycle. In the present study, we obtained complete datasets of immune-related genes and TFs in cashmere goat after analyzing datasets from GO and Animal TFDB database. In addition, to further understand the molecular mechanism of immunity in the regulation of wool growth in cashmere goats, we analyzed the expression patterns of immune-related genes in blood and skin of cashmere goats.
Focusing on interaction between the TFs and immune-related genes, we identified melatonin-dependent gene regulatory networks that might serve as a basis for investigating mechanisms related to wool and cashmere growth function. Using GCNs and random forest machine learning approach, we identified 4 biomarker immune-related genes that were differentially expressed in the melatonin-implanted and control cashmere goats (|log2fc| > 1, p < 0.05), including RAB7B, CFD, SLC11A1, and EMP2. A recent study regarding the influence of the senescent cells on matrix production reported that the secretion of CFD promotes degradation of nearby fibroblasts in the dermal layer (Ezure et al., 2019). Evidence supports that SLC11A1, encoding a transport protein, is expressed in macrophage phagosome membranes, contributing to the ion transport (Vidal et al., 1993; Govoni et al., 1997; Canonne-Hergaux et al., 1999). Codon 169 of SLC11A1 gene encodes Gly in the R allele and ASP in the S allele. The two alleles are related to the sensitivity and drug resistance in the enteric typhoid infection by Salmonella (Malo and Skamene, 1994; Cellier et al., 1995; Vidal et al., 1995; Cellier et al., 1996; Cellier et al., 1997). The mutation of S allele can change the ion concentrations in the phagosome (Hackam et al., 1998; Zaharik et al., 2004). SLC11A1 gene is critical to regulating macrophage activation, oxidative and nitrosamine bursts (Fritsche et al., 2003), and expression of TNF-a, IFN-c, IL-1 (Lalmanach et al., 2001), and MHC II (Wojciechowski et al., 1999). RAB7B is a paralogous gene of RAB7A and is involved in membrane tracking in the late endosomal/lysosomal pathway (Osińska et al., 2011). A study reported that EMP2 may also be involved in cell proliferation and cell–cell interactions (Wadehra et al., 2003). In summary, although there is overwhelming information demonstrating the relationship between the immune and melatonin, many questions related to TFs and the mechanisms of other aspects are unanswered. We suspect that many additional mechanisms underlying the formation and maintenance of HFs and growth of wool and cashmere have not yet been discovered. We only see the “tip of the iceberg.”
Although SNPs found in RNA sequence data have limitations, the results are encouraging. Remarkably, we observed a significant genetic variation in the treatment group. The regulatory network between 991 genes and 8 biomarker genes, such as LOC10218951, LOC10650189, LOC102185917, MYADM, LOC108633278, Trib3, TIRAP, and SUDS3, was analyzed based on the PCC in the treatment group (Figure 7B). We observed that this network includes 24 TFs for target gene sets in the treatment group, and these 24 TFs simultaneously regulated the melatonin-dependent networks (Figure 7C) between the TFs and immune-related genes in cashmere goats. In our previous study, 47 high-sulfur protein genes exhibited interaction with TFs and cofactors with ASE in the skin of cashmere goat (Chai et al., 2021). These TFs and cofactors were inhibited after melatonin implantation in the skin of cashmere goat, and 13 genes in the interaction network were functionally associated with the function of hair growth and WNT pathway (Chai et al., 2021). Our data are particularly alarming because overlapping TFs, such as PHF5A, REXO4, STRAP, JUNB, GATAD2A, ZNF710, and VDR, were already observed in 3 networks (Figure 7) and were also expressed in blood and skin tissue of cashmere goat. In the study by Wang et al., VDR and JUNB gene specific expression at embryonic day 120 are thought to play an important role in HF and WNT pathway (Wang et al., 2020). The transcriptional activity of β-catenin-TCF4 was determined by Trib3 and β-catenin synergism, which then affects the initiation and progression of CRC (Hua et al., 2018). Therefore, Trib3 was selected as the candidate gene. The major types of networks between the immune-related genes, high-sulfur protein genes, and TFs in blood and skin after melatonin implantation were constructed. The network contained 7 bridge genes in the cashmere goat of blood and skin. To investigate whether the network gene in blood affects the high-sulfur protein genes in skin, we analyzed the interaction between Trib3 and high-sulfur protein in the skin (size, 1 cm2), brain, ductus deferens, epididymis, kidney, liver, testis, and spleen tissues of rats with Trib3 gene knockout. Notably, 8 high-sulfur protein, KRT, and KRTAP genes were specifically expressed (Supplementary Table S6) in the rat skin tissue. These genes were not expressed in adult male homozygous (−/−) skin tissue, and the expression of high-sulfur protein genes in the wild-type (+/+) Wistar rat skin was significantly increased. However, we observed that the decreased expression became conspicuous in adult male heterozygous (+/−) rats.
Mature wool mainly consists of three layers: cortical, cuticle, and wool central core and medulla (Powell and Rogers, 1997a). Studies have reported that cortical keratin consists of microfibers (keratin intermediate filaments) of 8-nm diameter embedded in the matrix. Keratin-related proteins were located the matrix between microfibers (Plowman, 2003). A classification system was proposed based on whether the wool fiber contains keratin intermediate filaments (KIF). KRTm. nxpL is used to classify intermediate filaments, whereas other proteins with high-sulfur, ultra-high-sulfur and high glycine–tyrosine are used as keratin-associated proteins (abbreviated as KAPm.nxpL) (Rogers and Powell, 1993; Powell and Rogers, 1997b). Overall, 90% of KIF and KAP constitute the main body of wool fibers, and the spatial structure of KIFs and chemical bond of KAP in the matrix are believed to largely determine the physical properties of fibers (Gong et al., 2016). While discussing these aspects, the importance of their expression for cashmere fiber and length traits was emphasized, that is, through the current understanding of the mechanism by which keratin and keratin-associated proteins affect wool traits and the expression of specific genes as part of these mechanisms.
Compared with other species, in goats, keratin-associated protein and keratin can also be classified into its gene families based on the similarity of their amino acid sequences. Sequence analyses suggest that the KAP and KRT sequences identified to date probably represent different family members. The KAP1 family of goat fiber proteins, such as KAP1.1 and KAP1.4, are primarily described in recent years, which is a HS-KAP family (Zhang et al., 2011; Shah et al., 2013). Liu et al. used PCR-SSCP method to reveal that the KAP8.2 gene may be a potential molecular marker for determining the cashmere fiber diameter of cashmere goats, and the genotype polymorphism has no effect on the length of cashmere (Liu et al., 2007). KAP8.2 is strongly expressed in the cortical tissue of primary and secondary HFs of embryonic and adult goat skin with high homology at the nucleotide and amino acid levels in goats of different breeds and from different regions (Hang et al., 2008). The study on KAP6-1.2 in goats indicated that this gene was strongly located on the cortical layers of primary and secondary HFs in the skin of embryonic and adult goats (Zhang et al., 2007). The study by JIN et al. reported that KAP7.1 genes are specifically expressed in the HF, which can effectively regulate the HF growth and development in cashmere goats and may even affect the wool fiber diameter; the expression of secondary follicles was significantly higher in primary HFs (Jin et al., 2011). KAP8.2 is strongly expressed in the cortical tissue of primary and secondary HFs of embryonic and adult goat skin with high homology at the nucleotide and amino acid levels in goats of different breeds and from different regions (Zhang et al., 2008). KAP24-1 is one of the most recently discovered high-sulfur KAP gene in goats. Wang reported that this KAP gene domain was on goat chromosome 1 (Wang et al., 2019). Prior to the publication of Wang, KAP13.1 gene, were reported (Fang et al., 2005). The discovery of the polymorphisms of goat KAP13.1 genes has suggested that polymorphic variants might be relevant in determining fiber diameter. Not fully unexpectedly, the patterns of KAP11-1 expression in the caprine skin, to a certain degree, that of fiber diameter size. Real-time PCR results revealed that KAP11.1 exhibited higher expression in the catagen than in the anagen phase in the primary HFs. However, in the secondary HFs, KAP11.1 exhibited a significantly higher expression in the anagen than in the catagen phase. Moreover, KAP11.1 gene exhibited a strong expression in the inner root sheath and hair matrix and a lower expression in the hair bulb (Jin et al., 2017). A search for similar sequences in the caprine genome using the human KAP20-2 gene (KRTAP20-2) revealed a homologous sequence on chromosome 1, which contains an open reading frame of 189 bp, and the nucleotide homology with human KAP20-2 is 75% (Wang et al., 2017).
Even in the likely event that the new focus in immune-related gene research on the still insufficiently explored potential regulated genes of the network that we strongly advocate here. However, this focus is guaranteed to produce a rich immunological harvest in terms of improving our general understanding of the comprehensive regulation of wool and cashmere growth associated with melatonin, immune-related genes, high-sulfur protein genes, ASE, and SNPs.
Data availability statement
The datasets analyzed for this study can be found in the Genome Sequence Archive in BIG Data Center, Beijing Institute of Genomics (BIG), Chinese Academy of Sciences, under accession number CRA007385 and accession number CRA004598 that are publicly accessible at https://ngdc.cncb.ac.cn/gsa.
Ethics statement
All samples were collected in accordance with the International Guiding Principles for Biomedical Research Involving Animals. The animal study was reveiewed and approved by the Special Committee on Scientific Research and Academic Ethics of Inner Mongolia Agricultural University, responsible for the approval of biomedical research ethics of Inner Mongolia Agricultural University (Approval No. [2020] 056). We confirm that all methods were carried out in accordance with relevant guidelines and regulations. We confirm that the authors complied with the Animal Research Reporting of In Vivo Experiments (ARRIVE) guidelines.
Author contributions
We thank study participants for their contribution. YC: writing-original draft preparation and writing-reviewing and editing. WZ: supervision and project administration. SF, BL, LG, ZL, LD, and YS: experimental sample and data analysis. WZ: Writing-reviewing and editing, CL and TL: funding acqusision.
Funding
The reported work was supported by the National Natural Science Foundation of China (no.31672385, no.31402052), Inner Mongolia Natural Science foundation (no. 2019MS03078), and from the Special Open Project on the Construction of Mongolia Pharmacy Doctoral Station for National Service with Special Needs (no. MDMYBJ2019004) and the Doctoral Scientific Research Foundation of Inner Mongolia University for Nationalities (no.BS308). This work is also supported by the Doctoral Scientific Research Foundation of Inner Mongolia Medical University (no. YKD2019BSJJ009). The funding played a role in the design of the study and the collection, analysis, and interpretation of data.
Conflict of interest
BL, was empolyed by Nei Mongol BioNew Technology Co. Ltd.
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.
Publisher’s note
All claims expressed in this article are solely those of the authors and do not necessarily represent those of their affiliated organizations, or those of the publisher, the editors and the reviewers. Any product that may be evaluated in this article, or claim that may be made by its manufacturer, is not guaranteed or endorsed by the publisher.
Supplementary material
The Supplementary Material for this article can be found online at: https://www.frontiersin.org/articles/10.3389/fgene.2022.967402/full#supplementary-material
Supplementary Figure S1 | ROC curve discrimination between the control and treatment group.
Supplementary Figure S2 | Consensus weighted gene coexporession network analysis. Modules were identified in the resulting dendrogram using the Dynamic Tree Cut algorithm. A total of 4 modules were identified. Modules were distinguished from each other by assigning different colors.
Supplementary Figure S3 | Network heatmap plot, all genes.
Supplementary Figure S4 | “Genes of the regulatory networks target-biological process-signaling pathway” network.
Supplementary Table S1 | Summary of RNA-sequencing results and quality data output.
Supplementary Table S2 | The positive and negative cutoff values for the control group and the melatonin-implanted group.
Supplementary Table S3 | Numbers of co-expression relationships in each developmental stage and GCN.
Supplementary Table S4 | The module-associated SNP by eGWAS analysis.
Supplementary Table S5 | The blast search results.
Supplementary Table S6 | The expression of 8 high sulfur protein gene in rat tissues.
Supplementary Table S7 | The expression of 12 KRT gene and 19 KRTAP gene in rat skin tissue.
Abbreviations
HF, hair follicle; DEGs, differentially expressed genes; TFs, transcription factor; WGCNA, weight gene co-expression network analysis; GCN, gene co-expression network; PCC, Pearson correlation coefficient; ME, module eigenvalue; SNP, single nucleotide polymorphisms; eGWAS, Expression-based genome-wide association studie; AUC, receiver operating characteristic curve.
References
Ali, N., Zirak, B., Rodriguez, R., S., Pauli, M. L., Truong, H. A., Lai, K., et al. (2017). Regulatory T cells in skin facilitate epithelial stem cell differentiation. Cell 169 (6), 1119–1129. e11. doi:10.1016/j.cell.2017.05.002
Canonne-Hergaux, F., Gruenheid, S., Govoni, G., and Gros, P. (1999). The Nramp1 protein and its role in resistance to infection and macrophage function. Proc. Assoc. Am. Physicians 111 (4), 283–289. doi:10.1046/j.1525-1381.1999.99236.x
Castellana, D., Paus, R., and Perez-Moreno, M. (2014). Macrophages contribute to the cyclic activation of adult hair follicle stem cells. PLoS Biol. 12 (12), e1002002–e1002016. doi:10.1371/journal.pbio.1002002
Cellier, M., Belouchi, A., and Gros, P. (1996). Resistance to intracellular infections: Comparative genomic analysis of nramp. Trends Genet. 12 (6), 201–204. doi:10.1016/0168-9525(96)30042-5
Cellier, M., Prive, G., Belouchi, A., Kwan, T., Rodrigues, V., Chia, W., et al. (1995). Nramp defines a family of membrane proteins. Proc. Natl. Acad. Sci. U. S. A. 92 (22), 10089–10093. doi:10.1073/pnas.92.22.10089
Cellier, M., Shustik, C., Dalton, W., Rich, E., Hu, J., Malo, D., et al. (1997). Expression of the human NRAMP1 gene in professional primary phagocytes: Studies in blood cells and in HL-60 promyelocytic leukemia. J. Leukoc. Biol. 61 (1), 96–105. doi:10.1002/jlb.61.1.96
Chai, Y., Sun, Y., Jia, C., Bin, L., Guo, L., L., Zhou, L., et al. (2021). Role of sulfur metabolism gene and high-sulfur gene expression in wool growth regulation in the cashmere goat. Front. Genet. 12, 715526. doi:10.3389/fgene.2021.715526
Chang, Y., M., Lin, H., H., Liu, W. Y., Yu, C. P., Chen, H. J., Wartini, P. P., et al. (2019). Comparative transcriptomics method to infer gene coexpression networks and its applications to maize and rice leaf transcriptomes. Proc. Natl. Acad. Sci. U. S. A. 116 (8), 3091–3099. doi:10.1073/pnas.1817621116
Daniele, M., Isserlin, R., Stueker, O., Emili, A., and Bader, G. D. (2010). Enrichment Map: A network-based method for gene-set enrichment visualization and interpretation. Plos One 5 (11), e13984. doi:10.1371/journal.pone.0013984
Dong, Y., Xie, M., Jiang, Y., Xiao, N. Q., Du, X. Y., Zhang, W. G., et al. (2013). Sequencing and automated whole-genome optical mapping of the genome of a domestic goat (Capra hircus). Nat. Biotechnol. 31 (2), 135–141. doi:10.1038/nbt.2478
Ezure, T., Sugahara, M., and Amano, S. (2019). Senescent dermal fibroblasts negatively influence fibroblast extracellular matrix-related gene expression partly via secretion of complement factor D. Biofactors 45 (4), 556–562. doi:10.1002/biof.1512
Fang, Y., Liu, W., J., Zhang, F., Q., Shao, Y., and Yu, S. (2005). The polymorphism of a novel mutation of KAP13.1 gene and its associations with cashmere traits on xinjiang local goat breed in China. Asian J. Anim. Vet. Adv. 5 (1), 34–42. doi:10.3923/ajava.2010.34.42
Fischer, T., W., Slominski, A., Tobin, D., J., and Paus, R. (2010). Melatonin and the hair follicle. J. Pineal Res. 44 (1), 1–15. doi:10.1111/j.1600-079X.2007.00512.x
Fritsche, G., Dlaska, M., Barton, H., Theurl, I., Garimorth, K., and Weiss, G. (2003). Nramp1 functionality increases inducible nitric oxide synthase transcription via stimulation of IFN regulatory factor 1 expression. J. Immunol. 171 (4), 1994–1998. doi:10.4049/jimmunol.171.4.1994
Gong, H., Zhou, H. T., Forrest, R. H. J., Li, S. B., Wang, J. Q., Dyer, J. M., et al. (2016). Wool keratin-associated protein genes in sheep—a review. Genes 7 (6), 24. doi:10.3390/genes7060024
Govoni, G., Gauthier, S., Billia, F., Iscove, N. N., and Gros, P. (1997). Cell-specific and inducible Nramp1 gene expression in mouse macrophages in vitro and in vivo. J. Leukoc. Biol. 62 (2), 277–286. doi:10.1002/jlb.62.2.277
Guerrero, J., M., and Reiter, R., J. (2002). Melatonin-immune system relationships. Curr. Top. Med. Chem. 2 (2), 167–179. doi:10.2174/1568026023394335
Hackam, D., J., Rotstein, O., D., Zhang, W., Gruenheid, S., Gros, P., and Grinstein, S. (1998). Host resistance to intracellular infection: Mutation of natural resistance-associated macrophage protein 1 (Nramp1) impairs phagosomal acidification. J. Exp. Med. 188 (2), 351–364. doi:10.1084/jem.188.2.351
Hang, J., X., Wang, L., and Yin, J., (2008). Analysis of cloning and expression for goat KAP8.2 gene. Biotechnology.
Hashizume, H., Tokura, Y., and Takigawa, M. (1994). Increased number of dendritic epidermal T cells associated with induced anagen phase of hair cycles. J. Dermatol. Sci. 8 (2), 119–124. doi:10.1016/0923-1811(94)90005-1
He, N., Su, R., Wang, Z., Zhang, Y., and Li, J. (2020). Exploring differentially expressed genes between anagen and telogen secondary hair follicle stem cells from the Cashmere goat (Capra hircus) by RNA-Seq. PLoS ONE 15 (4), e0231376. doi:10.1371/journal.pone.0231376
Hua, F., Shang, S., Yang, Y., W., Zhang, H. Z., Xu, T. L., Yu, J. J., et al. (2018). TRIB3 interacts with β-catenin and TCF4 to increase stem cell features of colorectal cancer stem cells and tumorigenesis. Gastroenterology 156 (3), 708–721. e15. doi:10.1053/j.gastro.2018.10.031
Jameson, J., M., Sharp, L., L., Witherden, D. A., and Havran, W. L. (2004). Regulation of skin cell homeostasis by gamma delta T cells. Front. Biosci. 9, 2640–2651. doi:10.2741/1423
Jiao, Q., Yin, R., H., Zhao, S., J., Wang, Z. Y., Zhu, Y. B., Wang, W., et al. (2019). Identification and molecular analysis of a lncRNA-HOTAIR transcript from secondary hair follicle of cashmere goat reveal integrated regulatory network with the expression regulated potentially by its promoter methylation. GENE 688, 182–192. doi:10.1016/j.gene.2018.11.084
Jin, M., Cao, Q., Wang, R., Piao, J., Zhao, F., and Piao, J. (2017). Molecular characterization and expression pattern of a novel Keratin-associated protein 11.1 gene in the Liaoning cashmere goat (Capra hircus). Asian-Australas. J. Anim. Sci. 30 (3), 328–337. doi:10.5713/ajas.16.0078
Jin, M., Wang, L., Li, S., Xing, M. X., and Zhang, X. (2011). Characterization and expression analysis of KAP7.1, KAP8.2 gene in Liaoning new-breeding cashmere goat hair follicle. Mol. Biol. Rep. 38 (5), 3023–3028. doi:10.1007/s11033-010-9968-6
Joost, S., Annusver, K., Jacob, T., Sun, X., Dalessandri, T., Sivan, U., et al. (2020). The molecular anatomy of mouse skin during hair growth and rest. Cell stem Cell 26 (3), 441–457. e7. doi:10.1016/j.stem.2020.01.012
Kim, D., Langmead, B., and Salzberg, S., L. (2015). Hisat: A fast spliced aligner with low memory requirements. Nat. Methods 12 (4), 357–360. doi:10.1038/nmeth.3317
Kobayashi, H., Kromminga, A., Dunlop, T., W., Tychsen, B., Conrad, F., Suzuki, N., et al. (2005). A role of melatonin in neuroectodermal-mesodermal interactions: The hair follicle synthesizes melatonin and expresses functional melatonin receptors. Faseb J. 19 (12), 1710–1712. doi:10.1096/fj.04-2293fje
Lalmanach, A., C., Montagne, A., Menanteau, P., and Lantier, F. (2001). Effect of the mouse Nramp1 genotype on the expression of IFN-gamma gene in early response to Salmonella infection. Microbes Infect. 3 (8), 639–644. doi:10.1016/s1286-4579(01)01419-8
Langfelder, P., and Horvath, S. (2007). Eigengene networks for studying the relationships between co-expression modules. BMC Syst. Biol. 1, 54. doi:10.1186/1752-0509-1-54
Langfelder, P., and Horvath, S. (2008). Wgcna: an R package for weighted correlation network analysis. Bmc Bioinforma. 9 (1), 559. doi:10.1186/1471-2105-9-559
Langfelder, P., Luo, R., Oldham, M. C., and Horvath, S. (2011). Is my network module preserved and reproducible? PLoS Comput. Biol. 7 (1), e1001057. doi:10.1371/journal.pcbi.1001057
Liu, H., Ning, L., Jia, C., Zhu, X., and Jia, Z. (2007). Effect of the polymorphisms of keratin associated protein 8.2 gene on fibre traits in inner Mongolia cashmere goats. Asian-Australas. J. Anim. Sci. 20 (6), 821–826. doi:10.5713/ajas.2007.821
Malo, D., and Skamene, E. (1994). Genetic control of host resistance to infection. Trends Genet. 10 (10), 365–371. doi:10.1016/0168-9525(94)90133-3
Nagao, K., Kobayashi, T., Moro, K., Ohyama, M., Adachi, T., Kitashima, D. Y., et al. (2012). Stress-induced production of chemokines by hair follicles regulates the trafficking of dendritic cells in skin. Nat. Immunol. 13 (8), 744–752. doi:10.1038/ni.2353
Nocelli, C., Cappelli, K., Capomaccio, S., Pascucci, L., Mercati, F., Pazzaglia, I., et al. (2020). Shedding light on cashmere goat hair follicle biology: From morphology analyses to transcriptomic landascape. BMC Genomics 21 (1), 458. doi:10.1186/s12864-020-06870-x
Osińska, M., Wiejak, J., Wypych, E., Bilski, H., Bartosiewicz, R., and Wyroba, E. (2011). Distinct expression, localization and function of two Rab7 proteins encoded by paralogous genes in a free-living model eukaryote. Acta Biochim. Pol. 58 (4), 597–607. PMID: 22030555. doi:10.18388/abp.2011_2230
Paterson, A., M., Foldes, A., et al. (1994). Melatonin and farm animals: Endogenous rhythms and exogenous applications. J. Pineal Res. 16 (4), 167–177. doi:10.1111/j.1600-079x.1994.tb00097.x
Paus, R., Christoph, T., and Muller-Rover, S. (1999). Immunology of the hair follicle: A short journey into terra incognita. J. Investig. Dermatol. Symp. Proc. 4 (3), 226–234. doi:10.1038/sj.jidsp.5640217
Paus, R., Hofmann, U., Eichmuller, S., and Czarnetzki, B. M. (1994). Distribution and changing density of gamma-delta T cells in murine skin during the induced hair cycle. Br. J. Dermatol. 130 (3), 281–289. doi:10.1111/j.1365-2133.1994.tb02922.x
Paus, R., van der Veen, C., Eichmuller, S., Kopp, T., HagEn, E., Muller-Rover, S., et al. (1998). Generation and cyclic remodeling of the hair follicle immune system in mice. J. Invest. Dermatol. 111 (1), 7–18. doi:10.1046/j.1523-1747.1998.00243.x
Pazzaglia, I., Mercati, F., Antonini, M., Capomaccio, S., Cappelli, K., Dall'Aglio, C., et al. (2019). PDGFA in cashmere goat: A motivation for the hair follicle stem cells to activate. Animals. 9 (2), 38. doi:10.3390/ani9020038
Pertea, M., E., 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 (3), 290–295. doi:10.1038/nbt.3122
Plowman, J., E. (2003). Proteomic database of wool components. J. Chromatogr. B Anal. Technol. Biomed. Life Sci. 787 (1), 63–76. doi:10.1016/s1570-0232(02)00211-8
Powell, B., C., and Rogers, G., E. (1997). The role of keratin proteins and their genes in the growth, structure and properties of hair. EXS 78, 59–148. doi:10.1007/978-3-0348-9223-0_3
Powell, B., C., and Rogers, G., E. (1997). The role of keratin proteins and their genes in the growth, structure and properties of hair. Inf. Struct. Hum. Hair 78, 59–148. doi:10.1007/978-3-0348-9223-0_3
Ritchie, M. E., Belinda, P., Wu, D., Hu, Y. F., Law, C. W., Shi, W., et al. (2015). Limma powers differential expression analyses for RNA-sequencing and microarray studies. Nucleic Acids Res. 43 (7), e47. doi:10.1093/nar/gkv007
Rogers, G., E., and Powell, B., C. (1993). Organization and expression of hair follicle genes. J. Investigative Dermatology 101, S50–S55. ∼55. doi:10.1016/0022-202x(93)90500-h
Shah, R., M., Ganai, T., Sheikh, F., D., Shanaz, S., ShabirM., , and Khan, H. M. (2013). Characterization and polymorphism of keratin associated protein 1.4 gene in goats. Gene 518 (2), 431–442. doi:10.1016/j.gene.2012.12.021
Statistics, L., B., and Breiman, L. (2001). Random forests. Mach. Learn 45 (1), 5–32. doi:10.1023/A:1010950718922
Subramanian, A., Tamayo, P., Mootha, V. K., Mukherjee, S., Ebert, B. L., Gillette, M. A., et al. (2005). Gene set enrichment analysis: A knowledge-based approach for interpreting genome-wide expression profiles. Proc. Natl. Acad. Sci. U. S. A. 102 (43), 15545–15550. doi:10.1073/pnas.0506580102
Vidal, S., M., Malo, D., Vogan, K., Skamene, E., and Gros, P. (1993). Natural resistance to infection with intracellular parasites: Isolation of a candidate for bcg. Cell 73 (3), 469–485. doi:10.1016/0092-8674(93)90135-d
Vidal, S., Tremblay, M., L., Govoni, G., Gauthier, S., Sebastiani, G., Malo, D., et al. (1995). The ity/lsh/bcg locus: Natural resistance to infection with intracellular parasites is abrogated by disruption of the Nramp1 gene. J. Exp. Med. 182 (3), 655–666. doi:10.1084/jem.182.3.655
Wadehra, M., Su, H., Gordon, K., L., Goodglick, L., and Braun, J. (2003). The tetraspan protein EMP2 increases surface expression of class I major histocompatibility complex proteins and susceptibility to CTL-mediated cell death. Clin. Immunol. 107 (2), 129–136. doi:10.1016/s1521-6616(03)00048-2
Wang, E., and Higgins, C., A. (2020). Immune cell regulation of the hair cycle. Exp. Dermatol. 29 (3), 322–333. doi:10.1111/exd.14070
Wang, J., Q., Che, L. J., Hickford, J. g. H., Hao, Z., Luo, Y., Hu, J., et al. (2017). Identification of the caprine keratin-associated protein 20-2 (KAP20-2) gene and its effect on cashmere traits. Genes 8 (11), 328. doi:10.3390/genes8110328
Wang, J. L. (2020). Effects of Trib3 gene knockout on spermatogenesis in male rats, 8. Inner Mongolia Medical University, 23. (in Chinese).
Wang, J., Zhou, H., Luo, Y., Zhao, M., Gong, H., Hao, Z., et al. (2019). Variation in the caprine KAP24-1 gene affects cashmere fibre diameter. Animals. 9 (1), E15–E12. doi:10.3390/ani9010015
Wang, S., Li, F., Liu, J., Zhang, Y., Zheng, Y., Ge, W., et al. (2020). Integrative analysis of methylome and transcriptome reveals the regulatory mechanisms of hair follicle morphogenesis in cashmere goat. Cells 9 (4), 969. doi:10.3390/cells9040969
Wojciechowski, W., DeSanctis, J., Skamene, E., and Radzioch, D. (1999). Attenuation of MHC class II expression in macrophages infected with Mycobacterium bovis bacillus Calmette–Guerin involves class II transactivator and depends on the Nramp1 gene. J. Immunol. 163 (5), 2688–2696. PMID: 10453010.
Zaharik, M., L., Cullen, V., L., Fung, A., M., Libby, S., J., Kujat, Choy., S., L., Coburn, B., et al. (2004). The Salmonella enterica serovar typhimurium divalent cation transport systems MntH and SitABCD are essential for virulence in an Nramp1G169 murine typhoid model. Infect. Immun. 72 (9), 5522–5525. doi:10.1128/IAI.72.9.5522-5525.2004
Zhang, J., X., Wang, L., Yin, J., and li, J., Q. (2007). Analysis of cloning and expression for goat KAP6 - 1. 2 gene. BIOTECHNOLOGY 17 (3), 20–23. (in Chinese).
Zhang, J., X., Wang, L., Yin, J., and li, J., Q. (2008). Analysis of cloning and expression for goat KAP8-2 gene. BIOTECHNOLOGY 18 (4), 15–18. (in Chinese).
Keywords: cashmere goat, immune genes, melatonin, wool growth, high sulfur protein genes
Citation: Chai Y, Liu Z, Fu S, Liu B, Guo L, Dai L, Sun Y, Zhang W, Li C and Liu T (2022) Effects of exogenous melatonin on expressional differences of immune-related genes in cashmere goats. Front. Genet. 13:967402. doi: 10.3389/fgene.2022.967402
Received: 12 June 2022; Accepted: 30 September 2022;
Published: 24 October 2022.
Edited by:
Ran Di, Institute of Animal Sciences, Chinese Academy of Agricultural Sciences, ChinaReviewed by:
Wenzhu Yang, Agriculture and Agri-Food Canada (AAFC), CanadaKarim Hasanpur, University of Tabriz, Iran
Copyright © 2022 Chai, Liu, Fu, Liu, Guo, Dai, Sun, Zhang, Li 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: Wenguang Zhang, atcgnmbi@aliyun.com; Chun Li, lichun1985@126.com; Taodi Liu, taodiliu@163.com