Edited by Scott V. Edwards, Harvard University, Cambridge, MA, and approved December 17, 2018 (received for review August 7, 2018)
It is a tenet of modern biology that species adapt through natural selection to cope with the ever-changing environment. By comparing genetic variants between the island and mainland populations of a passerine, we inferred the related age of genetic variants across its entire genome and suggest that preexisting standing variants played the predominant role in local adaptation. Our findings not only resolve a long-standing fundamental problem in biology regarding the genetic sources of adaptation, but imply that the evolutionary potential of a population is highly associated with its preexisting genetic variation.
What kind of genetic variation contributes the most to adaptation is a fundamental question in evolutionary biology. By resequencing genomes of 80 individuals, we inferred the origin of genomic variants associated with a complex adaptive syndrome involving multiple quantitative traits, namely, adaptation between high and low altitudes, in the vinous-throated parrotbill (Sinosuthora webbiana) in Taiwan. By comparing these variants with those in the Asian mainland population, we revealed standing variation in 24 noncoding genomic regions to be the predominant genetic source of adaptation. Parrotbills at both high and low altitudes exhibited signatures of recent selection, suggesting that not only the front but also the trailing edges of postglacial expanding populations could be subjected to environmental stresses. This study verifies and quantifies the importance of standing variation in adaptation in a cohort of genes, illustrating that the evolutionary potential of a population depends significantly on its preexisting genetic diversity. These findings provide important context for understanding adaptation and conservation of species in the Anthropocene.
Adaptation lies at the heart of Darwinian evolution. It requires the presence of advantageous alleles that are either new mutations (1, 2) or preexisting standing variants (3, 4). Population genetic theory conventionally considers novel mutations to be the genetic source of adaptation (1, 5), and their significance has been highlighted in various empirical studies (6, 7). However, a high level of standing variation may allow a faster response to environmental changes than waiting for appropriate mutations to arise. Recent genomic studies also provide evidence for the role of standing variation in adaptation (8, 9). Nevertheless, the relative contribution of new mutations and standing variation to adaption has rarely been evaluated. This is fundamental to the identification of the major driver underpinning this critical evolutionary process and to the forecast of the fate of species in the Anthropocene, as the ability to adapt is central to species’ survival in changing environments (10).
Evaluating the relative contribution of standing and new genetic variation to adaptation requires assessing a cohort of genes that are under directional selection. Screening for the genetic basis underlying a complex adaptive syndrome is likely to yield a set of such genes. Here, we screened for genetic variants that underlie adaptation between high and low altitudes that likely involves many quantitative traits (11) of the vinous-throated parrotbill (Sinosuthora webbiana), one of the most widely distributed nonmigrant passerines in East Asia (12). In Taiwan, this bird is found from sea level up to 3,100 m (13). Individual parrotbills only disperse within a limited distance (14). Field data suggests that it exhibits no seasonal altitudinal migration and may only range within 100 m in altitude (15).
To evaluate the relative contribution of different sources of adaptive genetic variants, we inferred the relative age of genetic variants across the entire genome. As a continental island, Taiwan was periodically connected with the Asian mainland when land bridges surfaced during the Pleistocene glacial periods (16), allowing Taiwan and mainland populations to resume gene flow. Consequently, genetic variants shared by the two populations are likely to represent standing variation inherited from their most recent common ancestors. In contrast, new mutations which have arisen since the latest episode of gene flow would only be found in either one of the populations.
Since the last glacial period, temperatures in Taiwan have risen by about 9.6 °C (17), likely driving species to shift or expand their ranges from lowland glacial refugia toward higher altitudes (18, 19). While populations that expanded to higher altitudes faced challenges such as lower temperature and oxygen partial pressure, those that remained in the lowlands would also have found themselves in a warmer environment, an aspect that has been little investigated. Knowledge of how species adapted to global warming in the past may help us understand the challenges they face in a continuously warming planet.
We first sequenced and assembled de novo a reference genome of a male Taiwanese vinous-throated parrotbill. The assembled genome, 1.06 Gb in length (genomic coverage = 110×, SI Appendix, Table S1), contains 6,508 scaffolds (2,034 scaffolds > 10 kb) with N50 size of 1.94 Mb. This draft genome includes a total of 72 Mb (6.8%) of repetitive sequences (SI Appendix, Table S2). We found 20,149 predicted genes, 9,581 of which could be functionally annotated. In total, we obtained 278 Mb of genes (26% of total genome) consisting of 253 Mb of introns and 25 Mb of coding regions.
We resequenced the genomes of 40 parrotbill individuals that were collected from four populations (10 individuals per population; Fig. 1A) with an average depth of 5.8× genome coverage (SI Appendix, Table S3). They formed two high- (>2,000 m) and low- (<100 m) altitudinal comparison pairs, one from the east side of the Central Mountain Range (CMR) in Taiwan and one from the west side. The altitudinal gradients between the two high- and low-altitude populations are steep (43 km and 55 km apart for the east and west population pairs, respectively). The two low-altitude localities were much warmer (over 10 °C higher annual mean temperature) and drier (50% lower annual precipitation) than the two high-altitude localities (SI Appendix, Table S3).
(A) The vinous-throated parrotbill and four sites (red dots) on an east–west section of central Taiwan at which vinous-throated parrotbills were sampled and (B) the distribution of the FST and ΔFST in the east and west high-/low-altitude local population pairs of each 10-kb nonoverlapping genomic window that was aligned with the published genome of the zebra finch. Red dots on top of and within each panel represent candidate regions on the genome (n = 24); red horizontal lines indicate the top 1% of FST and ΔFST. EH, high-altitude population east of CMR; EL, low-altitude population east of CMR; WH, high-altitude population west of CMR; WL, low-altitude population west of CMR.
Comparing with the reference genome, 20,017,403 polymorphic sites (S) within the 40 resequenced genomes were found. The levels of intrapopulation genetic diversity are similar (π and mean individual heterozygosity,
The population structure was revealed in an unrooted neighbor-joining population tree (SI Appendix, Fig. S3) based on genome-wide FST (fixation index, a genetic measurement of population differential) values among all four populations (SI Appendix, Table S6), which indicates that pairs of populations on the same side of the CMR are clustered together. Thus, the two population pairs can serve as independent samples for illustrating parallel adaptation along altitudinal gradients (22).
Bottleneck and other events could leave genetic signatures resembling those of selection (23). We used the diploid sequence of the reference genome and a high-coverage genome sequence (24.5×) of an additional male Taiwanese parrotbill to test for recent severe bottleneck events. Results of pairwise sequentially Markovian coalescent (PSMC) analysis (24) indicate that their Ne varied through time, although without evidence of a severe population bottleneck (SI Appendix, Fig. S4). Ne was ∼30,000–100,000 from the middle of the last glacial period (∼50–70 thousand years ago) to the present.
For the east population pair, we found 2,812 outlier regions (0.715 ≥ FST ≥ 0.271; the top 1% of FST values; mean genome-wide FST = 0.108, 95% CI = 0.108–0.109); for the west population pair, we found 3,051 outlier pairs (0.536 ≥ FST ≥ 0.140; mean genome-wide FST = 0.052, 95% CI = 0.051–0.052) (Fig. 1B). Of these regions, 152 are shared between the east and west populations. It is important to consider that both ecologically divergent selection and purifying selection would reduce genetic diversity on functionally conserved regions and lead to high FST in the linked genomic region. Whereas divergent selection should only affect populations in different environments, by contrast, purifying selection would affect all populations (25). For the Taiwanese parrotbills, FST of 10-kb windows was highly correlated between all population pairs (SI Appendix, Fig. S5), implying that the high-FST regions across genome is largely the result of purifying selection. Purifying selection on linked genomic regions is further supported by the negative correlation between FST and π values (SI Appendix, Fig. S6), positive correlation between FST values and the level of linkage disequilibrium, r2 (SI Appendix, Fig. S7), and negative correlation between π and r2 values (SI Appendix, Fig. S8). Therefore, we employed ΔFST analysis to reduce the effect of purifying selection on genetic differentiation of genomic regions.
Only 24 10-kb nonoverlapping autosomal regions in the top 1% for both FST and ΔFST values (FST ≥ 0.464 and 0.344 in the east and west population pairs, respectively; and ΔFST ≥ 0.164 and 0.101 in the east and west population pairs, respectively) were considered candidate regions for adaptation to altitudes along both sides of the CMR (Fig. 1B). These regions have high levels of genetic differentiation between different altitudes and are less likely to be affected by purifying selection. Only six of the candidate regions are concatenated into two longer divergent genomic islands that were 20 and 40 kb in length.
Seventeen candidate genes were identified (SI Appendix, Table S8). Among them, four are involved in oxygen cascade: VAV3 and COL15A1 are related to angiogenesis; IGF2 to respiratory system phenotype; and TPPP to hemoglobin content. Three candidate genes are related to aspects of metabolism that could be involved in thermoregulation: SUPT7L (26) and HBP1 are related to lipid metabolism; and OLA1 to ATP hydrolysis. High‐altitude mammals have physiological adaptations in their circulatory, hematologic, respiratory, and thermoregulatory systems as well as in their metabolism in comparison with their sea-level counterparts (27, 28). Therefore, the genes we identified may be associated with adaptation between high and low altitudes. Among them, VAV3 (29), TPPP (30), IGF2 (31), and COL1A1 (32) have been reported to be associated with altitudinal adaption or hypoxia syndrome in other organisms. We also found that HBP1 (33) and OLA1 (34) are associated with temperature acclimation and heat shock phenotypes, respectively, in other vertebrates. The phenogenotypic associations of these genes in other systems not only validate our inferences, but also support the view that multiple quantitative traits, including oxygen cascade and thermal regulation, were involved in the adaptation to different altitudes.
Although adaptation could have been achieved through regulatory evolution of individual genes (e.g., those mentioned above), functionally related genes could also have changed expression simultaneously in the adaptive process. To examine this possibility, we performed an enrichment analysis on gene ontology (GO) terms. Among the 23 GO terms in which the 17 candidate genes were enriched [false discovering rate (FDR) < 0.05; P < 0.0019], many were found to be associated with glycogen metabolism (7 GO terms) and epigenetic control of gene expression (11 GO terms), especially histone acetylation (8 GO terms) for the latter group (SI Appendix, Table S9). This observation is consistent with previous studies showing reduced activity of glycogen metabolism accompanied by increased metabolism of fatty acids in organisms moving to a higher altitude (35, 36). This result also implies that epigenetic control, especially histone modification involving acetylation, played a critical role in the adaptation of vinous-throated parrotbills to different altitudes and their ability to cope with rapid environmental changes (37). Consistently with that possibility, the involvement of chromatin modifications in high-altitude adaptation has been observed in several mammals (38, 39).
Within the 24 candidate regions, 35 (average FST = 0.767; range: 0.631–1.000) and 32 (average FST = 0.628; range: 0.461–0.832) SNPs with the highest FST values in the east and west altitudinal pairs, respectively, were considered as candidate SNPs (SI Appendix, Tables S10–S13). The FST values of these candidate SNPs are significantly higher than those for all SNPs in both altitudinal pairs (Welch two-sample t test: t = 19.795 in the east pair and 21.592 in the west pair, P < 1 × 10−15). Furthermore, we found that most of the candidate SNPs are fixed or nearly fixed (frequency of major allele ≥ 0.9) in one of the four populations (34 and 27 of the candidate SNPs for the eastern and western populations, respectively; SI Appendix, Tables S10–S13), suggesting that they may have had major effects on the associated phenotypes. The higher polymorphism levels of the other candidate SNPs suggest that they are probably under soft sweep or polygenic selection and thus unable to drive these advantageous genetic variants to fixation (40). Noticeably, among these candidate SNPs, only 10 (29.4%) and 4 (14.8%) of them were fixed or nearly fixed in both high and low altitudes in the east and west populations, respectively. This implies that the intensity of divergent selection could be asymmetric at different altitudes for these candidate SNPs.
Intriguingly, none of these candidate SNPs is located within coding regions: eight are located in intronic regions of the HBP1 and AAK1 genes, and the other 59 are located in intergenic regions. The disproportionate number of candidate SNPs located in intergenic regions (Fisher’s exact test: P = 0.008) supports the predominant role of regulatory regions in adaptation (41, 42). Alleles in regulatory regions are often codominant in their effect, have mild pleiotropic consequences, and modify only the expression of genes rather than their functions (43).
The prerequisite for inferring standing variation using the comparative approach is recent coancestry between sampled taxa. We found that the mainland and Taiwan parrotbill populations diverged early in the last glacial period with estimated Nes of 482,065 for Taiwan and 1,527,174 for the mainland (SI Appendix, Table S14). Their recent divergence is also supported by the low genome-wide net genetic distance (Da) between them, 0.0007 (Dxy = 0.0051 between Taiwan and mainland populations; πs of Taiwan and mainland populations are 0.0040 and 0.0048, respectively). This implies a divergence time (t = Da/2μt; t, divergence time; μ, mutations per site per year) of about 152,000 y, meaning that genetic variants shared between the extant Taiwan and mainland parrotbill populations are unlikely to have arisen independently. However, the FST values between the mainland population and the low-altitude populations in Taiwan are lower (mean = 0.122) than those between the mainland and the high-altitude ones (mean = 0.133). This supports the view that the high-altitude populations are the leading edge of a postglacial population that possesses lower genetic diversity and a smaller Ne than the low-altitude refuge population.
Assuming genetic panmixia across the vinous-throated parrotbill’s mainland range (44), we further resequenced the genomes of 40 parrotbill individuals (SI Appendix, Table S4) to represent the species’ genetic diversity in the mainland. We found ∼1.3 times as many polymorphic sites in the mainland (S = 27,879,100) as in the Taiwan population (S = 20,017,403). About 66.2% of polymorphic sites in Taiwan (S = 13,259,879 sites) were also polymorphic in the mainland. These shared polymorphic sites likely represent common standing SNPs (minimum allele frequency, MAF, of SNPs ≥ 1.25%) in both the mainland and Taiwan populations, while the private genetic variants of each population are more likely to be rare standing SNPs (MAF < 1.25%) that our sampling scheme failed to detect in the other population, together with some new mutations that occurred since the populations split. About 94% of the candidate SNPs in the Taiwan populations are also polymorphic (or shared) in the mainland populations (Fig. 2A), significantly more than the 66.2% of all Taiwan SNPs (Fisher’s exact test: P = 1.6 × 10−5 and P = 0.002 for the east and west altitudinal pairs, respectively; Fig. 2A). The genomic regions around functional elements including genes should experience different evolutionary trajectories from those around nonfunctional regions because of the distinct densities of selection targets within them (45). Therefore, we compared the proportion of shared SNPs in a range 1 Mb upstream and downstream of all genes in the genome with the proportion of the candidate SNPs themselves that were shared. The proportion of candidate SNPs that was shared was still significantly higher (Fisher’s exact test: P < 0.001 and P = 0.004; Fig. 2A), confirming that shared standing genetic variation was the predominant genetic source of adaptation.
(A) The proportion of polymorphic SNPs shared with the mainland population is significantly higher for candidate SNPs inferred in both high-and low-altitude population pairs than that of the entire genome and SNPs within 1 Mb downstream and upstream regions of all genes in the genome (within 1 mb) of the Taiwan population (East, east population pair; West, west population pair; Fisher’s exact test, candidate SNPs: P = 1.6 × 10−5 and P = 0.002 for the east and west high-/low-altitude population pairs, respectively; all of the SNPs within 1 Mb regions: P < 0.00001 for both east and west high-/low-altitude population pairs). (B) The minimum allele frequency (MAF) of candidate SNPs in the mainland population is significantly higher than that of noncandidate SNPs (Welch two-sample t test: t = −10.524, df = 63, P = 1.644 × 10−15). (C) The MAF of shared variants is significantly higher than private variants (Welch two-sample t test, t = −2,909.8, df = 19,750,000, P < 2.2e-16). PV, private variants; SV, shared variants.
In the mainland parrotbill population, the MAF of candidate SNPs (average MAF = 0.300, σ2 = 0.021) was significantly higher than that of other SNPs (average MAF = 0.107, σ2 = 0.016; Welch two-sample t test: t = −10.524, df = 63, P = 1.644 × 10−15; Fig. 2B). Assuming that the mainland population with its large long-term Ne (SI Appendix, Table S14) is more likely to retain the ancestral allelic frequencies, our results support the view that standing variants with high initial frequencies facilitate a swift response to selection posed by changing environments (46, 47). Because the Ne of parrotbills is relatively large in Taiwan, genetic drift should only have a limited effect on the dynamics of allelic frequency in the population. Because the average MAF of private SNPs in Taiwan (average MAF = 0.038, σ2 = 0.002) is significantly lower than that of shared ones (average MAF = 0.158, σ2 = 0.019; Welch two-sample t test, t = −2909.8, df = 19,750,000, P value < 2.2e-16; Fig. 2C), these private SNPs are likely to have preexisted in the ancestral population in low frequency. The low proportion of private SNPs in candidate SNPs implies that rare genetic variants in the ancestral population are less likely to enable response to rapid environmental changes as predicted by the theory (4).
The significant role of shared standing genetic variation in adaptation revealed in this study (Fig. 2) has important implications for our understanding of adaptation and biological conservation. Our results suggest that to cope with environmental changes quickly, a population may benefit significantly from the short waiting time needed to exploit standing genetic variation, especially for variants with high initial frequencies. Our results are consistent with the accumulating evidence that species with high levels of genetic diversity can persist in changing environments (48) and that they are more likely to colonize novel environments (49). Therefore, a species’ evolutionary potential can be constrained by a lack of genetic variation (50). This underscores the need to preserve the genetic diversity of species in the face of accelerating environmental changes. Conversely, species with low genetic variation could fail to cope with accelerating anthropogenic environmental changes, which warrants more attention to their conservation status.
During the course of postglacial range expansion, parrotbills at the leading expansion edge are thought to have encountered different environmental challenges from those at the trailing edge, and high altitudes are conventionally considered to have harsher environments than low altitudes (11). We would thus expect high-altitude parrotbills to be under more stringent recent selection. To test the “harsh high altitudes” hypothesis, we examined the polymorphic patterns of the candidate SNPs at different altitudes and found similar numbers of fixed or nearly fixed candidate SNPs on both sides of the island (Fisher’s exact test P = 0.809 and 0.803 for the east and west pairs, Fig. 3). Assuming that fixation probability depends on the strength of directional selection in a given environment, this suggests that the strength of selection pressures at different altitudes was comparable. Evidence of strong selection pressure in the low-altitude populations also arose from the results of the Tajima’s D test that we used to detect recent selective sweep [<0.1 Ne generations (51, 52)], in 10-kb sliding regions across the genome. Recent selective sweeps that occurred since the last glacial period as indicated by the lowest 1% of Tajima’s D value in each local population were detected for 6 (Tajima’s D ranges from: −1.3569 to −1.6825) out of 24 candidate regions, with similar frequencies in population pairs on both sides of the CMR (none and four for the east high- and low-altitude populations, respectively; one and two for the west high- and low-altitude populations, respectively; Fisher’s exact test: P = 0.112 and P = 1.000 for the east and west population pairs, respectively).
The proportion of candidate SNPs that are fixed or nearly fixed (90% ≤ the major allele frequency < 100%) in the low- and high-altitude populations.
Our results suggest that during postglacial expansion, the population remaining at low altitudes (the trailing edge of the expansion) could also have experienced strong selection pressure. These low-altitude selection pressures may have been caused by increased temperatures, whereas those at the high-altitude leading edge of the expansion arose from the lower oxygen partial pressure. This finding calls for more research on adaptation at low altitudes or latitudes where populations at the trailing edge of postglacial expansion are located and are considered as long-term reservoirs of species’ genetic diversity and cradles of speciation (53).
To assemble a reference genome, we constructed two paired-end and mate-paired libraries for Illumina short-read sequencing from the DNA of a male Taiwan vinous-throated parrotbill (Sinosuthora webbiana bulomachus). The draft genome was assembled with ALLPATHS-LG (54) version-44099.
Repeated DNA sequences were masked with RepeatMasker (55) version open-4.0.5. We followed the procedure in Ellegren et al. (56) to remove potential contaminating DNA in the draft genome. Finally, 6,508 scaffolds were left for use in the draft parrotbill genome assembly. Then we used Satsuma (57) to align the parrotbill draft genome assembly to the well-annotated and assembled zebra finch (Taeniopygia guttata) genome. Furthermore, we used Augustus (58) for gene prediction and the ALDB database (59) for identifying lncRNAs.
Forty parrotbills (34 males and 6 females) were collected from two high-altitude and two low-altitude locations in Taiwan (SI Appendix, Table S3). Six females were sampled from the high-altitude local population in Hualien County; due to the low population density at high altitudes (60), other samples were all males. Climatic information (annual temperature and annual mean precipitation) for each locality was extracted from global climate data (Worldclim v1.4, ref. 61).
To infer whether SNP sites in Taiwan’s population are shared or private variants, we also sampled male parrotbills from four mainland populations with 10 individuals from each local population (SI Appendix, Table S4).
A whole-genome resequencing library was constructed for each individual. The average coverage of the parrotbill population from Taiwan and the mainland was 5.8× and 12.1×, respectively (SI Appendix, Tables S3 and S4).
We used the algorithm BWA-MEM to map the raw reads on the reference genome. Then we used Samtools (62) for variant calling and Vcf-tools (63) to generate a consensus sequence for each individual, and then we randomly assigned one allele from heterozygous genotypes to two putative haplotypes.
The demographic history of the vinous-throated parrotbill in Taiwan was estimated by a hidden Markov model (HMM) implanted in the pairwise sequentially Markovian coalescent (PSMC) method. Effective population sizes were inferred from autosomes of two individuals with average genome coverage of 24.5×.
We used G-PhoCS v1.3 (64) to estimate long-term Ne and the divergent time (τ) between the Taiwan and mainland parrotbill populations, based on genome sequences of four individuals (two from mainland and two from Taiwan) with high coverages (>20×).
The PopGenome package (65) in R with its sliding window method was used to calculate the following statistics in the Taiwan populations: interpopulation differentiation, FST and Dxy, intrapopulation diversity, π, and Tajima’s D with 10-kb nonoverlapping sliding windows. To evaluate linkage statistics, R2, within each 10-kb window, we used beagle version 4.1 (66) to phase all of the 40 individuals from Taiwan.
An unrooted neighboring-join tree based on the pairwise comparison of FST values between the four populations in Taiwan was constructed using MEGA7 (SI Appendix, Fig. S3) (67). Based on this result, four local populations were assigned to two high-/low-altitudinal pairs from the east and west side of the Central Mountain Range (CMR) in Taiwan.
To identify outlier regions that are likely associated with divergent selection between different altitudes we calculated the z-transformed FST value (zFST) of the two high-/low-altitudinal local population pairs of each 10-kb window. The 10-kb windows with zFST over 2.33 (the top 1%) were arbitrarily defined as outlier regions. Because there is a difference of about 580 m in altitude between the two high-altitude populations, which might result in some altitudinal effects, ΔFST was obtained by subtracting FST(low/low) from FST(high/low), and then z-transforming (zΔFST). Regions with higher ΔFST values are more likely to be under divergent than linked selection. Only regions with the top 1% values of both FST and ΔFST are referred to as candidate regions.
Genetic variants with the highest FST value within each candidate region were referred to as candidate SNPs. If these candidate SNPs were also found to be polymorphic in the mainland population, they were regarded as standing variants; otherwise, they were considered to be private variants.
Expression levels of a gene can be modulated by enhancers that are up to 1 Mb upstream or downstream of its start codon (68). We therefore defined genes within 1 Mb of the candidate region as candidate genes under divergent selection between different altitudes.
Protein sequences and GO annotations of zebra finch genes annotated by Ensembl were retrieved from BioMart. Orthologs of the vinous-throated parrotbill and zebra finch were identified using the InParanoid (v4.1, ref. 69) algorithm using default parameters). As a result, 8,635 vinous-throated parrotbill genes with at least one ortholog in the zebra finch genome were identified. The Ensembl-annotated GO terms of the zebra finch genes and the corresponding upstream GO terms retrieved using the R package “GO.db” (70) were assigned to the vinous-throated parrotbill genes. Fisher’s exact test (R function: fisher.test; null hypothesis no enrichment) was used to test the statistical significance of each GO term; the false discovering rate (FDR) was determined by R function: p.adjust.
We thank J. Q. Li, C. T. Wei, and L. W. Wang for providing parrotbill samples from China; Y. C. Zhang, T. C. Chu, and K. H. Chao for bioinformatics assistance; A. Watson for improving the manuscript; Y. Liu, F. Dong, W. N. Bai, and J. H. Hu for commenting on the early version of manuscript; and C. Y. Chen for illustration of parrotbill. This project was funded by a grant from Taiwan Ministry of Science and Technology (NSC 99 2321-B-003-004-MY3).
Author contributions: K.L. and S.-H.L. designed research; Y.-T.L., W.L., Y.-C.H., C.-T.Y., and L.D. performed research; W.L., Y.-C.H., C.-T.Y., and L.D. contributed new reagents/analytic tools; Y.-T.L., E.-L.P., Y.H., B.-Y.L., H.-F.C., B.-W.Z., C.-F.Y., C.-M.H., and M.-Y.Y. analyzed data; and Y.-T.L., C.K.L.Y., K.E.O., B.-Y.L., B.-W.Z., C.-M.H., H.-Y.H., M.-Y.Y., and S.-H.L. wrote the paper.
The authors declare no conflict of interest.
This article is a PNAS Direct Submission.
Data deposition: The sequence reported in this paper has been deposited in the National Center for Biotechnology Information Sequence Read Archive (accession nos. SRX5087906 to SRX5087990; BioProject ID PRJNA504683).
This article contains supporting information online at www.pnas.org/lookup/suppl/doi:10.1073/pnas.1813597116/-/DCSupplemental.
This open access article is distributed under Creative Commons Attribution-NonCommercial-NoDerivatives License 4.0 (CC BY-NC-ND).
