Front. Genet. Frontiers in Genetics Front. Genet. 1664-8021 Frontiers Media S.A. 1239979 10.3389/fgene.2023.1239979 Genetics Original Research Genetic diversity, tissue-specific expression, and functional analysis of the ATP7A gene in sheep Li et al. 10.3389/fgene.2023.1239979 Li Hao 1 2 Du Xiaolong 2 Li Xinyue 2 Feng Pingjie 2 Chu Mingxing 2 * Jin Yi 1 * Pan Zhangyuan 2 * 1 Engineering Research Center of North-East Cold Region Beef Cattle Science & Technology Innovation, Ministry of Education, College of Agriculture, Yanbian University, Yanji, China 2 Key Laboratory of Animal Genetics, Breeding and Reproduction of Ministry of Agriculture and Rural Affairs, Institute of Animal Science, Chinese Academy of Agricultural Sciences, Beijing, China

Edited by: Angela Cánovas, University of Guelph, Canada

Reviewed by: Poliana Fernanda Giachetto, EMBRAPA Agricultural Informatics, Brazil

Nagarjun Vijay, Indian Institute of Science Education and Research, Bhopal, India

 *Correspondence: Mingxing Chu, mxchu@263.net; Yi Jin, yijin@ybu.edu.cn; Zhangyuan Pan, zhypan01@163.com 20 09 2023 2023 14 1239979 14 06 2023 06 09 2023 Copyright © 2023 Li, Du, Li, Feng, Chu, Jin and Pan. 2023 Li, Du, Li, Feng, Chu, Jin and Pan

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.

In humans, variation of the ATP7A gene may cause cranial exostosis, which is similar to “human horn,” but the function of the ATP7A gene in sheep is still unknown. Tissue expression patterns and potential functional loci analysis of the ATP7A gene could help understand its function in sheep horn. In this study, we first identified tissue, sex, breed, and species-specific expression of the ATP7A gene in sheep based on the RNA-sequencing (RNA-seq) data. Second, the potential functional sites of the ATP7A gene were analyzed by using the whole genome sequencing (WGS) data of 99 sheep from 10 breeds. Last, the allele-specific expression of the ATP7A gene was explored. Our result showed the ATP7A gene has significantly higher expression in the big horn than in the small horn, and the ATP7A gene has high expression in the horn and skin, suggesting that this gene may be related to the horn. The PCA results show that the region around the ATP7A can distinguish horned and hornless groups to some extent, further indicating that the ATP7A may be related to horns. When compared with other species, we find seven ruminate specific amino acid sites of the ATP7A protein, which can be important to the ruminate horn. By analyzing WGS, we found 6 SNP sites with significant differences in frequency in horned and hornless populations, and most of these variants are present in the intron. But we still find some potential functional sites, including three missenses, three synonymous mutations, and four Indels. Finally, by combining the RNA-seq and WGS functional loci results, we find three mutations that showed allele-specific expression between big and small horns. This study shows that the ATP7A gene in sheep may be related to horn size, and several potential functional sites we identified here can be useful molecular markers for sheep horn breeding.

 ATP7A sheep horn tissue-specific function section-at-acceptance Livestock Genomics 1 Introduction

Sheep are an important economic resource for humans since they provide a variety of products including meat, fur, skin, and milk (Simoes et al., 2021; Zeder, 2008). Similar to other ruminants, sheep normally have horns which is a bony cranial extension on their craniums (Simon et al., 2022; Wang et al., 2019). Horns serve a dual purpose in the natural world, functioning as a means of self-defense against predators and a primary mechanism for intra-species competition for mates and territory (Buntgen et al., 2018; Nasoori, 2020). However, in the process of breeding, the existence of horns is easy to cause harm to livestock and personnel, which is unfavorable to animal husbandry production management and has great safety risks. At the same time, the removal of horns by artificial means is not only contrary to animal welfare but also increases the workload, and is inherently undesirable. Therefore, farmers invest a lot of effort in selecting hornless animal breeds. Using gene-editing technology to edit specific genes, cows can be hornless, so that calves do not have to go through the painful process of being burned off their horns after birth (Mueller et al., 2021; Wang et al., 2022). But the hornless mutation used in gene-editing breeding in sheep has not yet been confirmed.

The ATP7A gene is an important gene located on the X chromosome, and the ATP7A protein is widely involved in the transport of metal ions and is responsible for the regulation of copper ions (Bitter et al., 2022). The ATP7A protein is mainly distributed in tissues such as the liver and intestines (Gao et al., 2021; Telianidis et al., 2013). The protein is involved in the physiological processes of copper uptake, transport, and storage (Sudhahar et al., 2019), and is an essential component of many important enzymes such as copper-dependent enzymes and cytochrome C oxidases (Guthrie et al., 2020; Horn and Wittung-Stafshede, 2021). In recent years, more and more studies have shown that mutations in the ATP7A gene are closely related to the occurrence of many important diseases. The disease caused by the ATP7A gene defect in the human body is collectively called an ATP7A gene-related disorder, which is mainly manifested as a disorder of copper metabolism (Chen et al., 2020). Mutations in the ATP7A gene cause Menkes disease. Menkes is an infantile, fatal, inherited copper deficiency characterized by progressive neurological damage that eventually causes death. A mild form of Menkes disease is called occipital angle syndrome (OHS) (Moller et al., 2021; Yasmeen et al., 2014), which showed a symptom of exophytic bone warts (De Feyter et al., 2023) that is similar to the bony core of the horn. So we suspect ATP7A can have an important role in the ovine horn. However, currently, there is no research on the ATP7A gene in sheep yet. In this study, RNA-seq data was used to character the expression of the ATP7A gene in different tissues of sheep, for getting the basic function of the ATP7A gene. Then some horn-relative functional loci of the ATP7A gene have been explored by combining WGS data from 10 sheep populations with different horn types. This result will provide some valuable molecular makers in sheep horn breeding.

 2 Materials and methods 2.1 Ethical statement

All experiments involving animals were ethically approved by the Animal Ethics Committee of the Institute of Animal Science, Chinese Academy of Agricultural Sciences (IAS-CAAS) (Beijing, China) on 28 April 2020 (No. IAS 2020-63).

 2.2 Animal and sample collection

The samples used in the RNA-seq study were collected from Tibetan sheep in Dangxiong, Tibet, China. This population consists of four scurred (0–12 cm) and four spiral and horizontally extended horned (>12 cm, SHE) (Supplementary Table S1). The soft-horned tissue of these eight sheep was collected, placed in cryopreservation tubes, and stored in liquid nitrogen.

 2.3 RNA extraction, library construction, and sequencing

The total RNA of tissues was extracted using the Trizol (Thermo Scientific, Wilmington, United States) kit. Electrophoresis using agarose gels at a concentration of 1% to detect RNA degradation and contamination. RNA purity testing was performed using the Kaio K550 (CAIO, Beijing) spectrophotometer. The Qubit® RNA Detection Kit in the Qubit® 2.0 Fluorometer (Life Technologies, California, United States) is used to determine RNA concentration. RNA nano 2100 assay kit from the Agilent Bioanalyzer 6000 System (Agilent Technologies California, United States) for RNA integrity testing.

Take 3 μg of total RNA to construct the library. First, Ribo-ZeroTM GoldKits (Epicentre, United States) was used to remove ribosomal RNA (rRNA). In addition, RNA sequencing libraries were generated according to the manufacturer’s instructions for the Illumina NEB Next Ultra Directional RNA Library Preparation Kit (NEB, Ipswich, United States). Sequencing was then performed on the Hiseq X (Illumina, San Diego, CA, United States) platform.

 2.4 Sequencing data filtering, comparison, assembly and processing

We removed adaptors and low-quality reads using Trim Galore (v.0.6.7). We mapped clean reads to the ARS-UI_Ramb_v2.0 sheep reference genome using paired mapping modules of STAR (v.2.7.7a) with parameters of “--chimSegmentMin 10” and “--outFilterMismatchNmax 3” (Davenport et al., 2022; Dobin et al., 2013). A total of 8 high-quality RNA-seq clean data were obtained for subsequent analysis, with unique mapping reads > 85% and the number of clean reads > 20,000,000 (Supplementary Table S2). We then obtained TPM (Transcripts Per Kilobase of exon model per Million mapped reads) and FPKM (Fragments Per Kilobase of exon model per Million mapped fragments) of genes using StringTie (v.2.1.5) (Kovaka et al., 2019; Pertea et al., 2016; Pertea et al., 2015; Shumate et al., 2022), and extracted raw read counts of them with featureCounts (v.2.0.1) (Liao et al., 2013; Liao et al., 2014). We followed the recommended best practice pipeline in the genome analysis toolkit (GATK) (v.4.2.5.0) and used default settings to perform SNP site prediction (McKenna et al., 2010; Van der Auwera et al., 2013). We conducted allele-specific expression (ASE) analysis using the GATK ASEReadCounter tool. We used Subread_to_DEXSeq (https://github.com/vivekbhr/Subread_to_DEXSeq) of dexseq_prepare_annotation2.py script to format the annotation (GTF) file of the genome, then used load_SubreadOutput.R in Rstudio reads the formatted GTF file and the counts’ matrix output by featureCounts, constructs the DEXSeqDataSetFromFeatureCounts (dds) object, and performs exon difference analysis.

 2.5 Analysis of the expression profile of the <italic>ATP7A</italic> gene

To explore the tissue, sex, breed, and species-specific expression of the ATP7A gene, we used the ggplot2 package of Rstudio to draw the specific expression based on RNA-seq data. The other sheep data was 3889 high-quality RNA seq sequencing obtained after downloading and filtering from NCBI (https://www.ncbi.nlm.nih.gov/) and EBI (https://www.ebi.ac.uk/) (Supplementary Table S3). The data for pigs (2,652 samples) (Supplementary Table S4), cows (4,501 samples) (Supplementary Table S5), and humans (9,810 samples) (Supplementary Table S6) were obtained from their GTEx project (Ardlie et al., 2015; CONSORTIUM, 2020; Liu et al., 2022; Teng et al., 2022). The pig normalized gene expression (TPM) data were obtained in http://piggtex.farmgtex.org/. The cattle TPM data were obtained in https://cgtex.roslin.ed.ac.uk/. The human TPM data were obtained in https://gtexportal.org/home/datasets. After downloading, the mean values were calculated separately for each tissue. We used T-test to calculate the significant difference between scurred and SHE. The Wilcoxon test was used to calculate the significant differences between ram and ewe (Supplementary Table S7). The significant differences between breeds were calculated using variance tests (Supplementary Table S8).

 2.6 Analysis and three-dimensional structure prediction of the ATP7A protein

To investigate the differences in ATP7A protein among different species, we used MEGA (v.5.0) (Kumar et al., 2008) to analyze the homologous evolution tree. Amino acid differences of ATP7A protein were analyzed by the ClustalW algorithm. The amino acid sequences of the ATP7A gene in sheep and other species were obtained from the NCBI gene database (Supplementary Table S9). The 3D structure of the ATP7A gene (TrEMBL: W5QAF8) in sheep is predicated on the UniProt (https://www.uniprot.org). We predicted and labeled the key binding site of sheep ATP7A protein based on the characteristic sequence DKTGT of human P-type ATPase(Sazinsky et al., 2006). Based on the 7si3.1.A P-type Cu(+) transporter model, it is located in amino acids 1044–1048.

 2.7 Whole genome sequencing analysis

To explore the potential functional sites of the ATP7A gene and its 10,000 upstream and downstream base pairs (bp) based on our previous WGS data (Pan et al., 2018), we used plink2 (v2.00a2.3) (Chang et al., 2015) and vcftools (0.1.16) (Danecek et al., 2011) to analyze the allele frequency of each breeds, and analyze the F-statistics (Fst) between horned and hornless breeds. We used the ggplot2 package of Rstudio to draw the principal component analysis (PCA) of allele frequency and the scatter plot of Fst. The data comes from 99 sheep from 10 breeds, including 6 Mongolian breeds: STH (Small Tail Han sheep), WZ (Wuzhumuqin sheep), T (Tan sheep), H (Hu sheep), CB (Cele Black sheep), and BY (Bayinbuluke sheep); 3 Tibetan breeds: VT (Valley Tibetan sheep), PT (Prairie Tibetan sheep), OL (Oula sheep) and a European-originated breed: AM (Australian Merino sheep). Among them, all H and AM have no horns (hornless breeds). Some of STH, T, CB, and BY have horns (horned breeds). VT, PT, and OL have horn (horned breeds).

 3 Results 3.1 The expression of <italic>ATP7A</italic> in sheep with different horn types

The ATP7A (ATPase copper transporting alpha) gene is a protein-coding gene that belongs to the cation transport ATPase (P-type) (TC 3.A.3) family. It is part of the type IB subfamily. Based on gene annotation of ARS-UI_Ramb_v2.0 (GCF_016772045.1) version, ATP7A gene is located on the X chromosome and has two transcripts. The length of the above two transcripts is consistent, both are 153,239 bp. Compared to transcript XM_004022208.5, transcript XM_042242494.1 is missing exon 2. We will use transcript XM_004022208.5 to represent the ATP7A gene in this paper, which has the largest number of exons. From the RNA-seq results (Figure 1A), it can be seen that the FPKM of the ATP7A gene in the SHE group is significantly higher (p < 0.05) than in the scurred group. Exon expression results (Figure 1B) showed that the expression of exons of the ATP7A gene in the SHE group was higher than in the scurred group. We specifically checked the gene expression, GC percent, and repeat of the ATP7A gene (Figure 1C). We observed that the gene had readouts in the exons of each gene annotation, but there were still some expression readouts in regions without gene annotations, indicating that the gene still had some splicing without annotations. For example, there is a long Intron region between Exon 1 and 2. There were some expression readouts located this region in the scurred group and skin, fewer in the SHE group, almost none in muscles, kidneys, and lungs. The same situation exists between Exon 2 and 3. The expression of the ATP7A gene in horn tissues and skin tissues was higher, and in the SHE group is higher than the scurred group. ATP7A gene expression is usually low in muscles, kidneys, and lungs. This suggests that the ATP7A gene has some tissue specificity and may be related to horn function and epithelial function. Interestingly, we observed higher readings within about 5,000 bp (possibly poly-A tail) downstream of the ATP7A gene.

(A) Differential expression of ATP7A gene in scurred and SHE groups, the Expression is the fitted expression estimates from the glm regression, E001–E023 represents the number of Exon (B) differential expression of ATP7A gene exons in scurred and SHE groups (C) differential expression of ATP7A gene in horn tissues and other tissues.

 3.2 The tissue-specific function of the <italic>ATP7A</italic> gene in sheep

The tissue-specific heat map and tau values of the ATP7A gene for different species (Figure 2A) shows that the ATP7A gene has strong tissue specificity. The expression of the ATP7A gene in sheep skin tissue is much higher than in other tissues. At the same time, the gene is highly expressed in pig lung tissue. In humans and cattle, the tissue specificity of this gene is not obvious. Tissue-specific results for ewes and rams (Figure 2B) show sex differences in the expression of the ATP7A gene. The expression of this gene in the pituitary is the highest, and in the ewe is the most significantly higher (p < 0.0001) than the ram. Secondly, the gene is also highly expressed in skin, pituitary, adipose tissue, brain, and kidney. In the skin, the expression of the ewe is significantly higher than that of the ram (p < 0.05). The expression level of ram in the lung was extremely and significantly higher (p < 0.001) than that of the ewe. There were also significant differences between ewes and rams in blood immune, brain, kidney, and liver. The expression of the ATP7A gene in the skin of different breeds of sheep (Figure 2C) showed that the expression level of this gene is the lowest in Assaf sheep and the highest in Tan sheep. The expression of this gene in the hornless breed is significantly lower than that in the horned breed. As a breed mainly used for producing leather, although Hu sheep have no horns, they have thicker skin. This may be the reason for the high expression of the ATP7A gene in Hu sheep.

(A) The difference of RNA-seq expression of ATP7A gene in different species and tissues, the value is the mean of multiple replicates samples. (B) The difference of RNA-seq expression of ATP7A gene in tissues of rams and ewes (C) the difference of RNA-seq expression of ATP7A gene in skin tissues of horned and hornless sheep.

 3.3 The potential function of ATP7A protein

We used MEGA to establish a phylogenetic tree of 19 animal ATP7A proteins using the adjacency method (Figure 3A). The results showed that Artiodactyla is clustered together, and Bovidae and Cervidae are closer (both have horn). At the same time, we found some interesting amino acid sites in the comparison results of amino acid sequences (Figure 3B). The 3, 168, 235, 304, 340, 416, 437, 484, and 1152 amino acids of the ATP7A protein in Bovidae and Cervidae (pink region) have specific amino acid sites, which are different from other animals. The above loci exist in Artiodactyla with horn, rather than in all Artiodactyla (blue region). These amino acid sites may be related to horn. Finally, we predicted the three-dimensional structure and key positions of the ATP7A protein (Figure 3C).

(A) Phylogenetic tree of ATP7A homologous protein sequences (B) the amino acid sites of seven ATP7A proteins unique to ruminants (C) three-dimensional structure and key positions of ATP7A protein.

 3.4 The potential function mutations in the <italic>ATP7A</italic> gene

Based on the results of principal component analysis (PCA) (Figure 4A), hornless Australian merino sheep (AM) and Hu sheep (H) clustered together. Small-tailed sheep (STH) ewes have small horns or no horns, cluster in horned populations, and are relatively close to hornless populations. Among the Tibetan sheep breeds, the Prairie Tibetan sheep (PT) is far from the others, probably because the breed mostly has horns. This result indicates that although the functional sites in the ATP7A gene region cannot completely separate horned or hornless breeds, they can to some extent separate populations with horn differences, indicating that the functional sites in this region may have some regulatory effects on the horn. We screened 16 potential functional sites (Figure 4B, Tables 1, 2 and Supplementary Table S10). Indels4 is located in the upstream region of the promoter, SNP12 and Indels3 are located in the downstream region of the gene, and SNP7, SNP8, SNP9, SNP10, SNP11, Indels1, and Indels2 are located in the intron region. The Fst values of these functional sites are high, indicating that these functional sites differ significantly in horned and hornless sheep populations, which may have some key regulatory effects on horns. The exon region of the ATP7A gene has 6 SNP sites. There are 3 synonymously mutated SNP sites, among which the differentiation of SNP1 and SNP3 is not obvious, while the differentiation of SNP2 is large. There were also 3 SNP sites with missense mutations, among which SNP4 (p.A222E) and SNP5 (p.L33Q) differentiation were not obvious, and SNP6 (p.V10I) differentiation was more obvious, which may play a role in certain groups.

(A) Principal component analysis (PCA) of ATP7A region. Blue: hornless or partially hornless, red: most individuals horned, black: horned Tibetan breeds (B) differential distribution (Fst) of functional sites in the ATP7A region.

Potential SNP functional loci of ATP7A gene.

INFO Mutation type Position Fst Ref Alt Ref allele frequency
AM H BY CB OL PT STH T VT WZ
SNP1 Synonymous 69092941 0.025 C T 1 1 1 1 1 0.9 1 1 1 0.9
SNP2 Synonymous 69092872 0.1875 T C 0.55 0.25 0.7 0.8 0.85 0.5 0.6 0.25 0.6 0.4
SNP3 Synonymous 69075762 0.025 G A 1 1 1 1 1 0.9 1 1 1 0.9
SNP4 Non-synonymous 69053701 0.025 G T 1 1 1 1 1 1 1 0.9 1 0.9
SNP5 Non-synonymous 69031423 0.01875 A T 1 1 1 1 1 1 1 1 0.95 0.9
SNP6 Non-synonymous 69031353 0.125 C T 0.75 0.9 0.9 1 0.95 1 1 1 0.85 0.9
SNP7 Intronic 69080490 0.59375 A C 0.25 0 1 1 1 0.45 1 1 0.2 0.1
SNP8 Intronic 69042600 0.6125 C T 0 0.75 1 1 1 1 1 1 1 0.9
SNP9 Intronic 69027147 0.8625 C T 0 0 0 1 1 1 1 1 1 0.9
SNP10 Intronic 69009769 0.5125 C T 1 1 0 0 1 1 0 0 1 0.9
SNP11 Intronic 68976839 0.6875 C T 1 1 1 0 0 0.5 0 1 0 0
SNP12 Downstream 69112510 0.3125 T C 0.4 0.25 0.8 0.8 0.85 0.55 0.65 0.3 0.65 0.5

Potential Indels functional loci of ATP7A gene.

INFO Mutation type Position Fst Ref allele frequency
AM H BY CB OL PT STH T VT WZ
Indels1 Intronic 69024119 0.51875 1 1 0.05 0 0.15 0.5 0.25 1 1 0.9
Indels2 Intronic 68981839 0.5125 0.55 1 0.25 0.2 0.15 0.45 0.15 0.15 0.4 0.35
Indels3 Downstream 69112032 0.3125 1 0.1 0.9 1 1 1 1 1 1 0
Indels4 Upstream 68950083 0.28125 0.45 0.9 1 1 1 1 1 1 1 0.65
3.5 Allele-specific expression in <italic>ATP7A</italic> gene

We identified three AES (allele-specific expression) loci from RNA-seq data, with ASE1 (chrX: 69092873) located in exon 16 of the ATP7A gene (Figure 5A), and ASE2 (chrX: 69107640) and ASE3 (chrX: 69107842) located downstream of the ATP7A gene (Figure 5B, C). ASE1, ASE2, and ASE3 have both refs (reference) and alt (alternative) alleles in the scurred group. However, ASE1 only has the alt allele in the SHE group, while ASE2 and ASE3 almost only have the ref allele in the SHE group. These three ASE loci may be closely related to horns. Combining WGS and RNA-seq data, we have discovered an interesting phenomenon. ASE1 (chrX: 69092873) and SNP5 (chrX: 69092872) jointly encode amino acid 1040 of ATP7A protein. When these two sites undergo individual mutations, they are both synonymous mutation sites. But when they undergo mutations simultaneously, their encoded Val (valine) will become Ala (alanine). We also found that the ASE3 (chrX: 69107842) was also detected in WGS data. But this site did not show much difference in WGS data.

(A) The ref allele frequency of SNP5 in horned and hornless populations, the allele frequency of ASE1, and the 1040 amino acid coding of the ATP7A gene. (B) Count of alleles in ASE2. (C) Count of alleles in ASE3.

 4 Discussion

The horn is a hollow paired structure located on the frontal bone (Nasoori, 2020). The sheep horn consists of two main parts: the outer layer is the keratinized skin tissue, the horn sheath, and the inner layer is the bony horn core (Huang et al., 2017). Mutations in the ATP7A gene are associated with Menkes disease and occipital angle syndrome (Beyens et al., 2019; Lenartowicz et al., 2015). It is clinically characterized by connective tissue abnormalities and skeletal changes. Occipital angle syndrome, also known as X-linked knife flaccid disease, begins early in adolescence and is characterized by a wedge-shaped calcified nodule at the junction of the trapezius, sternocleidomastoid muscle, and occipital angle (Kaler and Distasio, 1993), X-ray radiographs reveal occipital bone exostoses and hammer-shaped clavicular heads (Kaler, 2014). This abnormal hyperplasia of bone is very similar to the bone core formation process at the horn base. Therefore, we speculate that the ATP7A gene is involved in the formation of the bone core of the horn. From the RNA-seq results, we found that the expression of the ATP7A gene in the large horn group (SHE) is higher than that in the small horn group (scurred) at both the gene level and the exon level, suggesting that the ATP7A gene may also be involved in the regulation of horn size. On the other hand, the ATP7A gene is expressed more in horned cultivars than in hornless breeds, which further proves that ATP7A genes may affect the formation and even size of horns. We found that the ATP7A gene is specifically highly expressed in soft horns and skin tissues, but almost not expressed in other internal organs, combined with previous studies (Panichsillaphakit et al., 2022), indicating that the ATP7A gene is closely related to the function of the skin, may play a role in the keratinization of the skin, and participate in the formation process of horn sheaths. In four animals, humans, cattle, pigs and sheep, the ATP7A gene is only highly expressed in the skin of horned animals, sheep. The ATP7A gene was highly expressed in the skin of both rams and ewes, and the differences were not significant. We also identified some more interesting amino acid variants that were only present in three horned ruminants, sheep, goats and cattle. Taken together, we suggest that the ATP7A gene may be involved in both the formation of the horn bone matrix and the formation of the horn sheath. This suggests that it plays an important role in regulating both the presence and size of the horn.

Based on the analysis of the resequencing data, it was found that to some extent we were able to distinguish between groups with or without horns very clearly. Past studies have found that genes such as RXFP2 and HOXD can control horn phenotype in sheep (Pan et al., 2018; Ren et al., 2016), but these genes only partially explain the phenomenon, suggesting that horned and hornless are not a single point or single gene controlled traits, while the present study found that the ATP7A gene may also control the presence or absence of horns to some extent. We identified seven mutation sites of Fst > 0.5 in the intron region of the ATP7A gene, and three mutation sites of Fst >0.25 within 10,000 bp upstream and 10,000 bp downstream of the gene. These loci can significantly distinguish the horned and hornless populations and may play a huge role in the presence or absence of horns. We identified six SNP loci located in the exonic region of the ATP7A gene, resulting in three synonymous and three missense mutations. SNP6 (p.V10I), is a missense mutation with a high Fst value that may affect the structure and function of ATP7A protein and may play a controlling role in certain populations.

Analyzing the RNA-seq data, we likewise found three more interesting loci. Among them, the SNP located at chrX: 69092873 in the CDS region causes a synonymous mutation, and SNP5 (chrX: 69092872), which co-encode the same amino acid with it, is also a synonymous mutation; both mutations by themselves do not produce changes in the structure of the protein, but when both are mutated at the same time, they change the valine formed by their normal encoding (Val) to alanine (Ala). SNPs located at chrX: 69107640 exhibited allele-specific expression in the scurred and SHE groups, reference alleles could be found in the scurred group, but only alternative alleles existed in the SHE group. The SNP located at chrX: 69107640 also showed the above problem. These two SNP sites regulate the size of the horn. The SNP at chrX: 69107640, which was analyzed in the RNA-seq data, also appeared in the WGS data and was the only SNP site where the RNA-seq results and WGS results intersected.

 5 Conclusion

This study found that the ATP7A gene has the highest expression level in sheep skin tissue and is also highly expressed in soft horn tissue. The ATP7A gene may be related to differences in horn size between ram and ewe and the presence or absence of horns. We have identified multiple potential functional sites, including SNP, Indels, and ASE sites, which may be related to horns. We have discovered seven special amino acid sites of ATP7A protein that may be important to the ruminate horn. This study lays the foundation for future research on the functional mechanism of the ATP7A gene.

 Data availability statement

The datasets presented in this study can be found in online repositories. The names of the repository/repositories and accession number(s) can be found below: https://www.ncbi.nlm.nih.gov/, PRJNA1003277.

 Ethics statement

The animal study was approved by Animal Ethics Committee of the Institute of Animal Science, Chinese Academy of Agricultural Sciences (IAS-CAAS) (Beijing, China). The study was conducted in accordance with the local legislation and institutional requirements.

 Author contributions

HL, YJ, and ZP: study concept, design, and data interpretation. HL, XL, and PF: data acquisition and analysis. HL: experiment and manuscript writing. HL, XD, MC, YJ, and ZP: preparation and manuscript editing, MC, YJ, and ZP: resources. All authors contributed to the article and approved the submitted version.

 Funding

This work was supported by National Key R&D Program of China (2022YFF1000103), National Natural Science Foundation of China (31802031 and 31960659), Agricultural Science and Technology Innovation Program of China (CAAS-ZDRW202106 and ASTIP-IAS13), and China Agriculture Research System of MOF and MARA (CARS-38).

 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.

 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.2023.1239979/full#supplementary-material

References Ardlie K. G. Deluca D. S. Segrè A. V. Sullivan T. J. Young T. R. Gelfand E. T. (2015). Human genomics. The genotype-tissue expression (GTEx) pilot analysis: multitissue gene regulation in humans. Science 348 (6235), 648–660. 10.1126/science.1262110 Beyens A. Van Meensel K. Pottie L. De Rycke R. De Bruyne M. Baeke F. (2019). Defining the clinical, molecular and ultrastructural characteristics in occipital horn syndrome: two new cases and review of the literature. Genes. 10 (7), 528. 10.3390/genes10070528 Bitter R. M. Oh S. Deng Z. Rahman S. Hite R. K. Yuan P. (2022). Structure of the Wilson disease copper transporter ATP7B. Sci. Adv. 8 (9), eabl5508. 10.1126/sciadv.abl5508 Buntgen U. Galvan J. D. Mysterud A. Krusic P. J. Hulsmann L. Jenny H. (2018). Horn growth variation and hunting selection of the Alpine ibex. J. Anim. Ecol. 87 (4), 1069–1079. 10.1111/1365-2656.12839 Chang C. C. Chow C. C. Tellier L. C. Vattikuti S. Purcell S. M. Lee J. J. (2015). Second-generation PLINK: rising to the challenge of larger and richer datasets. GigaScience 4, 7. 10.1186/s13742-015-0047-8 Chen J. Jiang Y. Shi H. Peng Y. Fan X. Li C. (2020). The molecular mechanisms of copper metabolism and its roles in human diseases. Pflugers Arch. 472 (10), 1415–1429. 10.1007/s00424-020-02412-2 Consortium T. G. (2020). The GTEx consortium atlas of genetic regulatory effects across human tissues. Science 369 (6509), 1318–1330. 10.1126/science.aaz1776 Danecek P. Auton A. Abecasis G. Albers C. A. Banks E. DePristo M. A. (2011). The variant call format and VCFtools. Bioinformatics 27 (15), 2156–2158. 10.1093/bioinformatics/btr330 Davenport K. M. Bickhart D. M. Worley K. Murali S. C. Salavati M. Clark E. L. (2022). An improved ovine reference genome assembly to facilitate in-depth functional annotation of the sheep genome. GigaScience 11, giab096. 10.1093/gigascience/giab096 De Feyter S. Beyens A. Callewaert B. (2023). ATP7A-related copper transport disorders: A systematic review and definition of the clinical subtypes. J. Inherit. Metab. Dis. 46 (2), 163–173. 10.1002/jimd.12590 Dobin A. Davis C. A. Schlesinger F. Drenkow J. Zaleski C. Jha S. (2013). Star: ultrafast universal RNA-seq aligner. Bioinformatics 29 (1), 15–21. 10.1093/bioinformatics/bts635 Gao W. Huang Z. Duan J. Nice E. C. Lin J. Huang C. (2021). Elesclomol induces copper-dependent ferroptosis in colorectal cancer cells via degradation of ATP7A. Mol. Oncol. 15 (12), 3527–3544. 10.1002/1878-0261.13079 Guthrie L. M. Soma S. Yuan S. Silva A. Zulkifli M. Snavely T. C. (2020). Elesclomol alleviates Menkes pathology and mortality by escorting Cu to cuproenzymes in mice. Science 368 (6491), 620–625. 10.1126/science.aaz8899 Horn N. Wittung-Stafshede P. (2021). ATP7A-regulated enzyme metalation and trafficking in the Menkes disease puzzle. Biomedicines 9 (4), 391. 10.3390/biomedicines9040391 Huang W. Zaheri A. Jung J. Y. Espinosa H. D. Mckittrick J. (2017). Hierarchical structure and compressive deformation mechanisms of bighorn sheep (ovis canadensis) horn. Acta Biomater. 64, 1–14. 10.1016/j.actbio.2017.09.043 Kaler S. G. DiStasio A. T. (1993). “ATP7A-related copper transport disorders,” in GeneReviews (Seattle, Washington, DC, USA: University of Washington). Kaler S. G. (2014). Translational research investigations on ATP7A: an important human copper ATPase. Ann. N.Y. Acad. Sci. 1314, 64–68. 10.1111/nyas.12422 Kovaka S. Zimin A. V. Pertea G. M. Razaghi R. Salzberg S. L. Pertea M. (2019). Transcriptome assembly from long-read RNA-seq alignments with StringTie2. Genome Biol. 20 (1), 278. 10.1186/s13059-019-1910-1 Kumar S. Nei M. Dudley J. Tamura K. (2008). Mega: A biologist-centric software for evolutionary analysis of DNA and protein sequences. Brief. Bioinform. 9 (4), 299–306. 10.1093/bib/bbn017 Lenartowicz M. Krzeptowski W. Lipinski P. Grzmil P. Starzynski R. Pierzchala O. (2015). Mottled mice and non-mammalian models of Menkes disease. Front. Molec. Neurosci. 8, 72. 10.3389/fnmol.2015.00072 Liao Y. Smyth G. K. Shi W. (2014). featureCounts: an efficient general purpose program for assigning sequence reads to genomic features. Bioinformatics 30 (7), 923–930. 10.1093/bioinformatics/btt656 Liao Y. Smyth G. K. Shi W. (2013). The Subread aligner: fast, accurate and scalable read mapping by seed-and-vote. Nucleic Acids Res. 41 (10), e108. 10.1093/nar/gkt214 Liu S. Gao Y. Canela-Xandri O. Wang S. Yu Y. Cai W. (2022). A multi-tissue atlas of regulatory variants in cattle. Nat. Genet. 54 (9), 1438–1447. 10.1038/s41588-022-01153-5 McKenna A. Hanna M. Banks E. Sivachenko A. Cibulskis K. Kernytsky A. (2010). The genome analysis toolkit: A MapReduce framework for analyzing next-generation DNA sequencing data. Genome Res. 20 (9), 1297–1303. 10.1101/gr.107524.110 Moller L. B. Mogensen M. Weaver D. D. Pedersen P. A. (2021). Occipital horn syndrome as a result of splice site mutations in ATP7A. No activity of ATP7A splice variants missing exon 10 or exon 15. Front. Molec. Neurosci. 14, 532291. 10.3389/fnmol.2021.532291 Mueller M. L. Cole J. B. Connors N. K. Johnston D. J. Randhawa I. Van Eenennaam A. L. (2021). Comparison of gene editing versus conventional breeding to introgress the polled allele into the tropically adapted Australian beef cattle population. Front. Genet. 12, 593154. 10.3389/fgene.2021.593154 Nasoori A. (2020). Formation, structure, and function of extra-skeletal bones in mammals. Biol. Rev. 95 (4), 986–1019. 10.1111/brv.12597 Pan Z. Li S. Liu Q. Wang Z. Zhou Z. Di R. (2018). Whole-genome sequences of 89 Chinese sheep suggest role of RXFP2 in the development of unique horn phenotype as response to semi-feralization. GigaScience 7 (4), giy019. 10.1093/gigascience/giy019 Panichsillaphakit E. Kwanbunbumpen T. Chomtho S. Visuthranukul C. (2022). Copper-histidine therapy in an infant with novel splice-site variant in the ATP7A gene of Menkes disease: the first experience in south east asia and literature review. BMJ Case Rep. 15 (4), e247937. 10.1136/bcr-2021-247937 Pertea M. Kim D. Pertea G. M. Leek J. T. Salzberg S. L. (2016). Transcript-level expression analysis of RNA-seq experiments with HISAT, StringTie and Ballgown. Nat. Protoc. 11 (9), 1650–1667. 10.1038/nprot.2016.095 Pertea M. Pertea G. M. Antonescu C. M. Chang T. C. Mendell J. T. Salzberg S. L. (2015). StringTie enables improved reconstruction of a transcriptome from RNA-seq reads. Nat. Biotechnol. 33 (3), 290–295. 10.1038/nbt.3122 Ren X. Yang G. L. Peng W. F. Zhao Y. X. Zhang M. Chen Z. H. (2016). A genome-wide association study identifies a genomic region for the polycerate phenotype in sheep (Ovis aries). Sci. Rep. 6, 21111. 10.1038/srep21111 Sazinsky M. H. Mandal A. K. Arguello J. M. Rosenzweig A. C. (2006). Structure of the ATP binding domain from the archaeoglobus fulgidus Cu+-ATPase. J. Biol. Chem. 281 (16), 11161–11166. 10.1074/jbc.M510708200 Shumate A. Wong B. Pertea G. Pertea M. (2022). Improved transcriptome assembly using a hybrid of long and short reads with StringTie. PLoS Comput. Biol. 18 (6), e1009730. 10.1371/journal.pcbi.1009730 Simoes J. Abecia J. A. Cannas A. Delgadillo J. A. Lacasta D. Voigt K. (2021). Review: managing sheep and goats for sustainable high yield production. Animal 15 (1), 100293. 10.1016/j.animal.2021.100293 Simon R. Drogemuller C. Luhken G. (2022). The complex and diverse genetic architecture of the absence of horns (polledness) in domestic ruminants, including goats and sheep. Genes. 13 (5), 832. 10.3390/genes13050832 Sudhahar V. Das A. Horimatsu T. Ash D. Leanhart S. Antipova O. (2019). Copper transporter ATP7A (copper-transporting p-type ATPase/menkes ATPase) limits vascular inflammation and aortic aneurysm development: role of microRNA-125b. Arterioscler. Thromb. Vasc. Biol. 39 (11), 2320–2337. 10.1161/ATVBAHA.119.313374 Telianidis J. Hung Y. H. Materia S. Fontaine S. L. (2013). Role of the p-type ATPases, ATP7A and ATP7B in brain copper homeostasis. Front. Aging Neurosci. 5, 44. 10.3389/fnagi.2013.00044 Teng J. Gao Y. Yin H. Bai Z. Liu S. Zeng H. (2022). A compendium of genetic regulatory effects across pig tissues. biorxiv. 10.1101/2022.11.11.516073 Van der Auwera G. A. Carneiro M. O. Hartl C. Poplin R. Del A. G. Levy-Moonshine A. (2013). From fastq data to high confidence variant calls: the genome analysis toolkit best practices pipeline. Curr. Protoc. Bioinforma. 43 (1110), 1–11. 10.1002/0471250953.bi1110s43 Wang H. Zhu H. Hu Z. Heng N. Gong J. Wang Y. (2022). Uncovering novel features of the pc locus in horn development from gene-edited holstein cattle by RNA-sequencing analysis. Int. J. Mol. Sci. 23 (20), 12060. 10.3390/ijms232012060 Wang Y. Zhang C. Wang N. Li Z. Heller R. Liu R. (2019). Genetic basis of ruminant headgear and rapid antler regeneration. Science 364 (6446), eaav6335. eaav6335. 10.1126/science.aav6335 Yasmeen S. Lund K. De Paepe A. De Bie S. Heiberg A. Silva J. (2014). Occipital horn syndrome and classical Menkes syndrome caused by deep intronic mutations, leading to the activation of ATP7A pseudo-exon. Eur. J. Hum. Genet. 22 (4), 517–521. 10.1038/ejhg.2013.191 Zeder M. A. (2008). Domestication and early agriculture in the mediterranean basin: origins, diffusion, and impact. Proc. Natl. Acad. Sci. U. S. A. 105 (33), 11597–11604. 10.1073/pnas.0801317105