Pigment regression due to genetic factors is commonly observed in vertebrates (1, 2). This can appear as a morphological trait of an entire population or species adapted to specific environments (e.g., darkness in caves or soil layers) (3, 4). More often, hypopigmentation is a minority feature appearing in members of a population or species carrying allele variants (5, 6). Typical examples are the pigmentary disorders occurring in humans, such as vitiligo, piebaldism, and albinism. This raises the question as to which internal factors regulate the pigment cells and pigmentation processes and thus govern hypopigmentation phenomena in vertebrates.

From a mechanistic perspective, hypopigmentation is always associated with mutations in genes that either participate in melanin synthesis (e.g., tyrosinase and melanosomal transmembrane protein) (4, 5) or are responsible for signal regulation in melanin synthesis and melanophore proliferation (e.g., melanocortin receptors and agouti) (3, 7). In the latter type of depigmented morphs, hypopigmentation is accompanied by additional physiological outcomes due to the potential crosstalk between regulation pathways (8, 9). Pigmentation is linked to the levels of many endocrine factors and the activation of related signal transduction pathways, including the melanocyte-stimulating hormones (MSHs), adrenocorticotropic hormone, steroid hormones, and prostaglandins (10). Among these factors, alpha-MSH is the best described; its precursor, proopiomelanocortin, is synthesized in the brain and pituitary primordium (11). Additionally, hypopigmentation can alternatively arise as a result of genetic variation in cellular processes indirectly related to melanin systems. One typical example of this type of case is vitiligo (12), an autoimmune skin disease characterized by patches of depigmentation caused by the destruction of melanophores (13). This disease is linked to genetic polymorphism of the MHC regions (14–17).

Depigmented animals are excellent models for investigation of the genetic determination of the corresponding morphological traits, and studies of such animals may provide evolutionary and mechanistic insights into animal and human chromatophore-related diseases (e.g., vitiligo, albinism, and melanoma). Although recent investigations have deepened our understanding of these diseases, many questions surrounding their initiation and progression remain to be answered (13), and proper animal models are important as a means of studying their pathogenesis and facilitating the discovery and evaluation of therapeutic interventions (18). The Chinese giant salamander (Cryptobranchidae: Andrias davidianus) is the largest extant amphibian species (19). Its ancestors diverged from other amphibians over 170 million years ago during the Jurassic Period (20), making it one of the oldest families on the amphibian tree of life. The artificial culture of A. davidianus is flourishing in China due to the great ecological and economic value of this animal. In practice, approximately 0.1% of individuals of the species express an albino phenotype of the entire body (Figure 1A), characterized by semitransparent skin with reduced skin pigmentation (Figures 1B, C). These albino individuals always suffer additionally from other physiological problems. For example, they grow more slowly than their typical siblings and do not live through their first year. The genetic basis of their albino phenotype has potential significance in the aquaculture breeding of A. davidianus and in the development of disease models. In this study, the gene transcriptional profiles of the dorsal and tail skin and brain are compared between albino and typical individuals to provide insight into the molecular basis of depigmentation and associated physiological abnormalities.

Figure 1 Morphological and transcriptional differences between A'lbino and typical A davidianus. (A, B) Appearance of albino individuals. (C) Comparison of dorsal skin pigment. (D, E) Differences between groups in overall length (D) and head width (E). (F) Transcriptome sampling scheme. (G) PCoA scatter plot presenting the similarity in transcriptional profiles between samples. Dorsal and tail skin tissues with the same labels are samples from the same individuals. ‘A’ and ‘B’ in the labels denote ‘Albino’ and ‘typical’, respectively. Note that dorsal and tail skins from the same individuals are always associated with more similar transcriptomes, indicating a high degree of technical stability in our RNA-seq data. (H) Statistics on DEGs in different tissues. (I) Transcriptional variations of genes in melanin synthesis. Asterisks denote significant differences between groups (at a threshold of p < 0.05, Student’s t test).

Artificially-bred A. davidianus were collected from a farm located in Hongya, Sichuan Province, China (103°10′05′′E, 29°52′36′′N). Fertilized eggs of this species take approximately 70 days to develop to the stage at which formation of the fifth digit in the hindlimb occurs (stage 46 according to the chronological table suggested by Shi and Boucaut (21)) (22). These larvae take more than 300 days to complete their metamorphosis (characterized by disappearance of the gills) after hatching at 15 °C. In this study, albino and typical individuals (n = 7 per group) were collected approximately 70–80 days after hatching. All the larvae were at stage 46 (21) when they were collected. The albino phenotype appeared sporadically across clutches, and thus the seven albino individuals all had different parents. Correspondingly, the typical individuals were randomly collected from different clutches. Their body weight, body length, head width, tail length, head length, orbital distance, nasal distance, and eye diameters were measured. After having been euthanized with MS-222, the larvae were dissected to collect the dorsal skin, the tail skin, and the brain. All samples were stored at -80 °C.

Transcriptome measurements were obtained from the brain and the dorsal and tail skin of six albino and three typical individuals. The brain transcriptome was included as the brain is a central organ in the regulation of pigmentation in amphibians. Previously-described protocols were followed for RNA extraction, purification, cDNA library construction, filtration, assembly, annotation, and gene expression quantification (23). RNA-seq was performed on an Illumina HiSeq 4000 platform by Annoroad (Beijing); paired-end reads were generated. A previously-reported multi-organ whole-length transcriptome was used as the reference genome for transcript identification, annotation, and gene expression quantification (24), and FPKM values were calculated for each unigene (see gene expression matrix in Table S1). This approach may improve the accuracy and reliability of the transcriptional quantification in comparison to a de novo assembly approach. Differently expressed genes (DEGs) were identified by Student’s t test, and functional enrichment analyses were conducted by querying DEGs against the KEGG database (based on KOBAS 3.0, with default parameters) (25).

Sequences of targeted genes were retrieved from Genbank or from our transcriptome database. Sequence alignment was performed using Clustal X2, and further edits were made using GeneDoc. A maximum likelihood tree was constructed using MEGA7 with default parameters.

3D models of MHC class I proteins were predicted on the SWISS-MODEL server (https://www.swissmodel.expasy.org/) with “marsupial MHC I (7edo.2.A)” as the model. Analysis of 3D models was performed using Swiss- Pdb Viewer.

Basic statistical analyses were conducted using IBM SPSS 22.0 (SPSS Inc., Chicago, USA). Kolmogorov–Smirnov tests were conducted to assess the deviation of the data from normal distribution. Inter-group differences in body length and width were analyzed via ANCOVA with body weight as a covariate. Inter-group differences in other body traits (using relative values) were analyzed via Student’s t test. Dissimilarity in transcriptomes was operationalized using the Bray–Curtis distance, which was calculated using the vegdist function of the R package Vegan. Subsequently, PERMANOVA (using the Adonis function of the Vegan package) was employed to identify differences in transcriptomes between the groups, and PCoA was conducted to present the differences between groups. Graphs were generated using GraphPad Prism 5 or the R package ggplot2 (26).

The albino individuals had comparable body width to typical individuals with similar body weight, but their bodies were smaller in length, resulting in a stockier overall body shape (Figures 1D, E, S1). Their relative tail length, head length, orbital distance, and eye diameters were comparable to those of typical individuals; however, their nasal distance was much smaller (Figure S1). Transcriptomics was conducted for the dorsal and tail skin samples, as well as the brain (Figure 1F). As expected, the two types of skin tissue shared similar transcriptional profiles, which were quite distinct from that of the brain (Figure 1G). PERMANOVA results suggested that there was a significant inter-group difference in brain transcriptome, while the transcriptional differences between groups in the skin samples were less significant (Figure S2). In total, 1,023, 1,183, and 840 DEGs between the groups were identified in the brain, dorsal skin, and tail skin, respectively, among which 299 DEGs were shared by the two skin tissue types (Figure 1H).

Three genes involved in melanin synthesis (namely, tyrosinase, tyrosinase B, and dopachrome tautomerase) were identified among the DEGs; these exhibited decreased transcription in the skin of albino individuals (Figure 1I). Interestingly, however, melanin synthesis did not account for the most prominent transcriptional changes in the albino larvae. Instead, an MHC class I gene was highlighted as exhibiting the highest degree of fold changes (in the albino vs. the typical group) in all three types of tissue (fold changes = 2,000, 144, and 340 in the brain, dorsal skin, and tail skin, respectively; Figures 2A, B). Correspondingly, KEGG enrichment analysis based on brain DEGs (at a threshold of p < 0.05) highlighted changes in immune-related processes (e.g., autophagy, endocytosis, and Fc gamma R-mediated phagocytosis), signal transduction (e.g., the AMPK signaling pathway and adrenergic signaling in cardiomyocytes), and metabolic pathways (e.g., thermogenesis and oxidative phosphorylation) (Figure 2C). The skin DEGs (meeting the threshold of p < 0.05 in both tissue types) highlighted changes in the autoimmune and cancer pathways, e.g., rheumatoid arthritis and bladder cancer (Figure 2D). Furthermore, several signaling pathways (e.g., adrenergic signaling in cardiomyocytes and the estrogen signaling pathway) were also enriched by DEGs, implying the potential presence of endocrine disorders in the albino larvae.

Figure 2 Identification of crucial DEGs and cellular processes. (A) The DEGs representing the greatest fold change in the brain. (B) Fold changes in DEGs (p < 0.05) in the dorsal (horizontal axis) and tail (vertical axis) skins. Red dots denote DEGs shared by all three tissue types; blue dots denote DEGs shared by the two skin tissue types. Note that most DEGs exhibit consistent trends in variation (i.e., an increase or decrease) between the two types of skin tissue, indicating the reliability of the DEGs. (C, D) Functional analysis based on DEGs in the brain (C) or those shared by the skin samples (D). Only the top 30 significant items are displayed (q < 0.01).

A total of eight MHC class I isotypes (labeled ISO 1–8) were identified in the whole-length transcriptome of A. davidianus (Figures 3A, S3). Isotypes ISO 1–3 and ISO 4–5 belong to different subclasses, whose members have been reported previously, while isotypes ISO 6–8 belong to new subclasses. ISO 1 is the dominating isotype expressed in the tissues of typical individuals. In albino individuals, the transcription of this isotype and β2 microglobulin (B2M), components of MHC I, exhibited dramatic downregulation. Meanwhile, they also exhibited notably increased transcription of ISO 2 and ISO 3 (Figure 3B). The sequence variances between ISO 1 and ISO 2–3 were mainly distributed in the α1 and α2 functional domains, which are responsible for antigen binding (Figures 3C, D). Within these domains, six neutral amino acid residues were replaced by charged ones in ISO 2/3. These substitutions increased the charges and reduced the amount of space in the antigen-binding pocket.

Figure 3 Sequence and expression analysis of MHC class I genes. (A) Maximum likelihood tree for MHC class I genes (see the sequences in Table S2). Red items denote MHC class I isotypes identified in this study; black items denote the known MHC class I genes of A davidianus and closely related species. (B) Transcriptional levels of MHC class I isotypes and B2M. Asterisks denote significant differences between groups (at a threshold of p < 0.05, Student’s t test). (C) Amino acid substitutions in the antigen- binding regions of MHC class I isotypes. Red and blue characters represent negative and positive charges. (D) 3D structure of MHC class I isotypes. The sites of amino acid substitution are highlighted by electron clouds (blue: positive residues; red: negative residues). Red arrows denote antigen- binding pockets. Note the reduced space and increased charges in the antigen- binding pockets of isotypes 2 and 3.

Our results suggest that the albino phenotype is associated with shorter body length and reduced relative nasal distance in A. davidianus. As the relative tail and head lengths are unchanged, the shorter body length can most likely be attributed to a systematic problem with growth. The biological significance of nasal distance has not been fully understood in animals, but it is a common morphological trait used in amphibian taxonomy, as it reflects the anatomic feature of the skull bone. These morphological changes may be explained by skeletal dysplasias, a class of disorders characterized by orthopedic complications and varying degrees of dwarfism or short stature (27). The morphological abnormalities observed in the albino phenotype in A. davidianus are not restricted to pigmentation problems. Similar morphological abnormalities have not been reported in white and albino Ambystoma mexicanum, which can grow and breed normally (28). The white phenotype of A. mexicanum is linked to genetic mutations preventing the differentiation and migration of all types of chromatophores, while the albino phenotype is due to the dysfunction of enzymes responsible for melanin synthesis (29, 30). Albino A. davidianus larvae do retain some melanophore and melanin throughout the entire body, particularly in the tail, whereas white and albino A. mexicanum are almost devoid of melanin. This suggests that different mechanisms should be expected to underlie the depigmentation phenomenon in A. davidianus and A. mexicanum.

The decreased transcription of tyrosinases and dopachrome tautomerase observed in the skin could explain the reduced pigmentation in albino A. davidianus (Figure 1). Interestingly, the transcriptional change in the brain was even more prominent than that occurring in the skin samples (Figure S2), although the latter contains the largest amount of melanin. In both the brain and the skin tissues, the transcription of an MHC class I gene accounted for the largest fold changes in albino individuals (Figures 2A, B). MHC class I molecules occur on all nucleated cells and present intracellular peptides to killer T cells in order to mediate cellular immunity, including the cytotoxic effect, phagocytosis, and inflammation. MHC class I gene variants are highly polymorphic, and their profiles and expression levels play an important role in autoimmune disorders, infectious diseases, and immunosurveillance (31, 32). In albino A. davidianus, the changes in the transcriptional profiles of MHC class I isotypes are consistent with the enrichment of DEGs in autophagy, phagocytosis, autoimmune, and cancer pathways in the brain and skin (Figures 2C, D). Functional enrichment based on DEGs also suggests the potential presence of endocrine disorders in albino individuals. However, this may be a result of autoimmune problems, despite the fact that we cannot exclude the potential role of the endocrine system in causing pigment regression, growth retardation, and increased mortality in albino A. davidianus. Overall, our results suggest that the albino phenotype encompasses far more than a pigmentation defect; rather, it represents a systematic problem associated with the immune system. This could explain the concomitant physiological defects, such as malformation, delayed growth, and mortality at early development (33), as immune activity can be linked to growth and development performance through both resource allocation (34) and crosstalk between regulation pathways (35).

The association between MHC class I and depigmentation has been well established in vitiligo, which is visible in the form of white spots and affects ∼1% of the world’s human population. Mounting strands of evidence support the theory that vitiligo is linked to genetic changes in MHC gene regions (14–17). High-risk MHC class I alleles can present many autoantigens, including some melanocyte proteins, and thus induce autoimmune responses against melanocytes (36). In addition to the change in antigen specificity produced by variations in coding regions, the altered expression of MHC class I and II genes, caused by genetic variation in adjacent transcriptional regulatory regions, plays an even more important role in the activation of autoimmune responses (37–39). Albino A. davidianus is characterized by an inverse transcriptional change to MHC I isotypes 1 and 2/3. The albino-specific MHC I isotypes 2 and 3 exhibit increased charges and decreased space in their antigen- binding pockets, implying a drastic change in their antigen specificity. Accordingly, this transcriptional shift in MHC I isotypes may explain the depigmentation occurring in albino individuals. This suggests that the albino A. davidianus, which has uncovered skin and is readily accessible due to large-scale farming, constitutes a potential animal model for study of the pathogenesis of and therapeutic interventions for vitiligo or other autoimmune diseases.

On this basis, further investigations may focus on the molecular mechanisms underlying the associations between MHC I isotypes and morphological and physiological defects: for example, examining how the expression of MHC I isotypes 2 and 3 affects the autoimmune response and apoptosis in different tissue and cell types. Moreover, it is also important to study the genetic structure, polymorphism, and diversity of the MHC regions of A. davidianus at the population level, and to identify the regulatory regions that determine the expression profile of different MHC I isotypes.

The datasets presented in this study can be found in online repositories. The names of the repository/repositories and accession number(s) can be found in the article/Supplementary Material.

The Supplementary Material for this article can be found online at: https://www.frontiersin.org/articles/10.3389/fendo.2022.1053732/full#supplementary-material

1. Yamaguchi Y, Hearing VJ. Melanocytes and their diseases. Cold Spring Harbor Perspect Med (2014) 4(5):a017046. doi: 10.1101/cshperspect.a017046

CrossRef Full Text | Google Scholar

2. Jeffery WR, Ma L, Parkhurst A, Bilandžija H. Chapter 8 - pigment regression and albinism in astyanax cavefish. In: Keene AC, Yoshizawa M, McGaugh SE, editors. Biology and evolution of the Mexican cavefish. Oval Road, London NW1 7DX, England: Academic Press. (2016). p. 155–73.

Google Scholar

3. Zhu W, Liu L-S, Wang X-G, Gao X-Y, Jiang J-P, Wang B. Transcriptomics reveals the molecular processes of light-induced rapid darkening of the non-obligate cave dweller Oreolalax rhodostigmatus (Megophryidae, anura) and their genetic basis of pigmentation strategy. BMC Genomics (2018) 19(1):422. doi: 10.1186/s12864-018-4790-y

PubMed Abstract | CrossRef Full Text | Google Scholar

4. Yang J, Chen X, Bai J, Fang D, Qiu Y, Jiang W, et al. The sinocyclocheilus cavefish genome provides insights into cave adaptation. BMC Biol (2016) 14(1):1. doi: 10.1186/s12915-015-0223-4

PubMed Abstract | CrossRef Full Text | Google Scholar

6. Utzeri VJ, Bertolini F, Ribani A, Schiavo G, Dall’Olio S, Fontanesi L. The albinism of the feral asinara white donkeys (Equus asinus) is determined by a missense mutation in a highly conserved position of the tyrosinase (TYR) gene deduced protein. Anim Genet (2016) 47(1):120–4. doi: 10.1111/age.12386

PubMed Abstract | CrossRef Full Text | Google Scholar

8. Voisey J, Daal AV. Agouti: from mouse to man, from skin to fat. Pigment Cell Res (2010) 15(1):10–8. doi: 10.1034/j.1600-0749.2002.00039.x

CrossRef Full Text | Google Scholar

10. Abdel-Malek ZA. Endocrine factros as effectors of integumental pigmentation. Dermatol Clinics (1988) 6(2):175–83. doi: 10.1016/S0733-8635(18)30664-8

CrossRef Full Text | Google Scholar

11. Hayes WP, Loh YP. Correlated onset and patterning of proopiomelanocortin gene expression in embryonic xenopus brain and pituitary. Development (1990) 110(3):747–57. doi: 10.1242/dev.110.3.747

PubMed Abstract | CrossRef Full Text | Google Scholar

12. Shi F-Y, Kong B-W, Song JJ, Lee JY, Dienglewicz RL, Erf GF. Understanding mechanisms of vitiligo development in smyth line of chickens by transcriptomic microarray analysis of evolving autoimmune lesions. BMC Immunol (2012) 13(1):18. doi: 10.1186/1471-2172-13-18

PubMed Abstract | CrossRef Full Text | Google Scholar

14. Jin Y, Birlea SA, Fain PR, Gowan K, Spritz RA. Variant of TYR and autoimmunity susceptibility loci in generalized vitiligo. New Engl J Med (2010) 362(18):1686–97. doi: 10.1056/NEJMoa0908547

PubMed Abstract | CrossRef Full Text | Google Scholar

15. Jin Y, Andersen G, Yorgov D, Ferrara TM, Ben S, Brownson KM, et al. Genome-wide association studies of autoimmune vitiligo identify 23 new risk loci and highlight key pathways and regulatory variants. Nat Genet (2016) 48(11):1418–24. doi: 10.1038/ng.3680

PubMed Abstract | CrossRef Full Text | Google Scholar

17. Cao L, Zhang R, Wang Y, Hu X, Yong L, Li B, et al. Fine mapping analysis of the MHC region to identify variants associated with Chinese vitiligo and SLE and association across these diseases. Front Immunol (2021) 12:758652. doi: 10.3389/fimmu.2021.758652

PubMed Abstract | CrossRef Full Text | Google Scholar

18. Singh VP, Motiani RK, Singh A, Malik G, Aggarwal R, Pratap K, et al. Water buffalo (Bubalus bubalis) as a spontaneous animal model of vitiligo. Pigment Cell melanoma Res (2016) 29(4):465–9. doi: 10.1111/pcmr.12485

PubMed Abstract | CrossRef Full Text | Google Scholar

19. Fei L, Ye CY, Jiang JP. Colored atlas of Chinese amphibians and their distributions. Chengdu, China: Sichuan Publishing House of Science & Technology (2012).

Google Scholar

21. Shi D-L, Boucaut J-C. The chronological development of the urodele amphibian Pleurodeles waltl (Michah). Int J Dev Bio (1995) 39:427–41.

Google Scholar

22. Fan W-B, Zhu W-B, Zhang M-H, Zhao T, Jiang J-P. Preliminary studies on the early development of Andrias davidianus. Sichuan J Zool (2022) 41(5):517–25. doi: 10.11984/j.issn.1000-7083.20210464

CrossRef Full Text | Google Scholar

23. Zhu W, Zhao C-L, Feng J-Y, Chang J, Zhu W-B, Chang L-M, et al. Effects of habitat river microbiome on the symbiotic microbiota and multi-organ gene expression of captive-bred Chinese giant salamander. Front Microbiol (2022) 13:884880. doi: 10.3389/fmicb.2022.884880

PubMed Abstract | CrossRef Full Text | Google Scholar

24. Zhu W, Zhao C, Zhao T, Chang L, Chen Q, Liu J, et al. Rising floor and dropping ceiling: organ heterogeneity in response to cold acclimation of the largest extant amphibian. Proc R Soc B (2022) 289:20221394. doi: 10.1098/rspb.2022.1394

PubMed Abstract | CrossRef Full Text | Google Scholar

25. Xie C, Mao X, Huang J, Ding Y, Wu J, Dong S, et al. KOBAS 2.0: a web server for annotation and identification of enriched pathways and diseases. Nucleic Acids Res 39(Web Server issue) (2011) 39:W316–22. doi: 10.1093/nar/gkr483

CrossRef Full Text | Google Scholar

26. Wickham H. Ggplot2: Elegant graphics for data analysis. New York, NY, USA: Springer Publishing Company, Incorporated (2009).

Google Scholar

28. Frost SK, Epp LG, Robinson SJ. The pigmentary system of developing axolotls. III. an analysis of the albino phenotype. J Embryol Exp Morphol (1986) 92:255–68. doi: 10.1242/dev.92.1.255

PubMed Abstract | CrossRef Full Text | Google Scholar

30. Woodcock MR, Vaughn-Wolfe J, Elias A, Kump DK, Kendall KD, Timoshevskaya N, et al. Identification of mutant genes and introgressed tiger salamander DNA in the laboratory axolotl, ambystoma mexicanum. Sci Rep (2017) 7(1):6. doi: 10.1038/s41598-017-00059-1

PubMed Abstract | CrossRef Full Text | Google Scholar

31. Matzaraki V, Kumar V, Wijmenga C, Zhernakova A. The MHC locus and genetic susceptibility to autoimmune and infectious diseases. Genome Biol (2017) 18(1):76. doi: 10.1186/s13059-017-1207-1

PubMed Abstract | CrossRef Full Text | Google Scholar

34. Soler JJ, L.d. N, Pérez–Contreras T, Soler M, Sorci G. Trade-off between immunocompetence and growth in magpies: an experimental study. Proc R Soc London Ser B: Biol Sci (2003) 270(1512):241–8. doi: 10.1098/rspb.2002.2217

CrossRef Full Text | Google Scholar

36. Jin Y, Ferrara T, Gowan K, Holcomb C, Rastrou M, Erlich HA, et al. Next-generation DNA re-sequencing identifies common variants of TYR and HLA-a that modulate the risk of generalized vitiligo via antigen presentation. J Invest Dermatol (2012) 132(6):1730–3. doi: 10.1038/jid.2012.37

PubMed Abstract | CrossRef Full Text | Google Scholar

37. Hayashi M, Jin Y, Yorgov D, Santorico SA, Hagman J, Ferrara TM, et al. Autoimmune vitiligo is associated with gain-of-function by a transcriptional regulator that elevates expression of HLA-A*02:01 in vivo. Proc Natl Acad Sci USA (2016) 113(5):1357–62. doi: 10.1073/pnas.1525001113

PubMed Abstract | CrossRef Full Text | Google Scholar

38. Cavalli G, Hayashi M, Jin Y, Yorgov D, Santorico SA, Holcomb C, et al. MHC class II super-enhancer increases surface expression of HLA-DR and HLA-DQ and affects cytokine production in autoimmune vitiligo. Proc Natl Acad Sci USA (2016) 113(5):1363–8. doi: 10.1073/pnas.1523482113

PubMed Abstract | CrossRef Full Text | Google Scholar

39. Jin Y, Roberts GHL, Ferrara TM, Ben S, van Geel N, Wolkerstorfer A, et al. Early-onset autoimmune vitiligo associated with an enhancer variant haplotype that upregulates class II HLA expression. Nat Commun (2019) 10(1):391. doi: 10.1038/s41467-019-08337-4

PubMed Abstract | CrossRef Full Text | Google Scholar