Skip to main content

Genome-wide run of homozygosity analysis reveals candidate genomic regions associated with environmental adaptations of Tibetan native chickens



In Tibet, the two most important breeds are Tibetan chicken and Lhasa white chicken, and the duo exhibit specific adaptations to the high altitude thereby supplying proteins for humans living in the plateau. These breeds are partly included in the conservation plans because they represent important chicken genetic resources. However, the genetic diversity of these chickens is rarely investigated. Based on whole-genome sequencing data of 113 chickens from 4 populations of Tibetan chicken including Shigatse (SH), Nyemo (NM), Dagze (DZ) and Nyingchi (LZ), as well as Lhasa white (LW) chicken breed, we investigated the genetic diversity of these chicken breeds by genetic differentiation, run of homozygosity (ROH), genomic inbreeding and selection signature analyses.


Our results revealed high genetic diversity across the five chicken populations. The linkage disequilibrium decay was highest in LZ, while subtle genetic differentiation was found between LZ and other populations (Fst ranging from 0.05 to 0.10). Furthermore, the highest ROH-based inbreeding estimate (FROH) of 0.11 was observed in LZ. In other populations, the FROH ranged from 0.04 to 0.06. In total, 74, 111, 62, 42 and 54 ROH islands containing SNPs ranked top 1% for concurrency were identified in SH, NM, DZ, LZ and LW, respectively. Genes common to the ROH islands in the five populations included BDNF, CCDC34, LGR4, LIN7C, GLS, LOC101747789, MYO1B, STAT1 and STAT4. This suggested their essential roles in adaptation of the chickens. We also identified a common candidate genomic region harboring AMY2A, NTNG1 and VAV3 genes in all populations. These genes had been implicated in digestion, neurite growth and high-altitude adaptation.


High genetic diversity is observed in Tibetan native chickens. Inbreeding is more intense in the Nyingchi population which is also genetically distant from other chicken populations. Candidate genes in ROH islands are likely to be the drivers of adaptation to high altitude exhibited by the five Tibetan native chicken populations. Our findings contribute to the understanding of genetic diversity offer valuable insights for the genetic mechanism of adaptation, and provide veritable tools that can help in the design and implementation of breeding and conservation strategies for Tibetan native chickens.

Peer Review reports


The Tibetan plateau is the largest high-altitude area on earth with an average altitude exceeding 4000 m, representing 25% of the landmass of China. This high-altitude environment becomes a habitat for many unique animal genetic resources. Tibetan chicken is one of the chicken breeds found in this area dating back to the seventh century A.C [1]. The chickens are widely distributed in farming areas of Tibet, including Shigatse, Lhasa, Lhoka and Nyingchi. Beside Tibetan chickens, Lhasa white chicken is another widespread native breed reared on the plateau, which is bred for several decades by intercrossing male White Leghorn and female Tibetan chickens. The basis of physiological and genetic adaptations to the extreme environmental conditions of the Tibetan plateau have recently been partly investigated in Tibetan chickens [1,2,3], and the genomic analysis indicated that the chickens might be a composite of multiple distinct populations [2, 4]. However, population genomic analysis was rarely conducted to explore the diversity of the chicken populations reared in Tibet and to guide genetic resource conservation and utilization efforts.

The availability of high-throughput affordable sequencing techniques has allowed genome-wide analysis of the genetic structure and relationships in animal populations. Large scale omics data have opened new perspectives for a more accurate animal genetic analysis. In addition to SNP and gene expression data, runs of homozygosity (ROH) has joined the assembly of omics data available in biological databases for livestock gene discovery and diversity assessment. ROH is a kind of long continuous homozygous stretches in the genome that are formed by the combination of two identical haplotypes in an individual [5]. These regions are ubiquitous even in outbred populations, and are usually considered to be the index of autozygosity. Long homozygous regions throughout the genome result from demographic events, mating between close relatives (population bottleneck), reduction in population size (leading to be more likely exposed to genetic drift), selection (breeding) and small inversions that suppress recombination events. Thus, population demography, structure and diversity can be explored based on the distribution and location of ROH regions of the genome. Moreover, previous studies have shown that ROH-based inbreeding estimates provide a better measure of individual autozygosity than pedigree-based estimates of overall inbreeding if kinship between founder animals are not accounted [6]. The occurrence of ROH in individual genomes has also facilitated our understanding of the genetic basis of complex phenotypes. ROH analyses enabled the study of genomic regions with high incidence of homozygosity across individuals, which were first referred to as ROH islands [7] within a population. In addition, the ROH islands can provide important insights into the population signatures of positive selection due to the linkage disequilibrium (LD) [8]. Hence, overlapping ROH islands across populations and species are valuable in comparative genomic studies and may reveal critical genetic regions for complex traits. Consequently, ROH analyses are becoming complementary to genome-wide association studies in the detecting population-specific major genes in humans and animals [9].

In chickens, analysis of ROH has been deployed to assess the genome-wide diversity in local and imported genetic resources. Results from such analysis help in the design and review of effective conservation program for endangered populations [10,11,12,13]. The long-term selection molded the presence of ROHs and their associated genomic regions resulting in unique population adaptation to environment-imposed challenges in broilers, suggesting that ROHs might also be a product of selection events rather than demographic and population history [14]. Moreover, multiple candidate genes involved in growth, egg production, disease resistance and behavior were associated with the ROH islands in different chicken breeds [15, 16].

Within-breed diversity and the population structure analyses in livestock species are fundamental for understanding environmental adaptation, implementing conservation programs and designing selection plans [10, 17, 18]. While a huge effort was expended to study cosmopolitan breeds in the past, a growing attention had been shifted to the local breeds which are important genetic resources for their potential to solving problems in agriculture related to environmental changes [19, 20]. Local chickens in Tibet have evolved over centuries under extreme natural conditions. They may serve as a great reservoir of the genetic pool for the identification of genes under natural and artificial selection particularly those harboring putative signatures of environmental adaptation. Herein, four Tibetan chicken populations and a local cultivated breed (Lhasa white) were sampled in this study during 2016 and 2017. All the Tibetan chicken populations were raised traditionally by local farmers, while Lhasa white was managed at the Institute of Animal Husbandry and Veterinary, Tibetan Academy of Agricultural and Animal Husbandry Sciences. The objectives of the present study were to (i) evaluate the genetic diversity of Tibetan chickens reared in different areas of the Tibetan plateau using whole-genome sequencing data, (ii) detect ROH within each chicken population and evaluate the genomic inbreeding and (iii) reveal the genomic regions of ROH islands that may influence the adaptation of Tibetan native chicken to high altitude. The results are expected to provide valuable information for the gene-phenotype association, as well as for the conservation of chicken genetic resources in Tibet.


Summary of the genetic diversity parameters

The genetic diversity for Shigatse (SH), Neymo (NM), Dagze (DZ), Nyingchi (LZ) and Lhasa white (LW) chicken populations was evaluated by observed heterozygosity (Ho), expected heterozygosity (He) and multiple-locus heterozygosity (MLH) using eligible SNPs under Hardy–Weinberg equilibrium. Ignoring the minor allele frequency (MAF) of SNPs, the Ho ranged from 0.15 to 0.18. Using filtered SNPs with MAF ≥ 5%, Ho increased ranging from 0.27 to 0.31. The estimates of He were similar to Ho for each population. MLH ranged from 0.22 to 0.31 when all SNPs were used, and it went up from 0.27 to 0.31 when only SNPs with MAF ≥ 5% were used. The MAF ranged from 0.11 to 0.13 when its filter was set at 0, and it went up from 0.21 to 0.22 when only SNPs with MAF ≥ 5% were used. The diversity indices of the five populations were shown in Table 1. The linkage disequilibrium (LD) decay pattern was different in the five populations (Fig. S1a). Strong LD was observed between SNPs within a short range in all five populations. For SNPs up to 5 kb apart, the average r2 values were 0.08, 0.07, 0.08, 0.20 and 0.09 for SH, NM, DZ, LZ and LW populations, respectively. This indicated a stronger LD between SNPs and a clear difference in LZ compared to the other four chicken populations. Furthermore, the analysis of population relatedness (Table S1) revealed that Fst values were small for all pairwise comparisons. The lowest value of 0.001 was observed in SH vs. LW while the highest value of 0.095 was observed in DZ vs. LZ. Concordantly, the LZ population was moderately distant from other populations with a larger Fst value ranging from 0.052 to 0.095, while the LW population was a synthesized line which seemed to be genetically closest to all Tibetan chicken populations. All eligible SNPs filtered for HWE and call rate were used for principal component analysis (PCA), which revealed a cluster separation between LW and the other 4 Tibetan chicken population. SH, LZ and DZ were also clearly separated by PC1 and PC2, except for a few scattered individuls. Chickens from NM population mixed together with DZ, occupying an intermediate position between SH and DZ populations (Fig. S1b).

Table 1 Observed (Ho), expected (He) heterozygosity and minor allele frequency (MAF) for each chicken population

Runs of homozygosity within the population

The current study identified 1269, 2438, 1366, 1284 and 1342 ROHs in SH, NM, DZ, LZ and LW chicken population, respectively. ROHs were identified in each sampled bird. Population-wise, the average number of ROH segments was highest in the LZ (128.4 ROH/bird) compared to the other populations. The lowest and highest average length of ROHs were observed in DZ (54.64 Mb) and LZ (102.54 Mb), respectively. The number of SNPs harbored in the ROHs varied between the studied populations, and the maximum number (22,386) was located on GGA1 of the SH population. Similarly, the minimum number of SNPs in the ROHs (50) was located on GGA1 found in the LZ chicken population (Table 2). As shown in Fig. 1a and b, ROHs identified in 106 birds were mainly distributed across the GGA1 to GGA15, GGA17 to GGA28 and GGA33 of the chicken genome. However, the majority of ROH regions were clustered in the macrochromosomes (GGA1 ~ GGA9). The classification of ROHs based on size shown that the short ROHs (< 1 Mb) were predominant, accounting for 79.75 to 86.09% of ROHs across all populations (Fig. 1c). For a better view, we plotted the ROH number against their size for each bird (Fig. 1d). We observed a consistently high correlation in NM (r = 0.87), DZ (r = 0.92), LZ (r = 0.91) and LW (r = 0.84), while correlation considerably varied in SH (r = 0.71) chickens. Moreover, the bird with an extremely long ROH (277.112 Mb) was found in the SH population, while the bird with the shortest ROH (0.831 Mb) belonged to the DZ population.

Table 2 Descriptive statistics of run of homozygosity for five chicken populations
Fig. 1
figure 1

Descriptive graphics of run of homozygosity (ROH) in 5 Tibetan native chicken populations. a The average number of ROHs per chromosome (bars) and the average percentage of each chromosome covered by ROHs (lines) for all chickens. b The distribution and density of ROHs in the whole genome for all chickens. The different color represents different number of ROHs within 18 Mb window size. c The percentage of total ROH within each ROH length category, including short (< 1 Mb), medium (1–3 Mb), and long (> 3 Mb) per chicken population. d Total number of ROHs and total length of genome (Mb) covered by ROH segments per birds for each chicken population

Genomic inbreeding coefficients

Genomic inbreeding was evaluated by the proportion of the genome within ROH (FROH), genomic SNP-by-SNP inbreeding coefficient (FGRM), excess of homozygosity (FHOM) and correlation between uniting gametes (FUNI). As shown in Fig. 2a, the four genomic inbreeding coefficients were small and varied across the five populations. These coefficients were consistent within each population except for LZ, which showed high FROH (0.110 ± 0.016) and low FHOM (− 0.029 ± 0.055), low FGRM (− 0.027 ± 0.046) and low FUNI (− 0.027 ± 0.052). We further analyzed the correlation among the four inbreeding coefficients which was 0.55 between FROH and FHOM, 0.51 between FROH and FUNI, 0.91 between FHOM and FUNI, 0.39 between FGRM and FUNI, respectively (Fig. 2b). However, the correlations between FROH and FGRM, and between FHOM and FGRM were not significant. Moreover, birds in the LZ population had consistent inbreeding coefficients for each index, whereas large inter-bird variability existed for each index in the LW population.

Fig. 2
figure 2

Genomic inbreeding and their correlation in 5 Tibetan native chicken populations. a Genomic inbreeding coefficients inferred from the proportion of the genome within ROH (FROH), genomic SNP-by-SNP inbreeding coefficient (FGRM), excess of homozygosity (FHOM) and correlation between uniting gametes (FUNI). b The correlation between each of 2 genomic inbreeding coefficients across birds. The scatter plot was distinguished by chicken population

Gene annotation of ROHs

We detected 74, 111, 62, 42 and 54 ROH islands that ranked top 1% in SH, NM, DZ, LZ and LW chicken population, respectively (Fig. 3). Annotated genes within the ROH islands were retrieved from the Ensembl genome browser, resulting in 316, 491, 259, 197 and 166 genes in SH, NM, DZ, LZ and LW, respectively. Gene ontology (GO) analysis revealed that these genes were significantly enriched in the biological processes including positive regulation of synapse assembly, positive regulation of I-kappaB kinase/NF-kappaB signaling, osteoblast differentiation, cellular response to amino acid stimulus, cell adhesion and endodermal cell differentiation (Table S2). Among these genes, 11 were common to all the five populations and were located on GGA5 and GGA7. These genes include BDNF, CCDC34, LGR4, LIN7C, GLS, LOC101747789, MYO1B, STAT1 and STAT4 (Fig. 3). In addition, the top 1% ROH islands were mapped to the chicken QTL database and the ROH islands overlapped with 26, 29, 21, 19 and 17 QTLs in SH, NM, DZ, LZ and LW population, respectively (Table S3). Common QTLs that overlapped with the top 0.1% ROH islands in the five populations were for comb and ileum weight. Moreover, the ROH islands that harbored QTL for the ovary weight and percentage and wattle weight and length were specifically detected in the Tibetan chicken (Table 3).

Fig. 3
figure 3

Circular Manhattan plot incidence of each SNP in run of homozygosity (ROH) for 5 Tibetan native chicken populations. From inside to outside, circles denote Shigatse, Nyemo, Dagze, Ningychi and Lhasa white population, respectively. The outermost circle denotes the chromosome. The shared genes harbored in the top 1% ROH islands by five populations were shown in red. The y axis denotes the frequency (%) of SNP that occurred in ROH

Table 3 Top 0.1% ROH islands overlapped with reported QTLs for 5 populations

Selection signature analysis

Notably, the genomic region containing common QTLs spanned from 0.03 Mb to 1.13 Mb of GGA8, and harbored six top 0.1% ROH islands across the five populations. By focusing on GGA8, we calculated integrated haplotype homozygosity (iHS) for each population. In all the populations except for the DZ population, SNPs harbored in the top 0.1% ROH islands were strongly selected for (Fig. 4), and the average |iHS| values of SNPs in each ROH island were higher than the average value of the total SNPs on GGA8. There were 1, 1, 8 and 2 SNPs (P-value ranked top 0.1%) that harbored signatures of selection in SH, NM, LZ and LW population, respectively. Further mining of this region revealed that only three genes (AMY2A, NTNG1 and VAV3) were located within it. The information of ROH islands, SNPs and genes harbored within the studied region were listed in Table 4. In addition, we found similar results on GGA5, in which the share ROH island was also strongly selected (Fig. S2).

Fig. 4
figure 4

Chromosome-wide distribution of selection signatures detected by iHS on Chromosome 8 for 5 Tibetan native chicken populations. The red line represents the threshold levels of SNPs with iHS value ranked top 0.1%. The green dots represented SNPs located in the studied region (0.03 Mb to 1.13 Mb) and other dots located outside of the studied region were shown in grey. SH, NM, DZ, LZ and LW denote Shigatse, Nyemo, Dagze, Ningychi and Lhasa white chicken population, respectively

Table 4 The extended homozygosity of ROH island detected on chromosome 8


This study used whole genome sequences of four Tibetan chicken populations reared on the Tibetan plateau to analyze genetic diversity, run of homozygosity, genomic inbreeding and selection signatures. A composite Tibetan local breed, Lhasa white was also included in the analyses to compare results among populations. Lhasa white is a synthetic breed generated by crossing male White Leghorn and female Tibetan chickens, which has been reared for over sixty years on the Tibetan plateau.

Observed (Ho) and expected (He) heterozygosity values for SNPs with MAF ≥ 0.05 close to 0.3 were found in all populations. Similar values were reported in modern chicken populations using sequence data [15] but lower values were documented for Italian local chickens [12]. Moreover, when all SNPs were used in the analysis to avoid bias [21], we found lower Ho for all populations compared to that reported in Dutch local chickens genotyped on 60 K SNP arrays [10]. In our study, we observed slightly lower heterozygosity than expected in SH, NM, DZ and LW, suggesting subtle inbreeding in these populations. However, a little heterozygosity excess (Ho > He) was observed in the LZ population. This suggests a recent bottleneck or an isolate-breaking effect [22] which may likely be due to the recent domestication and selection process, but also due to the small sample size. Pair-wise Fst among SH, NM and DZ was less than 0.05, indicating little genetic differentiation according to Wright’s interpretation [23]. The LZ population was isolated from SH, NM and DZ populations, suggesting little or no admixture between them during its domestication. This suggestion is supported by the clear separation of LZ birds from SH, NM and DZ revealed by principal component analysis, which corroborated the earlier submission that two or more distinct Tibetan chicken populations live in the plateau [2, 4].

In our present study, the number and the distribution of ROHs identified in Tibetan native chickens were comparable to those reported in broiler chickens [14]. Furthermore, short ROHs were predominant in the Tibetan chicken populations, indicating that little deleterious inbreeding happened in Tibetan native chicken populations [24]. The relationship between the total number of ROHs and the total length of the genome covered by ROHs showed considerable variation among birds within and across populations. Similar distributions were also reported in commercial chickens [15] and other livestock species including sheep [17] and cattle [25]. Genomic data is the only reliable source to estimate the inbreeding and relatedness of populations in the absence of other data sources, such as pedigree. The proportion of the genome within ROHs (FROH), genomic SNP-by-SNP inbreeding coefficient (FGRM), excess of homozygosity (FHOM) and correlation between uniting gametes (FUNI) were commonly accepted indicators for inbreeding assessment [26]. Herein, the ROH-based genomic inbreeding coefficients of Tibetan chicken were similar to the estimates in other Chinese local chickens studied [11], but much smaller than those in modern chickens [15, 16] and Italian local chickens [12]. Therefore, Tibetan native chickens are diverse and less affected by inbreeding, suggesting that they have maintained their natural diversity in the plateau. The correlations between FROH and FHOM, as well as FUNI were significantly high, similar to those reported in cattle [27], pigs [28], horses [29] and modern chickens [15]. This further supports FROH as an accurate estimate of identity by descent (IBD).

Whole-genome sequencing allows detection of the ROH and its analysis enables reliable inference of the demographic history of animal populations. It also provides a new approach for the exploration and discovery of complex traits [5]. The ROH islands of Tibetan chickens and Lhasa white chickens harbored many QTLs and candidate genes controlling economically important traits, including conformation, production, egg and meat quality, digestion and absorption, reproduction and growth traits. The QTLs for the comb weight and the ileum weight harbored within the common ROH islands in the Tibetan native chicken populations might play an important role in their adaptation to high altitude. Regarding common genes located on GGA5, leucine-rich repeat-containing G protein-coupled receptor 4 (LGR4), enriched in the biological process of osteoblast in GO database, contributes to the regulation of energy metabolism including food intake and energy expenditure [30]. Brain-derived neurotrophic factor (BDNF), enriched in the positive regulation of synapse assembly is considered important for the temperature perception in chickens [31]. In rats, BDNF administration in the paraventricular nucleus reduced energy intake and decreased body weight [32]. STAT1 and STAT4 are members of Janus kinase (JAK)-signal transducer and activators of transcription (STAT) pathway that plays critical roles in facilitating various cellular reactions to cellular stress including hypoxia, ultraviolet light and hyperosmolarity [33]. Moreover, these genes were previously identified as ROH islands-associated genes in Mexican highland chickens [20], suggesting their critical roles in adaptation to highland. Additionally, metal ion binding was enriched in 34 of the identified genes. Although the process of how metal ion binding affects animal’s physiology and production is rarely reported, some genes enriched in the term including VAV3, NOS2, COL3A1 and PRKD1 were putative candidate genes associated with highland adaptation [34], implying that the metal ion binding may be associated with highland adaptation.

ROH islands might be the representative genomic regions under natural and artificial selection [35]. The iHS approach appears to be powerful for detecting ongoing selection processes for which the target allele has a moderate to high frequency within a population [36]. Our iHS analysis revealed that the common genomic region with different ROH islands on GGA5 and GGA8 overlapped with putative selection signatures in SH, NM, LZ and LW populations, indicating ongoing selective forces. Commonly identified regions by both iHS and ROH analysis harbored AMY2A, NTNG1 and VAV3 genes on GGA8. AMY2A encodes a member of the alpha-amylase family of proteins, which is involved in carbohydrates and glycogen metabolism, affecting growth, carcass traits and feed intake efficiency in chicken [37]. The previous report showed that AMY2A was under selection for metabolism, energy availability and response to thermal stress in African chickens [38]. Similar to African village conditions, chicken feeding is mainly based on scavenging, household waste and some grain supplementation in the Tibetan plateau. Therefore, carbohydrate metabolism, energy generation and transport are important traits for feeding adaptation. NTNG1 is a responsible gene for axon and neurite growth [39]. It was also differentially expressed in a chicken hepatocellular cell line in response to stress [40]. VAV3 is a member of the VAV gene family that plays vital role as guanosine nucleotide exchange factors for Rho GTPases and signaling adaptors downstream of protein tyrosine kinases [41]. Specifically, VAV3 had been identified as a candidate gene associated with highland adaptation in Ethiopians [42] and Ethiopian sheep [43]. This probably resulted from its role played in the homeostasis of the cardiovascular system [44]. We therefore suggest that VAV3 also functions putatively in the adaptation to the high altitude of Tibetan native chicken. Given that different chicken populations were reared in the Tibetan plateau for many decades, we speculated that this genomic region and the candidate genes on GGA8 might be under natural or artificial selection for adaptation to the high-altitude environment.


In the present study, we used whole genome sequence data to characterize the genetic diversity and investigate the distribution of ROH across the genomes of five Tibetan indigenous chicken populations. Different LD, diversity and ROH patterns were observed in the five populations. Genetic diversity evaluated by observed heterozygosity was high in the five populations. The Nyingchi (LZ) population, which was distant from other populations had the highest proportion of long ROH fragments and ROH-based genomic inbreeding coefficient, reflecting recent inbreeding events. We identified a total of 343 ROH islands harboring 112 QTLs and 1429 genes. Five of such genes were involved in energy metabolism and STATs pathway. Specifically, ROH islands on GGA8 harbored genes including AMY2A, NTNG1 and VAV3. This region is suggested as a candidate genomic region for adaptation to the high-altitude environment, which should be further validated in the following studies. Our findings contribute to the understanding of genetic diversity, population inbreeding and the underlying genetic mechanism of the high-altitude adaptation, and may help in the design and implementation of breeding and conservation strategies for the chickens.


Ethics statement

All birds were handled following the guidelines established by the Council for Animal Welfare of China. The experimental protocols were approved by the Science Research Department of the Institute of Animal Sciences, Chinese Academy of Agricultural Sciences (CAAS).

Whole-genome sequencing and data processing

Genomic DNA was extracted from the blood of 114 female chickens. After purification and integrity verification of the DNA, whole-genome sequencing was performed using the Illumina HiSeq 2500 Sequencer (Illumina, Inc., SanDiego, CA, USA) to generate 150 bp paired-end reads. To minimize mapping errors, low-quality reads were removed using FastQC software following Yan et al. [45]. The clean reads from each bird were aligned to the chicken reference genome (Gallus gallus5.0) using the Burrows–Wheeler Alignment (BWA) tool [46] with the default parameters. Picard toolkit was subsequently used to filter out potential PCR duplicate reads. The resulting alignments were indexed using SAMtools [47] and processed according to the best practices for the Genome Analysis Toolkit (GATK) [48]. To obtain high-quality SNPs, we set a minimum quality score of 20 for both bases and mapped reads to call variants. Finally, the SNPs of each bird were combined to obtain a common set of SNP data, which was subjected to filtering with the following rigorous criteria using the GATK Variant Filtration module; (a) quality by depth > 5.0 (b) mapping quality score > 40.0 (c) FS < 60.0 (d) MQRankSum > − 12.5 (e) ReadPosRankSum > − 8.0 and (f) Filtering out any three SNPs clustered in a 10 bp window [49]. The qualified SNPs were annotated using the chicken reference genome. Finally, annotated SNP data was filtered using PLINK v1.90 software [50] with the following parameters: sample call rate > 90%, SNP call rate > 95%, and Hardy–Weinberg equilibrium p-value < 10− 5. After these quality control steps, a total of 20,385,015 SNPs distributed across 33 autosomes in 105 birds were obtained for subsequent analyses. Among these birds, 20, 32, 25, 10 and 18 belonged to Shigatse, Nyemo, Dagze, Nyingchi and Lhasa white population, respectively.

Genetic diversity analysis

The filtered SNPs were further pruned to obtain independent SNP markers using the ‘-indep-pairwise’ option, with a window size of 50 SNPs, a step of 10 SNPs, and r2 threshold of 0.2. Principal components (PC) analysis was carried out with the pruned data in PLINK, and the first 2 PCs were extracted and plotted using R software. The observed heterozygosity (Ho) and the expected heterozygosity (He) were estimated based on allele frequency of all eligible SNPs, SNPs with minor allele frequency (MAF) ≥0.01 and SNPs with MAF ≥ 0.05, respectively, using ‘-hardy’ option in PLINK. Multiple-locus heterozygosity (MLH) was also calculated for each bird with formula \(\frac{N-\mathrm{O}}{\mathrm{N}}\) [51], where O is the number of observed homozygotes, and N is the number of non-missing autosomal SNPs. The O and N were calculated using ‘-het’ option in PLINK. Further, PopLDdecay [52] was used to estimate linkage disequilibrium (LD) using SNPs with MAF ≥0.01, and the LD decay was plotted in R software. Moreover, fixation index (Fst) value of SNPs across the genomes were estimated using VCFtools [53], and then averaged to evaluate the relatedness among populations.

Runs of Homozygosity (ROH) detection

Prior to the ROH detection, the eligible SNPs with MAF ≥0.01 were separately filtered for each population. Long homozygous fragments were scanned in the pruned data using PLINK according to the following parameters: the minimum number of 50 SNPs in a ROH (−homozyg-snp 50), sliding windows of 50 SNPs (−homozyg-window-snp 50), allowance for not more than 5 missing SNPs (−homozyg-window-missing 5) and three heterozygous SNPs per window (−homozyg-window-het 3). The minimum length of an ROH segment was 300 kb (−homozyg-kb 300). The minimum SNP density was 1 SNP per 50 kb (−homozyg-density 50), and the maximum gap between two consecutive SNPs was 1000 kb (−homozyg-gap 1000). Finally, the rate in which a SNP was included in the total of sliding windows was at least 0.05 (−homozyg-window-threshold 0.05). After the run, the identified ROHs were classified according to their size into small (< 1 Mb), medium (1 ~ 3 Mb) and large (> 3 Mb) as previously delineated in chickens [10].

Genomic inbreeding analysis

Genomic inbreeding based on ROH (FROH) was estimated using PLINK according to previous methods [54]. The FROH for each bird was calculated as \(\frac{\sum_i{L}_{ROH_i}}{L_{aut}}\), where \({L}_{ROH_i}\) is the length of ROH of i th individual, and Laut is the genome length of autosome covered by the SNPs in the sequence data. The inbreeding coefficient for an individual based on homozygous SNPs (FHOM) was computed as \(\frac{\left(O-E\right)}{\left(L-E\right)}\), where O is the number of observed homozygotes, E is the number of expected homozygotes by chance, and L is the number of non-missing autosomal SNPs. Genomic SNP-by-SNP inbreeding coefficient (FGRM) and the correlation between uniting gametes (FUNI) were computed in GCTA software as previously reported [55]. Pair-wise correlations between these inbreeding coefficients were estimated by the Pearson method.

Identification of common ROH and gene annotation

To identify the genomic regions that harbored common ROHs across the five chicken populations, we estimated the occurrences of SNPs in ROHs by counting the number of times when the SNP was detected in those ROHs using detectRUN package [56] implemented in R. The genomic regions commonly associated with ROHs were screened by selecting the top 1% SNPs observed in ROHs. Adjacent SNPs that met this threshold were merged into genomic regions named ROH islands. Based on these consensus regions, we annotated QTL based on the chicken QTLdb using ‘-wa’ and ‘-wb’ options in BEDTools [57]. The ROH islands were also annotated with Gallus gallus5.0 genome assembly using the Ensembl BioMart [58] by extracting intersected and overlapped genes. Functional annotations and enrichments of the identified genes within the ROH islands were further carried out in the DAVID platform, and the chicken annotation file was set as background to identify significant (P < 0.05) GO terms and KEGG pathways.

Selection signatures analysis

To detect selection signatures in each ROH island, the integrated haplotype score (iHS) was calculated within each population. The iHS is a measure of the amount of extended haplotype homozygosity at a given SNP, that uses phased genotypes to identify putative regions of recent or ongoing positive selection in genomes [59]. Herein, the haplotype was phased using SHAPEIT [60] with recombination rate 0.01 as previously used for chicken genome [61]. The derived haplotypes were then analyzed using the rehh (v2.0) R package [62]. The iHS score was computed for each SNP and further standardized to a P value with the following formula piHS = −log10(1 − 2 | (iHS) − 0.5| ), where (iHS) represents the Gaussian cumulative distribution function, and piHS is the two sided P value associated with the neutral hypothesis of no selection [62]. The piHS higher than a threshold of 0.1% were considered as putative signatures of selection. Due to the limit of computing time and lack of accurate recombination rate of studied chickens, we only ran the selection signature analysis for candidate regions on GGA5 and GGA8 in the present study.

Availability of data and materials

The data and computing programs used in this manuscript are available from the corresponding authors on request.











Lhasa white


Run of homozygosity


ROH-based inbreeding coefficient


Linkage disequilibrium


Burrows–Wheeler Alignment


Genome Analysis Toolkit


Observed heterozygosity


Expected heterozygosity


Genomic SNP-by-SNP inbreeding coefficient


Excess of homozygosity


Correlation between uniting gametes


Minor allele frequency


Integrated haplotype score


Gene ontology

LGR4 :

Leucine-rich repeat-containing G protein-coupled receptor 4


Brain-derived neurotrophic factor


Amylase Alpha 2A


Netrin G1

VAV3 :

Vav Guanine Nucleotide Exchange Factor 3


  1. Wu C, Li N. Physiology and genetics of adaptability to highland hypoxia in Tibetan chicken. Beijing: China Agricultural University Press; 2012.

    Google Scholar 

  2. Wang M-S, Li Y, Peng M-S, Zhong L, Wang Z-J, Li Q-Y, et al. Genomic analyses reveal potential independent adaptation to high altitude in Tibetan chickens. Mol Biol Evol. 2015;32(7):1880–9.

    Article  CAS  PubMed  Google Scholar 

  3. Zhang Q, Gou W, Wang X, Zhang Y, Ma J, Zhang H, et al. Genome Resequencing identifies unique adaptations of Tibetan chickens to hypoxia and high-dose ultraviolet radiation in high-altitude environments. Genome Biol Evol. 2016;8(3):765–76.

    Article  PubMed  PubMed Central  Google Scholar 

  4. Li D, Che T, Chen B, Tian S, Zhou X, Zhang G, et al. Genomic data for 78 chickens from 14 populations. Gigascience. 2017;6(6):1–5.

  5. Ceballos FC, Joshi PK, Clark DW, Ramsay M, Wilson JF. Runs of homozygosity: windows into population history and trait architecture. Nat Rev Genet. 2018;19(4):220–34.

    Article  CAS  PubMed  Google Scholar 

  6. Peripolli E, Munari DP, Silva MVGB, Lima ALF, Irgang R, Baldi F. Runs of homozygosity: current knowledge and applications in livestock. Anim Genet. 2017;48(3):255–71.

    Article  CAS  PubMed  Google Scholar 

  7. Nothnagel M, Lu TT, Kayser M, Krawczak M. Genomic and geographic distribution of SNP-defined runs of homozygosity in Europeans. Hum Mol Genet. 2010;19(15):2927–35.

    Article  CAS  PubMed  Google Scholar 

  8. Qanbari S. On the extent of linkage disequilibrium in the genome of farm animals. Front Genet. 2020;10:1304.

    Article  PubMed  PubMed Central  CAS  Google Scholar 

  9. Gorssen W, Meyermans R, Janssens S, Buys N. A publicly available repository of ROH islands reveals signatures of selection in different livestock and pet species. Genet Sel Evol. 2021;53(1):2.

    Article  CAS  PubMed  PubMed Central  Google Scholar 

  10. Bortoluzzi C, Crooijmans RPMA, Bosse M, Hiemstra SJ, Groenen MAM, Megens H-J. The effects of recent changes in breeding preferences on maintaining traditional Dutch chicken genomic diversity. Heredity. 2018;121(6):564–78.

    Article  PubMed  PubMed Central  Google Scholar 

  11. Zhang M, Han W, Tang H, Li G, Zhang M, Xu R, et al. Genomic diversity dynamics in conserved chicken populations are revealed by genome-wide SNPs. BMC Genomics. 2018;19(1):598.

    Article  PubMed  PubMed Central  Google Scholar 

  12. Cendron F, Mastrangelo S, Tolone M, Perini F, Lasagna E, Cassandro M. Genome-wide analysis reveals the patterns of genetic diversity and population structure of 8 Italian local chicken breeds. Poult Sci. 2021;100(2):441–51.

    Article  CAS  PubMed  Google Scholar 

  13. Dementieva NV, Kudinov AA, Larkina TA, Mitrofanova OV, Dysin AP, Terletsky VP, et al. Genetic Variability in Local and Imported Germplasm Chicken Populations as Revealed by Analyzing Runs of Homozygosity. Animals. 2020;10(10):1887.

  14. Marchesi JAP, Buzanskas ME, Cantão ME, Ibelli AMG, Peixoto JO, Joaquim LB, et al. Relationship of runs of homozygosity with adaptive and production traits in a paternal broiler line. Animal. 2018;12(6):1126–34.

    Article  CAS  PubMed  Google Scholar 

  15. Talebi R, Szmatoła T, Mészáros G, Qanbari S. Runs of Homozygosity in Modern Chicken Revealed by Sequence Data. G3-Genes Genom Genet. 2020;10(12):4615–23.

    CAS  Google Scholar 

  16. Zhang J, Nie C, Li X, Ning Z, Chen Y, Jia Y, et al. Genome-wide population genetic analysis of commercial, indigenous, game, and wild chickens using 600K SNP microarray data. Front Genet. 2020;11:1154.

    Google Scholar 

  17. Mastrangelo S, Tolone M, Sardina MT, Sottile G, Sutera AM, Di Gerlando R, et al. Genome-wide scan for runs of homozygosity identifies potential candidate genes associated with local adaptation in Valle del Belice sheep. Genet Sel Evol. 2017;49(1):84.

    Article  PubMed  PubMed Central  CAS  Google Scholar 

  18. Malomane DK, Simianer H, Weigend A, Reimer C, Schmitt AO, Weigend S. The SYNBREED chicken diversity panel: a global resource to assess chicken diversity at high genomic resolution. BMC Genomics. 2019;20(1):345.

    Article  PubMed  PubMed Central  Google Scholar 

  19. Fleming DS, Koltes JE, Markey AD, Schmidt CJ, Ashwell CM, Rothschild MF, et al. Genomic analysis of Ugandan and Rwandan chicken ecotypes using a 600 k genotyping array. BMC Genomics. 2016;17(1):407.

    Article  CAS  PubMed  PubMed Central  Google Scholar 

  20. Strillacci MG, Vega-Murillo VE, Román-Ponce SI, López FJR, Cozzi MC, Gorla E, et al. Looking at genetic structure and selection signatures of the Mexican chicken population using single nucleotide polymorphism markers. Poult Sci. 2018;97(3):791–802.

    Article  CAS  PubMed  Google Scholar 

  21. Toro MA, Fernández J, Caballero A. Molecular characterization of breeds and its use in conservation. Livest Sci. 2009;120(3):174–95.

    Article  Google Scholar 

  22. Karamanlidis AA, Skrbinšek T, de Gabriel HM, Krambokoukis L, Munoz-Fuentes V, Bailey Z, et al. History-driven population structure and asymmetric gene flow in a recovering large carnivore at the rear-edge of its European range. Heredity. 2018;120(2):168–82.

    Article  CAS  PubMed  Google Scholar 

  23. Wright S. The interpretation of population structure by F-statistics with special regard to Systems of Mating. Evolution. 1965;19(3):395–420.

    Article  Google Scholar 

  24. Zhang Q, Guldbrandtsen B, Bosse M, Lund MS, Sahana G. Runs of homozygosity and distribution of functional variants in the cattle genome. BMC Genomics. 2015;16(1):542.

    Article  PubMed  PubMed Central  CAS  Google Scholar 

  25. Peripolli E, Stafuzza NB, Munari DP, Lima ALF, Irgang R, Machado MA. Panetto JCdC, Ventura RV, Baldi F, da Silva MVGB: assessment of runs of homozygosity islands and estimates of genomic inbreeding in Gyr (Bos indicus) dairy cattle. BMC Genomics. 2018;19(1):34.

    Article  PubMed  PubMed Central  CAS  Google Scholar 

  26. Alemu SW, Kadri NK, Harland C, Faux P, Charlier C, Caballero A, et al. An evaluation of inbreeding measures using a whole-genome sequenced cattle pedigree. Heredity. 2021;126(3):410–23.

    Article  CAS  PubMed  Google Scholar 

  27. Mastrangelo S, Tolone M, Di Gerlando R, Fontanesi L, Sardina MT, Portolano B. Genomic inbreeding estimation in small populations: evaluation of runs of homozygosity in three local dairy cattle breeds. Animal. 2016;10(5):746–54.

    Article  CAS  PubMed  Google Scholar 

  28. Schiavo G, Bovo S, Bertolini F, Tinarelli S, Dall’Olio S, Nanni Costa L, et al. Comparative evaluation of genomic inbreeding parameters in seven commercial and autochthonous pig breeds. Animal. 2020;14(5):910–20.

    Article  CAS  PubMed  Google Scholar 

  29. Polak G, Gurgul A, Jasielczuk I, Szmatoła T, Krupiński J, Bugno-Poniewierska M. Suitability of Pedigree Information and Genomic Methods for Analyzing Inbreeding of Polish Cold-Blooded Horses Covered by Conservation Programs. Genes. 2021;12(3):429.

  30. Li Z, Zhang W, Mulholland MW. LGR4 and its role in intestinal protection and energy metabolism. Front Endocrinol. 2015;6:131.

    Google Scholar 

  31. Lamont SJ, Coble DJ, Bjorkquist A, Rothschild MF, Persia M, Ashwell C, et al. Genomics of heat stress in chickens. Proceedings, 10th World Congress on Genetics Applied to Livestock Production. Vancouver; 2014.

  32. Wang C, Bomberg E, Billington C, Levine A, Kotz CM. Brain-derived neurotrophic factor in the hypothalamic paraventricular nucleus reduces energy intake. Am J Phys Regul Integr Comp Phys. 2007;293(3):R1003–12.

    CAS  Google Scholar 

  33. Dudley AC, Thomas D, Best J, Jenkins A. The STATs in cell stress-type responses. Cell Commun Signal. 2004;2(1):8.

    Article  PubMed  PubMed Central  CAS  Google Scholar 

  34. Witt KE, Huerta-Sánchez E. Convergent evolution in human and domesticate adaptation to high-altitude environments. Philos Trans R Soc Lond Ser B Biol Sci. 2019;374(1777):20180235.

    Article  CAS  Google Scholar 

  35. Liu J, Shi L, Li Y, Chen L, Garrick D, Wang L, et al. Estimates of genomic inbreeding and identification of candidate regions that differ between Chinese indigenous sheep breeds. J Anim Sci Biotechnol. 2021;12(1):95.

    Article  CAS  PubMed  PubMed Central  Google Scholar 

  36. Zhao F, McParland S, Kearney F, Du L, Berry DP. Detection of selection signatures in dairy and beef cattle using high-density genomic information. Genet Sel Evol. 2015;47:49.

    Article  PubMed  PubMed Central  CAS  Google Scholar 

  37. Zhang Z, Zhong H, Lin S, Liang L, Ye S, Xu Z, et al. Polymorphisms of AMY1A gene and their association with growth, carcass traits and feed intake efficiency in chickens. Genomics. 2021;113(2):583–94.

    Article  CAS  PubMed  Google Scholar 

  38. Elbeltagy AR, Bertolini F, Fleming DS, Van Goor A, Ashwell CM, Schmidt CJ, et al. Natural selection footprints among African chicken breeds and village ecotypes. Front Genet. 2019;10:376.

    Article  CAS  PubMed  PubMed Central  Google Scholar 

  39. Yamagishi S, Bando Y, Sato K. Involvement of Netrins and their receptors in neuronal migration in the cerebral cortex. Front Cell Dev Biol. 2021;8:1870.

    Article  Google Scholar 

  40. Sun L, Lamont SJ, Cooksey AM, McCarthy F, Tudor CO, Vijay-Shanker K, et al. Transcriptome response to heat stress in a chicken hepatocellular carcinoma cell line. Cell Stress Chaperon. 2015;20(6):939–50.

    Article  CAS  Google Scholar 

  41. Rodríguez-Fdez S, Bustelo XR. The Vav GEF Family: An Evolutionary and Functional Perspective. Cells. 2019;8(5):465.

  42. Scheinfeldt LB, Soi S, Thompson S, Ranciaro A, Woldemeskel D, Beggs W, et al. Genetic adaptation to high altitude in the Ethiopian highlands. Genome Biol. 2012;13(1):R1.

    Article  CAS  PubMed  PubMed Central  Google Scholar 

  43. Edea Z, Dadi H, Dessie T, Kim K-S. Genomic signatures of high-altitude adaptation in Ethiopian sheep populations. Genes Genom. 2019;41(8):973–81.

    Article  Google Scholar 

  44. Sauzeau V, Sevilla MA, Rivas-Elena JV, de Álava E, Montero MJ, López-Novoa JM, et al. Vav3 proto-oncogene deficiency leads to sympathetic hyperactivity and cardiovascular dysfunction. Nat Med. 2006;12(7):841–5.

    Article  CAS  PubMed  PubMed Central  Google Scholar 

  45. Yan L, Yang M, Guo H, Yang L, Wu J, Li R, et al. Single-cell RNA-Seq profiling of human preimplantation embryos and embryonic stem cells. Nat Struct Mol Biol. 2013;20(9):1131–9.

    Article  CAS  PubMed  Google Scholar 

  46. Li H, Durbin R. Fast and accurate short read alignment with burrows-wheeler transform. Bioinformatics. 2009;25(14):1754–60.

    Article  CAS  PubMed  PubMed Central  Google Scholar 

  47. Li H, Handsaker B, Wysoker A, Fennell T, Ruan J, Homer N, et al. Genome project data processing S: the sequence alignment/map format and SAMtools. Bioinformatics. 2009;25(16):2078–9.

    Article  PubMed  PubMed Central  CAS  Google Scholar 

  48. McKenna A, Hanna M, Banks E, Sivachenko A, Cibulskis K, Kernytsky A, et al. The genome analysis toolkit: a MapReduce framework for analyzing next-generation DNA sequencing data. Genome Res. 2010;20(9):1297–303.

    Article  CAS  PubMed  PubMed Central  Google Scholar 

  49. Li H, Ruan J, Durbin R. Mapping short DNA sequencing reads and calling variants using mapping quality scores. Genome Res. 2008;18(11):1851–8.

    Article  CAS  PubMed  PubMed Central  Google Scholar 

  50. Purcell S, Neale B, Todd-Brown K, Thomas L, Ferreira MAR, Bender D, et al. PLINK: a tool set for whole-genome association and population-based linkage analyses. Am J Hum Genet. 2007;81(3):559–75.

    Article  CAS  PubMed  PubMed Central  Google Scholar 

  51. Soh P, Hsu W, Khatkar MS, Williamson P. Evaluation of genetic diversity and management of disease in border collie dogs. Sci Rep. 2021;11(1):6243.

    Article  CAS  PubMed  PubMed Central  Google Scholar 

  52. Zhang C, Dong S-S, Xu J-Y, He W-M, Yang T-L. PopLDdecay: a fast and effective tool for linkage disequilibrium decay analysis based on variant call format files. Bioinformatics. 2019;35(10):1786–8.

    Article  CAS  PubMed  Google Scholar 

  53. Danecek P, Auton A, Abecasis G, Albers CA, Banks E, DePristo MA, et al. The variant call format and VCFtools. Bioinformatics. 2011;27(15):2156–8.

    Article  CAS  PubMed  PubMed Central  Google Scholar 

  54. McQuillan R, Leutenegger A-L, Abdel-Rahman R, Franklin CS, Pericic M, Barac-Lauc L, et al. Runs of Homozygosity in European populations. Am J Hum Genet. 2008;83(3):359–72.

    Article  CAS  PubMed  PubMed Central  Google Scholar 

  55. Yang J, Lee SH, Goddard ME, Visscher PM. GCTA: a tool for genome-wide complex trait analysis. Am J Hum Genet. 2011;88(1):76–82.

    Article  CAS  PubMed  PubMed Central  Google Scholar 

  56. Biscarini F, Cozzi P, Gaspa G. Marras G: detectRUNS: an R package to detect runs of homozygosity and heterozygosity in diploid genomes. Italy: University of Guelph; 2019.

    Google Scholar 

  57. Quinlan AR, Hall IM. BEDTools: a flexible suite of utilities for comparing genomic features. Bioinformatics. 2010;26(6):841–2.

    Article  CAS  PubMed  PubMed Central  Google Scholar 

  58. Smedley D, Haider S, Ballester B, Holland R, London D, Thorisson G, et al. BioMart – biological queries made easy. BMC Genomics. 2009;10(1):22.

    Article  PubMed  PubMed Central  CAS  Google Scholar 

  59. Voight BF, Kudaravalli S, Wen X, Pritchard JK. A map of recent positive selection in the human genome. PLoS Biol. 2006;4(3):e72.

    Article  PubMed  PubMed Central  Google Scholar 

  60. Delaneau O, Marchini J, Zagury J-F. A linear complexity phasing method for thousands of genomes. Nat Methods. 2012;9(2):179–81.

    Article  CAS  Google Scholar 

  61. Weng Z, Wolc A, Su H, Fernando RL, Dekkers JCM, Arango J, et al. Identification of recombination hotspots and quantitative trait loci for recombination rate in layer chickens. J Anim Sci Biotechnol. 2019;10(1):20.

    Article  PubMed  PubMed Central  Google Scholar 

  62. Gautier M, Klassmann A, Vitalis R. Rehh 2.0: a reimplementation of the R package rehh to detect positive selection from haplotype structure. Mol Ecol Resour. 2017;17(1):78–90.

    Article  CAS  PubMed  Google Scholar 

Download references


We thank Mmes. Xueying Ma and Jing Feng for providing Lhasa white chickens used in this study. We also sincerely thank Dr. Adamu Mani Isa of the Department of Animal Science, Usmanu Danfodiyo University, Sokoto for his careful revision of the manuscript.


This study was supported by Key R & D and transformation plan of Tibet Autonomous Region (CGZH2017000237), Agricultural Science and Technology Innovation Program (ASTIP-IAS04), National Key R & D Program (2021YFD1200302) and China Agriculture Research Systems (CARS-40).

Author information

Authors and Affiliations



JY contributed in conducting the study, the analysis of data and wrote the draft manuscript. ZS and SL contributed in the interpretation of result. SL, MZ, XL & ZY contributed in sample and data collection. SL, NY & JC designed the study, contributed in the interpretation of the result and also reviewed the manuscript. All authors read and approved the final manuscript.

Corresponding author

Correspondence to Jilan Chen.

Ethics declarations

Ethics approval and consent to participate

All birds were handled following the guidelines established by the Council for Animal Welfare of China. The experimental protocols were approved by the ethics committee of Science Research Department of the Institute of Animal Sciences, Chinese Academy of Agricultural Sciences.

Consent for publication

Not applicable.

Competing interests

The authors declare that they have no competing interests.

Additional information

Publisher’s Note

Springer Nature remains neutral with regard to jurisdictional claims in published maps and institutional affiliations.

Supplementary Information

Rights and permissions

Open Access This article is licensed under a Creative Commons Attribution 4.0 International License, which permits use, sharing, adaptation, distribution and reproduction in any medium or format, as long as you give appropriate credit to the original author(s) and the source, provide a link to the Creative Commons licence, and indicate if changes were made. The images or other third party material in this article are included in the article's Creative Commons licence, unless indicated otherwise in a credit line to the material. If material is not included in the article's Creative Commons licence and your intended use is not permitted by statutory regulation or exceeds the permitted use, you will need to obtain permission directly from the copyright holder. To view a copy of this licence, visit The Creative Commons Public Domain Dedication waiver ( applies to the data made available in this article, unless otherwise stated in a credit line to the data.

Reprints and permissions

About this article

Check for updates. Verify currency and authenticity via CrossMark

Cite this article

Yuan, J., Li, S., Sheng, Z. et al. Genome-wide run of homozygosity analysis reveals candidate genomic regions associated with environmental adaptations of Tibetan native chickens. BMC Genomics 23, 91 (2022).

Download citation

  • Received:

  • Accepted:

  • Published:

  • DOI: