Front. Vet. Sci. Frontiers in Veterinary Science Front. Vet. Sci. 2297-1769 Frontiers Media S.A. 10.3389/fvets.2023.1276673 Veterinary Science Original Research Effect of CTSS non-synonymous mutations on litter size in Qianbei Ma goats Zhang Yuan123 Chen Xiang123* Ruan Yong123 Guo Wei123 Chen Jiajing123 Tang Wen123 Ji Quan123 Fu Kaibin123 1Key Laboratory of Animal Genetics, Breeding and Reproduction in the Plateau Mountainous Region, Ministry of Education, Guizhou University, Guiyang, China 2Key Laboratory of Animal Genetics, Breeding and Reproduction of Guizhou Province, Guizhou University, Guiyang, China 3College of Animal Science, Guizhou University, Guiyang, China

Edited by: Sebastián Demyda-Peyrás, University of Córdoba, Spain

Reviewed by: Borhan Shokrollahi, Islamic Azad University, Sanandaj Branch, Iran; Pablo Corva, Universidad Nacional de Mar del Plata, Argentina

 *Correspondence: Xiang Chen, xchen2@gzu.edu.cn 27 11 2023 2023 10 1276673 12 08 2023 23 10 2023 Copyright © 2023 Zhang, Chen, Ruan, Guo, Chen, Tang, Ji and Fu. 2023 Zhang, Chen, Ruan, Guo, Chen, Tang, Ji and Fu

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.

Cathepsin S (CTSS) is a member of the cysteine protease family closely related to reproductive regulation in goats. However, its effect on litter size in goats remains unclear. In this study, the relationship between CTSS gene polymorphisms and litter size was revealed by analyzing the DNA sequence and mRNA expression of CTSS in the gonadal axis of Qianbei Ma goats. In addition, bioinformatics methods were used to evaluate the effect of non-synonymous mutations on CTSS protein structure and function. CTSS was expressed in all parts of the gonadal axis of Qianbei Ma goats, with the highest expression in the uterus in the multi-lamb group and in the fallopian tube in the single-lamb group. The sequencing results showed that four SNPs in CTSS, including g.7413C → T, g.8816A → T, g.9191 T → G and g.10193G → A, were significantly correlated with litter size (p < 0.05). All four analyzed mutation sites were in strong linkage disequilibrium (r2 > 0.33, D′ > 0.70). Additionally, the haplotype Hap1/2 had a significantly higher frequency than the other haplotypes (p < 0.05). g.7413C → T and g.8816A → T were non-synonymous mutations. The g.7413C → T mutation resulted in the substitution of serine 161 of the CTSS protein with phenylalanine (p.S161F), and the g.8816A → T mutation resulted in the substitution of aspartate 219 with tyrosine (p.N219Y). p.S161F was highly conserved across 13 species and that p.N219Y was relatively conserved in cloven-hoofed species. Mutations at two sites changed the local conformation of the CTSS protein, reduced its stability, and affected its function and goat breed evolution. These findings confirm that CTSS affects the lambing traits of goats and provide a theoretical basis for the regulatory mechanism of CTSS in affecting litter size.

 non-synonymous mutation litter size CTSS gene single nucleotide polymorphism goat section-at-acceptance Animal Reproduction - Theriogenology Introduction

The cathepsin (CTS) family of enzymes is widely expressed in various cells and tissues (1), playing important roles in catalyzing protein hydrolysis and regulating various normal biological processes such as cell death, proliferation, migration, cancer development and processing of antigens and antibodies (2, 3). The CTS family contains a variety of subtypes. Cathepsins mainly include cysteine cathepsins, serine cathepsins (cathepsins A, G) and aspartic cathepsins (cathepsins D, Sand E) (4). Activation of pregnancy-specific lysosomal function by CTS in blood leukocytes is highly correlated with interferon-τ (IFNT) expression during maternal–fetal recognition of pregnancy in pregnant cows (5). Cathepsin S (CTSS) is a lysosomal cysteine protease (6). Previous studies on this gene have focused on its role in the immune response, inflammatory response, cardiovascular disease progression and tumor progression (7–9). The activity of this gene is also closely associated with fibronectin degradation and obesity (10). CTSS regulates the secretion of progesterone and estradiol and the proliferation and apoptosis of ovarian granulosa cells in rabbits and is closely related to the regulation of early gestation in goats (11, 12). However, the mechanism underlying the effect of CTSS on litter size in goats is unclear.

Qianbei Ma goats are a unique goat breed raised in the Guizhou Plateau Mountain area of China. This breed is characterized by early sexual maturity, good adaptability, strong disease resistance and stable genetic performance (13) However, its low reproduction rate is a constraint to the development and utilization of this species. Litter size is an important index for quantifying the reproductive performance of female livestock (14). The average lambing rate of Qianbei Ma goats is approximately 207%. However, the breed comprises three groups: high reproductive rate, low reproductive rate and sterile (15). The low heritability of lambing traits in goats limits the traditional methods of selection for high reproductive performance groups. Therefore, it is important to study the expression of CTSS-encoding genes in the Qianbei Ma goat population in order to understand the relationship between CTSS gene polymorphisms and lambing traits and to screen for molecular markers associated with lambing traits to guide the breeding of Qianbei Ma goats.

 Materials and methods Experimental animals

The animals used in this study strictly comply with the guidelines of the Animal Welfare Committee of Guizhou University (EAE-GZU-2022-E030, 25th October, 2022). The Qianbei Ma goats used in the study were obtained from Fuxing Herding Co Ltd., Xishui County, Guizhou Province, China. Hundred and sixty healthy Qianbei Ma ewes with similar body weights were selected and the total number of births and live births of the first, second and third fetuses of the group were recorded. Four milliliters of blood were drawn through the jugular vein into EDTA anticoagulation tubes and stored in a refrigerator at −20°C. Three singleton pregnancy and three multiple pregnancy ewes were selected from 160 Qianbei Ma goats for which reproductive data were recorded, and were euthanized by carotid artery bloodletting after electrocution. Goat gonadal axis tissue samples (including hypothalamus, pituitary, ovaries, uterus and fallopian tubes) were collected within 20 min of euthanasia and washed with phosphate buffered saline solution (PBS). All samples were then rapidly frozen in liquid nitrogen and subsequently transferred to a − 80°C freezer for storage.

 RNA and DNA extraction and cDNA synthesis

Total RNA was extracted from the gonadal axis tissues of the singleton pregnancy and multiple pregnancy groups using TRIzol reagent (Invitrogen, Carlsbad, CA, USA) and a RNeasy RNA purification kit containing DNase treatment (Qiagen, Valencia, CA, USA) according to the manufacturer’s instructions. DNA extraction from the collected blood was performed strictly according to the instructions of the Blood DNA Extraction Kit (Beijing Tiangen Biochemical Technology Co., Ltd., Beijing, China). The concentration and purity of the RNA and DNA were measured using an ultramicro ultraviolet spectrophotometer (NanoDrop2000; Thermo Scientific, Waltham, MA, USA). The samples were tested for integrity on 1% agarose gels, and all samples were stored in a −20°C refrigerator.

 Primer design and synthesis

According to the RNA (accession number: XM_005677657) and DNA (accession number: NC_030810.1) sequences of goat CTSS as published in NCBI, Primer Premier 5.0 (PREMIER Biosoft International, Palo Alto, CA, USA) was used to design primers for amplification. Using β-actin as a fluorescent quantitative internal reference gene, primer sequences were sent to Beijing Tsingke Biotechnology Co., Ltd. for synthesis (Chongqing, China). The primer sequence information is shown in Table 1.

Primer information.

Genes Primer sequence (5′ → 3′) Productsize/bp Annealingtemperature∕°C
CTSS-Exon1,2 F:5’ AGGAAATCACGGAGGAAACCAG 3′R:5’ CCTCAGGATTGAAATATTCAAGCC 3’ 639 63
CTSS-Exon3 F:5’ GTAAAGTCCCTGCTTCCCTCAT 3′R:5’ CCAGGCTCCTATACTATCCATGAA 3’ 555 63
CTSS-Exon4 F:5’ AGAGGAAGAGTTAAGATTGGTGTGC 3′R:5’ GGAAAGTGGTCACAGTGTAGATCAA 3’ 483 63
CTSS-Exon5 F:5’ TCTTCTCCTTCCCGATGTCTGA 3′R:5’ CCTAAGGGACTATGAGATTCACTGC 3’ 457 59
CTSS-Exon6 F:5’ ATTAAAGTTAGACCTTGTTCCGGAG 3′R:5’ CGGCTTGGTGATAAGTTTAGTACAG 3’ 495 61
CTSS-Exon7 F:5’ TCCTCCGTTACTGGTGAAACATAG 3′R:5’ ACACAACTGAACAACAAGCACACA 3’ 642 63
CTSS-Exon8 F:5’ ATAGCATTGAGGGCAAAGAACC 3′R:5’ CTTATTGCTTGATTAGTTCTGGAGG 3’ 478 63
CTSS-Exon9 F:5’ CTCATTCTATGCAGAAGCAGGAGG 3′R:5’ TAATCTGGAGCAGGTGTGAGGAATA 3’ 1,160 63
q-CTSS F:5’ AAGTAGCACGGCGTCTCAT 3′R:5’ TGTCTCCCAGGTGGTTCAT 3’ 114 58
β-actin F:5’ TGATATTGCTGCGCTCGTGGT 3′R:5’ GTCAGGATGCCTCTCTTGCTC 3’ 189 58
PCR amplification and real-time fluorescence quantitative PCR analysis

Total PCR amplification system (20 μL): 10 μL 2× Taq PCR Master Mix (Beijing Tsingke Biological Co., Ltd., Beijing, China), 1 μL DNA template, 1 μL each forward and reverse primers (10 μmol/L), 7 μL deionized water (ddH2O). The PCR procedure was as follows: predenaturation at 98°C for 3 min, denaturation at 98°C for 10 s, annealing at 60°C for 10 s, and extension at 72°C for 15 s. After 35 cycles, the samples were stored at 4°C. After the PCR amplification products were tested by 1% agarose gel electrophoresis to check the expected fragment size, the PCR amplification products were sent to a biological company for sequencing.

The reaction system for fluorescence quantitative PCR (10 μL) contained 5 μL of 2 × UltraSYBR Mixture (Beijing Tsingke Biotechnology Co., Ltd., Beijing, China), 0.5 μL of cDNA, 0.5 μL each of the forward and reverse primers, and ddH2O to 10 μL. The reaction conditions were as follows (see Table 1 for details): 1 cycle at 95°C for 2 min, followed by 40 cycles at 95°C for 15 s, at the appropriate annealing temperature for 30 s, and at 72°C for 30 s. The melting curve was generated automatically by the machine (base temperature 65°C, increasing by 0.5°C every 5 s to 95°C). The annealing temperature of β-actin was the same as the annealing temperature for each experimental gene. The specificity of the PCR primers was confirmed by the presence of a single peak in the melting curve. Three biological replicates were established for each sample.

 Bioinformatics analysis

Evaluation of sequencing results and analysis of polymorphic loci in CTSS were performed using SeqMan and MegAlign in DNAStar (16, 17). Effect of non-synonymous single-nucleotide polymorphisms (nsSNPs) on CTSS protein function were predicted using PhD-SNP and SNAP2 (18, 19). I-Mutant 2.0 and MuPro were used to predict the effect of nsSNPs on protein stability (20, 21). Generation of multiple sequence comparisons based on CTSS amino acid sequences was performed using Cluster Omega for assessing CTSS sequence conservation (22). Sopma was used to analyze the secondary structure of CTSS proteins and AlphaFold2 was used to assess the tertiary structure of wild-type and mutant CTSS proteins (23–25).

 Statistical analysis

The sequencing results were analyzed, the peaks were plotted against one other using SeqMan software, and the identified SNP loci were analyzed statistically. The 2-ΔΔCt method was used to calculate the differential expression levels of the CTSS gene in the ovary, uterus, fallopian tube, pituitary and hypothalamus, and then the expression level of CTSS mRNA in tissues was analyzed by GraphPad Prism 6. Allele frequencies and genotype frequencies were calculated using Haploview 4.2. Population genetic indicators such as polymorphism information content (PIC), gene purity (Ho), effective allele number (Ne), and gene heterozygosity (He) were analyzed according to Chakraborty and Nei (26). Linkage disequilibrium (LD) analysis and haplotype analysis of SNP loci in CTSS were performed using the SHEsis platform.

The experimental data of different genotypes were analyzed using one-way ANOVA in PASW Statistics 18 software to identify associations between different genotypes and reproductive performance, and the analyzed data are expressed as the means ± standard deviations.

 Results Expression profile of <italic>CTSS</italic> in the gonadal axis

As shown in Figure 1, CTSS mRNA expression in the gonadal axis of the multi-lamb ewe population was significantly higher in the uterus than in other gonadal tissues (p < 0.01) and was significantly lower in the pituitary gland than in the ovary, hypothalamus and oviduct (p < 0.05). Analysis between the singleton pregnancy and multiple pregnancy groups showed that CTSS gene expression was significantly higher in the uterus of multiple pregnancy ewes than in the uterus of singleton pregnancy ewes (p < 0.01), and it was significantly higher in the pituitary gland of the single-lambing ewes group than in the pituitary gland of the multilambing ewes (p < 0.05); the expression of CTSS was similar among the remaining tissues, indicating that CTSS plays an important role in the regulation of lambing traits.

Differential expression analysis of CTSS gene in the gonadal axis of single and multi-lamb Qianbei Ma goat. "****" indicated that there were very significant differences among the same tissues of single and multi-lamb ewe (P < 0.0001). "A, B, C" means very significant difference between different tissues of multi-lamb ewes (P < 0.01), "a, b" means significant difference between different tissues of singletons ewes (P < 0.05), and the same letter means no significant differ.

 PCR amplification

In this study, the sizes of the PCR amplification products were consistent with the expected fragment sizes, and the bands appeared clear and bright, with no specific amplification and no obvious trailing phenomenon, confirming the good primer specificity and the ability to be used for direct sequencing (Figure 2).

CTSS gene amplification results. M: DL-2000Marker; P1: exon1 and exon2; P2-P8: exon 2 – exon 8.

 <italic>CTSS</italic> gene polymorphism analysis

The sequencing data were aligned against the CTSS (NC_030810.1) reference sequence using DNAStar software. Two SNP loci, g.7413C → T (exon 6) and g.8816A → T (exon 7), were identified in the exon 6 and 7 regions of the CTSS gene; two SNP loci, g.9191 T → G (intron 7) and g.10193G → A (intron 8), were found in the intron 7 and 8 regions. All of the four SNP loci listed above were present with two alleles and resulted in three genotypes. The sequencing chromatograms are shown in Figure 3.

Sequencing peaks of 4 snp of CTSS gene (n = 160). (A) g.7413C → T, (B) g.8816 A → T, (C) g.9191 T → G, (D) g.10193 G → A.

The identified sequences were aligned to GenBank reference sequences using MegAlign and compared using the Clustal W method. g.7413C → T is the non-synonymous nsSNP leading to the substitution of serine with phenylalanine, and g.8816A → T is the nsSNP leading to the substitution of aspartic acid with tyrosine; g.9191 T → G and g.10193G → A are synonymous mutations (Figure 4).

Amino acid sequence comparison of wild-type and mutant CTSS. Red boxes indicate non-synonymous SNPs.

 Population genetic analysis of <italic>CTSS</italic>

The four SNP loci were genetically characterized, and all four mutant loci had three genotypes. By chi-square test (χ2), the genotype distributions of g.7413C → T, g.8816A → T and g.10193G → A did not deviate from Hardy–Weinberg equilibrium (p > 0.05), while g.9191 T → G deviated from HWE (p < 0.05) (Table 2).

Genotype frequencies and gene frequencies of CTSS in Qianbei Ma goat.

SNPs Genotype frequency Gene frequency χ2 P
g.7413\u00B0C → T CC CT TT C T 0.847 0.357
0.081(13) 0.356(57) 0.563(90) 0.259 0.741
g.8816 A → T AA AT TT A T 3.039 0.081
0.094(15) 0.331(53) 0.575(92) 0.259 0.741
g.9191 T → G TT TG GG T G 4.444 0.035
0.094(15) 0.312(50) 0.594(95) 0.25 0.75
g.10193 G → A GG GA AA G A 3.441 0.064
0.081(13) 0.306(49) 0.613(98) 0.234 0.766

χ20.05 < 5.99; χ20.01 < 9.21, P > 0.05 indicates that the population is in Hardy–Weinberg equilibrium, with sample size in parentheses.

The effective number of alleles per SNP in the CTSS ranged from 1.560 to 1.624, the heterozygosity ranged from 0.359 to 0.384, and the purity ranged from 0.616 to 0.641 (Table 3). The polymorphism level in Qianbei Ma goats was intermediate, ranging from 0.295 to 0.310, indicating that the polymorphic loci were rich in genetic information.

Genetic diversity of CTSS gene SNPs loci in Qianbei Ma goat.

SNPs Effective allele numbers (Ne) Heterozygosity (He) Homozygosity (Ho) Polymorphism information content (PIC)
g.7413C → T 1.624 0.384 0.616 0.310
g.8816 A → T 1.624 0.384 0.616 0.310
g.9191 T → G 1.600 0.375 0.625 0.305
g.10193 G → A 1.560 0.359 0.641 0.295

PIC < 0.25 means low polymorphism, 0.25 < PIC < 0.50 means moderate polymorphism, PIC > 0.5 means high polymorphism.

 LD and haplotype analyses of <italic>CTSS</italic> gene SNPs

Analysis of LD was performed with the four SNPs of the CTSS gene. The results are shown in Figure 5. The SNP loci g.7413C → T, g.8816A → T, g.9191 T → G and g.10193G → A show strong LD (r2 > 0.33, D′ > 0.70).

Linkage disequilibrium analysis. The SPNsg.7413C → T, g.8816A → T, g.9191T → G and g.10193G → A show strong LD (r 2 > 0.33, D’ > 0.70).

Using the SHEsis platform, 2 CTSS haplotypes were identified: Hap1 (−TTGA-) and Hap2 (-CATG-); haplotypes with frequencies <3.% were not involved in the analysis. Hap2 had the highest frequency, accounting for 68.5% of all haplotypes, followed by Hap1 with 17.0% (Table 4). Random combination of two haplotypes produced three diploids, Hap1/1, Hap1/2 and Hap2/2.

CTSS gene haplotypes and frequencies.

Haplotype g.7413C>T g.8816 A>T g.9191 T>G g.10193 G>A Frequency
Hap1 T T G A 0.170
Hap2 C A T G 0.685
Relationship between <italic>CTSS</italic> polymorphisms and litter size

Association analysis of the CTSS gene SNP loci combined with the number of lambs produced in litters 1–3 of Qianbei Ma ewes was performed, and the results are shown in Table 5. The g.7413C → T, g.8816A → T, g.9191 T → G and g.10193G → A SNP loci were significantly correlated with the number of lambs born to Qianbei Ma goats. The frequency of the g.7413 TT genotype was significantly higher than that of CC and CT genotypes at the C → T locus in second births; the frequency of the g.8816 TT genotype was significantly higher than that of the AA genotype at the A → T locus in second births; the frequency of the g.9191 GG genotype was significantly higher than that of the TT genotype at the T → G locus in second births and the frequency of the TG genotype was significantly higher than that of the GG genotype in third births; the frequency of the g.10193 AA genotype was significantly higher than that of the GG genotype at the G → A locus in first births (all p < 0.05).

Association analysis between CTSS gene polymorphism and the number of lambs produced in 1 ~ 3 litters of Qianbei sheep.

SNPs Genotype First-born Second-born Third-born
g.7413C → T CC 1.846 ± 0.375 1.846 ± 0.533b 2.078 ± 0.494
CT 2.035 ± 0.597 2.070 ± 0.529b 2.230 ± 0.732
TT 2.089 ± 0.466 2.289 ± 0.503a 2.200 ± 0.690
g.8816 A → T AA 1.867 ± 0.352 1.933 ± 0.458b 2.133 ± 0.516
AT 1.981 ± 0.604 2.094 ± 0.628ab 2.226 ± 0.669
TT 2.054 ± 0.427 2.315 ± 0.490a 2.174 ± 0.689
g.9191 T → G TT 1.867 ± 0.352 2.000 ± 0.378b 2.133 ± 0.516ab
TG 1.980 ± 0.622 2.040 ± 0.605ab 2.320 ± 0.683a
GG 2.063 ± 0.480 2.290 ± 0.523a 2.116 ± 0.666b
g.10193 G → A GG 1.769 ± 0.439b 2.000 ± 0.577 2.231 ± 0.599
GA 2.000 ± 0.646ab 2.163 ± 0.553 2.225 ± 0.715
AA 2.061 ± 0.450a 2.214 ± 0.561 2.153 ± 0.648

Data in the table are compared in the same column within the same locus, different letters indicate significant differences between genotypes (P < 0.05) and the same letters indicate non-significant differences (P > 0.05).

Association analysis of diploidy with the number of lambs born in litters 1–3 revealed that the frequency of the Hap1/2 genotype was significantly higher than that of the Hap2/2 genotype in the third litter (p < 0.05) and that the Hap1/1 genotype was not significantly correlated with litter size (p > 0.05).

 CTSS bioinformatics analysis

The g.7413C → T mutation in the CTSS gene results in the substitution of serine by phenylalanine (p.S161F) and the g.8816A → T mutation results in the substitution of aspartic acid by tyrosine (p.N219Y). We analyzed the sequence conservation, function, stability, secondary structure and tertiary structure of CTSS proteins for mutant proteins.

Sequence conservation analysis of CTSS proteins

The conserved p.S161F and p.N219Y sites in the CTSS amino acid sequence were analyzed using Cluster Omega online software, and the results are shown in Figure 6. p.S161F is highly conserved in 13 species, and p.N219Y is relatively conserved in even-toed ungulate species. The more highly conserved a site is, the more likely the mutation will have an effect on the structure and function of the protein.

Protein sequence comparison of CTSS in 13 species.

 Effect of non-synonymous mutations on CTSS protein function

The effect of two nsSNPs on protein function was predicted to be neutral using PhD-SNP prediction software, with g.7413C > T scoring 3 on a scale of 0–9 and g.8816A > T scoring 2 on the same scale. According to SNAP2, is the SNPs are neutral, with scores of −38 and − 63, respectively, on a scale of −100 to 100.

 Effect of non-synonymous mutations on the stability of CTSS proteins

Analysis using I-Mutant 2.0 and MuPro showed that the p.S161F mutation increased the stability of the protein, the p.N219Y mutation decreased the stability of the protein, and the two mutations together reduced the stability of the CTSS protein (Tables 6–8).

Association analysis between combined genotypes and the number of lambs produced in 1 ~ 3 litters of Qianbei sheep.

Combined genotype First-born Second-born Third-born
Hap1/1 2.095 ± 0.501 2.257 ± 0.498 2.162 ± 0.683ab
Hap1/2 2.031 ± 0.740 2.125 ± 0.544 2.438 ± 0.669a
Hap2/2 1.714 ± 0.488 2.286 ± 0.488 2.000 ± 0.578b

Data in the table are compared in the same column within the same locus, different letters indicate significant differences between genotypes (P < 0.05) and the same letters indicate non-significant differences (P > 0.05).

Prediction of the influence on CTSS protein function.

Prediction software SNP locus Amino acid mutation locus Prediction results Score
PhD-SNP g.7413C>T p. S161F Neutral 3(0—9)
g.8816 A>T p. N219Y Neutral 2(0—9)
SNAP2 g.7413C>T p. S161F Neutral −38(−100 to 100)
g.8816 A>T p. N219Y Neutral −63(− 100 to 100)

A higher score indicates a greater impact on protein function.

Protein stability prediction.

SNP locus Amino acid mutation locus I-Mutant 2.0 MuPro
Free energy change (DDG)/(kJ▪mol−1) Free energy change (DDG)/(kJ▪mol−1)
g.7413C>T p. S161F 0.21 0.30
g.8816 A>T p. N219Y −0.43 −0.56

Delta delta G (DDG) < 0 means NSSNPS reduce protein stability, and delta delta G (DDG) > 0 indicates that NSSNPS protein stability.

 Effect of non-synonymous mutations on the secondary structure of CTSS proteins

The effects of p.S161F and p.N219Y mutations on the secondary structure of CTSS in goats were analyzed using Sopma, and the results are shown in Table 9. The secondary structures of the wild-type and mutant CTSS proteins contained four structures, namely, theα-helix, extended chain, β-turn and irregular curl, which accounted for 33.23, 17.52, 6.95, and 42.30% of the structures, respectively. The p.S161F and p.N219Y mutations resulted in a decreased proportion of extended chains and β-turns and an increased proportion of irregular curls.

Prediction of the secondary structure of the CTSS mutant.

Type Alpha helix Extended strand Beta turn Random coil
CTSS-wild 33.23% 17.52% 6.95% 42.30%
CTSS-mutant 33.23% 16.92% 6.34% 43.50%
Effect of non-synonymous mutations on the tertiary structure of CTSS proteins

The tertiary structure of a protein is closely related to its function. We used AlphaFold2 to compare the 3D models of wild-type and mutant CTSS proteins at the p.S161F and p.N219Y loci, The model has 100% of amino acid residues in the reasonable region, and the protein structure obtained by the construction has high reliability and can be used as a template for subsequent studies (Figure 7). The constructed model is shown in Figure 8, the p.S161F site contains a non-polar positively charged serine substituted with a non-polar positively charged phenylalanine, and the p.N219Y site contains a polar uncharged aspartic acid substituted with a polar uncharged tyrosine. These two mutations resulted in altered amino acid interactions near the corresponding sites, leading to changes in the structure and function of the mutated protein but little effect on the 3D structure of CTSS.

Predicted Raman profile of protein tertiary structure model reliability.

Simulation of the tertiary structure of CTSS protein, different colors represent different secondary structures, (A) indicates p.S161F site, (B) indicates p.N219Y site, (WT) wild-type, (MU) mutant type.

 Discussion

The lambing trait is one of important reproductive traits, and litter size is low heritability that is influenced by many factors, such as genetics, environment, management, and nutrition (27–29). Although there is a large body of research on the molecular basis of litter size in goats, the practical application of these findings is limited by the complexity of this quantitative trait (30, 31). CTSS polymorphisms are associated with acute atherosclerotic cerebral infarction (32). Mutations in the 5′-untranslated region of the CTSS gene were found to be strongly associated with feed conversion and average daily weight gain in Italian Large White pigs (33). In addition, CTSS is involved in the immune function pathway of high- and low-lambing rate populations in lake sheep, with critical effects on reproduction (34). Therefore, studying the CTSS gene will be beneficial to understand the variation in litter size in Qianbei Ma goats.

Qianbei Ma goats are resistant to adversity and disease, and retain the desirable traits of meekness, low odor, tender meat, and early sexual maturity (35). It is a valuable local breed in Guizhou Province, and the average litter size of it is lower than other domestic goat breeds (36, 37). To determine whether CTSS is a candidate gene for molecular breeding analysis in the high-reproductive performance Qianbei Ma goat population, we examined the expression of CTSS mRNA in the gonadal axis of single and multilamb ewes. In the single-lamb ewe population, CTSS expression was highest in the oviduct, whereas in the multilamb population, CTSS expression was highest in the uterus. CTSS is highly expressed in the sheep uterus; presumably, CTSS may be involved in endometrial remodeling and placenta formation in sheep (38, 39). CTSS induces increases in progesterone and estrogen levels in female rabbits to promote ovarian granulosa cell proliferation; estrogen and ovarian granulosa cells subsequently promote follicle development and ovulation (12, 39). The ovulation rate is an important determinant of litter size, while the uterus is critical for embryo implantation (15, 40, 41). Furthermore, CTSS expression underlies hormonal regulation in maternal tissues, is supportive of embryo implantation and is highly expressed in embryonic trophectoderm apposition sites and non-apposition sites (42, 43). Therefore, the upregulation of CTSS expression in the uterus of multilamb ewes improves lambing numbers in goats by affecting late embryo attachment. Therefore, CTSS may play an important role in reproduction in single and multilamb goats.

To investigate the regulatory mechanism of the CTSS gene in goat reproduction, we evaluated whether CTSS polymorphisms affect lambing traits in Qianbei Ma goats. After extraction of goat DNA, direct sequencing revealed that the genotype distribution of the SNP loci g.7413C → T, g.8816A → T and g.10193G → A did not deviate from Hardy–Weinberg equilibrium (HWE), while that of g.9191 T → G deviated from HWE. Further analysis revealed that all loci were moderately polymorphic (0.295 < PIC <0.310), which may be due to long-term artificial and natural selection (44). In addition, we analyzed the LD of the four mutant loci, and the analysis revealed that all of them were in strong LD (r2 > 0.33, D′ > 0.70). Correlation analysis showed that g.7413C → T, g.8816A → T, g.9191 T → G and g.10193G → A were all associated with litter size. In addition, the third litter number of lambs with Hap1/2 diploid (CTATTGGA) was significantly higher than that with Hap2/2 diploid (CCAATTGG). Therefore, CTSS expression may be closely related to the number of litters produced by goats. In the SNP analysis, we identified the g.7413C → T mutation in the CTSS gene as leading to the substitution of serine by phenylalanine at site 161 (p.S161F) and the g.8816A → T mutation as leading to the substitution of aspartic acid by tyrosine at site 219 (p.N219Y) of the CTSS protein. To elucidate the effect of this non-synonymous mutation on CTSS protein function, bioinformatics analysis revealed that compared with homologs in 12 other species, the p.S161F mutation was highly conserved in all 13 species, and p.N219Y was relatively conserved in even-toed ungulates and less conserved in other species (e.g., Homo sapiens, Canis lupus familiaris, Gallus gallus, and Mustela putorius furo). The mutation had a neutral effect on protein function; but, interestingly, the p.S161F mutation increased the stability of the protein, and the p.N219Y mutation decreased the stability of the protein. It has been shown that missense mutations that increase protein stability may also alter their function, that more stable proteins are more evolved and that mutations at the p.S161F and p.N219Y loci are consistent with speciation (45–47). Studies of the secondary and tertiary structures of the protein showed that mutations resulted in a decreased proportion of extended chains and β-turns and an increased proportion of irregular coiling. It has been found that a reduced β-turn angle and increased irregular coiling can improve a protein’s functional properties (48). Furthermore, our statistical analysis showed that the litter size was significantly greater in the g.7413C → T TT genotype group than in the corresponding CC and CT genotype group; At the g.8816 A → T mutation locus, the TT genotype gave birth to significantly more litter size than the AA genotype. Whether non-synonymous mutations in exons of this gene affect protein function by altering protein stability, thereby further affecting reproductive traits, needs to be determined by more in-depth studies.

 Conclusion

In conclusion, our study showed that uterine CTSS mRNA expression levels in the multilambing ewe population were significantly higher than those in the single-lambing ewe population. Four SNPs loci in the CTSS gene of Qianbei ma goat were significantly associated with litter size, and the g.7413C → T and g.8816A → T mutations were non-synonymous mutations resulting in the substitution of serine 161 with phenylalanine and aspartate 219 with tyrosine in the CTSS protein. Bioinformatic predictions indicated that the p.S161F mutation in CTSS is highly conserved across 13 species, and p.N219Y mutation is relatively conserved in even-toed species; these mutations may significantly reduce the stability of the CTSS protein. These results suggest that the CTSS gene may be closely related to litter size in Qianbei Ma goats. These findings may provide new approaches for the breeding of high-fertility populations of Qianbei Ma goats.

 Data availability statement

The original contributions presented in the study are included in the article/supplementary material, further inquiries can be directed to the corresponding author.

 Ethics statement

The animal studies were approved by Guizhou University Subcommittee of Experimental Animal Ethics. The animal handling procedures were in line with the Chinese Animal Welfare Guidelines and were approved by the Animal Protection and Use Committee of Guizhou University, Guiyang, China (Approval number: EAE-GZU-2022-E030). The studies were conducted in accordance with the local legislation and institutional requirements. Written informed consent was obtained from the owners for the participation of their animals in this study.

 Author contributions

YZ: Conceptualization, Data curation, Methodology, Software, Validation, Writing – original draft, Writing – review & editing. XC: Conceptualization, Funding acquisition, Methodology, Project administration, Writing – review & editing. YR: Methodology, Software, Writing – review & editing. JC: Investigation, Writing – review & editing. WT: Investigation, Validation, Writing – review & editing. QJ: Validation, Writing – review & editing. KF: Writing – review & editing. WG: Review, Editing & supervision.

 Funding

The author(s) declare financial support was received for the research, authorship, and/or publication of this article. This research was funded by the National Key Research and Development Plan of China (No. 2021YFD1200403), the Science and Technology Project of Guizhou Province (Qian Kehe Foundation—ZK [(2021)] General 151), and Guizhou High-Level Innovative Talents Project (Qian Kehe Platform Talents [(2022)] 021–1).

 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.

 References Man SM Kanneganti TD. Regulation of lysosomal dynamics and autophagy by ctsb/cathepsin b. Autophagy. (2016) 12:2504–5. doi: 10.1080/15548627.2016.1239679 Patel S Homaei A El-Seedi HR Akhtar N. Cathepsins: proteases that are vital for survival but can also be fatal. Biomed Pharmacother. (2018) 105:526–32. doi: 10.1016/j.biopha.2018.05.148, PMID: 29885636 Cheng XW Huang Z Kuzuya M Okumura K Murohara T. Cysteine protease cathepsins in atherosclerosis-based vascular disease and its complications. Hypertension. (2011) 58:978–86. doi: 10.1161/HYPERTENSIONAHA.111.180935, PMID: 21986502 Rossi A Deveraux Q Turk B Sali A. Comprehensive search for cysteine cathepsins in the human genome. Biol Chem. (2004) 385:363–72. doi: 10.1515/BC.2004.040, PMID: 15195995 Talukder M Balboula AZ Shirozu T Kim SW Kunii H Suzuki T . Activation of lysosomal cathepsins in pregnant bovine leukocytes. Reproduction. (2018) 155:515–28. doi: 10.1530/REP-18-0078, PMID: 29626104 Fei M Zhang L Wang H Zhu Y Niu W Tang T . Inhibition of cathepsin s induces mitochondrial apoptosis in glioblastoma cell lines through mitochondrial stress and autophagosome accumulation. Front Oncol. (2020) 10:516746. doi: 10.3389/fonc.2020.516746, PMID: 33425712 Kim N Bae KB Kim MO Yu DH Kim HJ Yuh HS . Overexpression of cathepsin s induces chronic atopic dermatitis in mice. J Invest Dermatol. (2012) 132:1169–76. doi: 10.1038/jid.2011.404, PMID: 22170489 Klinngam W Janga SR Lee C Ju Y Yarber F Shah M . Inhibition of cathepsin s reduces lacrimal gland inflammation and increases tear flow in a mouse model of Sjogren's syndrome. Sci Rep. (2019) 9:9559. doi: 10.1038/s41598-019-45966-7, PMID: 31267034 Lai CH Chang JY Wang KC Lee FT Wu HL Cheng TL. Pharmacological inhibition of cathepsin s suppresses abdominal aortic aneurysm in mice. J Vasc Surg. (2020) 59:990–9. doi: 10.1080/15548627.2016.1239679 Taleb S Cancello R Clement K Lacasa D. Cathepsin s promotes human preadipocyte differentiation: possible involvement of fibronectin degradation. Endocrinology. (2006) 147:4950–9. doi: 10.1210/en.2006-0386, PMID: 16825321 Fu K Chen X Guo W Zhou Z Zhang Y Ji T . Effects of n acetylcysteine on the expression of genes associated with reproductive performance in the goat uterus during early gestation. Animals (Basel). (2022) 58:978–86. doi: 10.3390/ani12182431 Song G Jiang Y Wang Y Song M Niu X Xu H . Modulation of cathepsin s (ctss) regulates the secretion of progesterone and estradiol, proliferation, and apoptosis of ovarian granulosa cells in rabbits. Animals (Basel). (2021) 385:363–72. doi: 10.3390/ani11061770 Tang W Xu QH Chen X Guo W Ao Z Fu K . Transcriptome sequencing reveals the effects of circrna on testicular development and spermatogenesis in qianbei ma goats. Frontiers in veterinary science. (2023) 10:1167758. doi: 10.3389/fvets.2023.1167758, PMID: 37180060 Sasaki Y Hayashi Y Murano S Kohigashi T. Quantitative relationship between the number of cross-fostering piglets and subsequent productivity of sows on commercial swine farms. Anim Sci J. (2022) 93:e13752. doi: 10.1111/asj.13752, PMID: 35726538 Zhou Z Chen X Zhu M Wang W Ao Z Zhao J . Cathepsin d knockdown regulates biological behaviors of granulosa cells and affects litter size traits in goats. J Zhejiang Univ Sci. (2021) 22:893–905. doi: 10.1631/jzus.B2100366, PMID: 34783220 Tsiodras S Gold HS Sakoulas G Eliopoulos GM Wennersten C Venkataraman L . Linezolid resistance in a clinical isolate of staphylococcus aureus. Lancet. (2001) 358:207–8. doi: 10.1016/S0140-6736(01)05410-1 Pal A Chakravarty AK Chatterjee PN. Polymorphism of growth hormone gene and its association with seminal and sexual behavioral traits in crossbred cattle. Theriogenology. (2014) 81:474–80. doi: 10.1016/j.theriogenology.2013.11.002, PMID: 24315684 Capriotti E Calabrese R Casadio R. Predicting the insurgence of human genetic diseases associated to single point protein mutations with support vector machines and evolutionary information. Bioinformatics. (2006) 22:2729–34. doi: 10.1093/bioinformatics/btl423, PMID: 16895930 Hecht M Bromberg Y Rost B. Better prediction of functional effects for sequence variants. BMC Genomics. (2015) 16:S1. doi: 10.1186/1471-2164-16-S8-S1, PMID: 26110438 Capriotti E Fariselli P Casadio R. I-mutant2.0: predicting stability changes upon mutation from the protein sequence or structure. Nucleic Acids Res. (2005) 33:W306–10. doi: 10.1093/nar/gki375, PMID: 15980478 Cheng J Randall A Baldi P. Prediction of protein stability changes for single-site mutations using support vector machines. Proteins. (2006) 62:1125–32. doi: 10.1002/prot.20810, PMID: 16372356 Sievers F Wilm A Dineen D Gibson TJ Karplus K Li W . Fast, scalable generation of high-quality protein multiple sequence alignments using clustal omega. Mol Syst Biol. (2011) 7:539. doi: 10.1038/msb.2011.75, PMID: 21988835 Geourjon C Deleage G. Sopma: significant improvements in protein secondary structure prediction by consensus prediction from multiple alignments. Comput Appl Biosci. (1995) 11:681–4. doi: 10.1093/bioinformatics/11.6.681, PMID: 8808585 Jumper J Evans R Pritzel A Green T Figurnov M Ronneberger O . Highly accurate protein structure prediction with alphafold. Nature. (2021) 596:583–9. doi: 10.1038/s41586-021-03819-2, PMID: 34265844 Varadi M Anyango S Deshpande M Nair S Natassia C Yordanova G . Alphafold protein structure database: massively expanding the structural coverage of protein-sequence space with high-accuracy models. Nucleic Acids Res. (2022) 50:D439–44. doi: 10.1093/nar/gkab1061, PMID: 34791371 Chakraborty R Nei M. Bottleneck effects on average heterozygosity and genetic distance with the stepwise mutation model. Evolution. (1977) 31:347–56. doi: 10.2307/2407757 Bajhau HS Kennedy JP. Influence of pre- and postpartum nutrition on growth of goat kids. Small Rumin Res. (1990) 3:227–36. doi: 10.1016/0921-4488(90)90040-D Mellado M Veliz FG Macias-Cruz U Avendano-Reyes L Garcia JE Rosales-Nieto CA. Effect of breed and management practices on reproductive and milking performance of rangeland goats. Trop Anim Health Prod. (2022) 54:193. doi: 10.1007/s11250-022-03193-9, PMID: 35650454 Naicy T Venkatachalapathy RT Aravindakshan TV Raghavan KC Mini M Shyama K. Relative abundance of tissue mRNA and association of the single nucleotide polymorphism of the goat ngf gene with prolificacy. Anim Reprod Sci. (2016) 173:42–8. doi: 10.1016/j.anireprosci.2016.08.009, PMID: 27576174 Ward MA Abuargob MO Hdud MI Elgusbi MT Ruban YS. Molecular genetic aspects of the pou 1 f 1 gene in Holstein reproductive traits. J Trend Sci Res Develop. (2018) 2:2567–70. doi: 10.1002/prot.20810 Raudsepp T Chowdhary BP. Chromosome aberrations and fertility disorders in domestic animals. Ann Rev Anim Biosci. (2016) 4:15–43. doi: 10.1146/annurev-animal-021815-111239, PMID: 26884101 Luo L Zhu M Zhou J. Association between ctss gene polymorphism and the risk of acute atherosclerotic cerebral infarction in Chinese population: a case-control study. Biosci Rep. (2018) 11:681–4. doi: 10.1042/BSR20180586 Fontanesi L Speroni C Buttazzoni L Scotti E Costa LN Davoli R . Association between cathepsin l (ctsl) and cathepsin s (ctss) polymorphisms and meat production and carcass traits in Italian large white pigs. Meat Sci. (2010) 85:331–8. doi: 10.1016/j.meatsci.2010.01.023, PMID: 20374908 Feng X Li F Wang F Zhang G Pang J Ren C . Genome-wide differential expression profiling of mrnas and lncrnas associated with prolificacy in hu sheep. Biosci Rep. (2018) 7:D439–D444. doi: 10.1042/BSR20171350, PMID: 29439142 Sui MX Wang HH Wang ZW. Molecular cloning, polymorphisms, and expression analysis of the rerg gene in indigenous chinese goats. Genet Mol Res. (2015) 14:14936–46. Cai HF Chen Z Luo WX. Associations between polymorphisms of the gfi1b gene and growth traits of indigenous Chinese goats. Genet Mol Res. (2014) 13:872–80. doi: 10.4238/2014.February.13.5, PMID: 24615051 Yao X Yang F El-Samahy MA Liu B Zhao B Gao X . Identification and characterization of unique and common lncRNAs and mRNAs in the pituitary, ovary, and uterus of hu sheep with different prolificacy. Genomics. (2022) 114:110511. doi: 10.1016/j.ygeno.2022.110511, PMID: 36283658 Song G Spencer TE Bazer FW. Cathepsins in the ovine uterus: regulation by pregnancy, progesterone, and interferon tau. Endocrinology. (2005) 146:4825–33. doi: 10.1210/en.2005-0768, PMID: 16099855 Hong L Chen X Zhu M Ao Z Tang W Zhou Z. Timp1 may affect goat prolificacy by regulating biological function of granulosa cells. Arch Anim Breed. (2022) 65:105–11. doi: 10.5194/aab-65-105-2022, PMID: 35320991 Pramod RK Sharma SK Singhi A Pan S Mitra A. Differential ovarian morphometry and follicular expression of bmp15, gdf9 and bmpr1b influence the prolificacy in goat. Reprod Domest Anim. (2013) 48:803–9. doi: 10.1111/rda.12165, PMID: 23581245 Reese J Das SK Paria BC Lim H Song H Matsumoto H . Global gene expression analysis to identify molecular markers of uterine receptivity and embryo implantation. J Biol Chem. (2001) 276:44137–45. doi: 10.1074/jbc.M107563200, PMID: 11551965 Baston-Buest DM Schanz A Buest S Fischer JC Kruessel JS Hess AP. The embryo's cystatin c and f expression functions as a protective mechanism against the maternal proteinase cathepsin s in mice. Reproduction (Cambridge, England). (2010) 139:741–8. doi: 10.1530/REP-09-0381, PMID: 20093401 Matsukawa S Sumigama S Kotani T Wang J Miki R Moriyama Y . Possible association between cathepsin v and the development of placenta accreta spectrum disorders. Gynecol Obstet Investig. (2019) 84:396–406. doi: 10.1159/000496609, PMID: 30759440 Islam R Liu Z Li Y Jiang L Ma Y. Conservation assessment of the state goat farms by using SNP genotyping data. Genes (Basel). (2020) 11:652. doi: 10.3390/genes11060652 Li X Yang J Shen M Xie XL Liu GJ Xu YX . Whole-genome resequencing of wild and domestic sheep identifies genes associated with morphological and agronomic traits. Nat Commun. (2020) 11:2815. doi: 10.1038/s41467-020-16485-1, PMID: 32499537 Bloom JD Labthavikul ST Otey CR Arnold FH. Protein stability promotes evolvability. Proc Natl Acad Sci U S A. (2006) 103:5869–74. doi: 10.1073/pnas.0510098103, PMID: 16581913 Sobitan A Edwards W Jalal MS Kolawole A Ullah H Duttaroy A . Prediction of the effects of missense mutations on human myeloperoxidase protein stability using in silico saturation mutagenesis. Genes (Basel). (2022) 13:1412. doi: 10.3390/genes13081412, PMID: 36011324 Zhang W Zhao P Li J Wang X Hou J Jiang Z. Effects of ultrasound synergized with microwave on structure and functional properties of transglutaminase-crosslinked whey protein isolate. Ultrason Sonochem. (2022) 83:105935. doi: 10.1016/j.ultsonch.2022.105935, PMID: 35114555