Difference in larval type explains patterns of nonsynonymous substitutions in two ancient paralogs of the histone H3 gene in sea stars

SUMMARY Paralogous genes frequently show differences in patterns and rates of substitution that are typically attributed to different selection regimes, mutation rates, or local recombination rates. Here, two anciently diverged paralogous copies of the histone H3 gene in sea stars, the tandem‐repetitive early‐stage gene and a newly isolated gene with lower copy number that was termed the “putative late‐stage histone H3 gene” were analyzed in 69 species with varying mode of larval development. The two genes showed differences in relative copy number, overall substitution rates, nucleotide composition, and codon usage, but similar patterns of relative nonsynonymous substitution rates, when analyzed by the d_N/d_S ratio. Sea stars with a nonpelagic and nonfeeding larval type (i.e., brooding lineages) were observed to have d_N/d_S ratios that were larger than for nonbrooders but equal between the two paralogs. This finding suggested that demographic differences between brooding and nonbrooding lineages were responsible for the elevated d_N/d_S ratios observed for brooders and refuted a suggestion from a previous analysis of the early‐stage gene that the excess nonsynonymous substitutions were due to either (1) gene expression differences at the larval stage between brooders and nonbrooders or (2) the highly repetitive structure of the early‐stage histone H3 gene.     

Evolutionary neutrality of mtDNA introgression: evidence from complete mitogenome analysis in roe deer

Introgressive hybridization offers a unique platform for studying the molecular basis of natural selection acting on mitogenomes. Most of the mtDNA protein‐coding genes are extremely conserved; however, some of the observed variations have potentially adaptive significance. Here, we evaluated whether the evolution of mtDNA in closely related roe deer species affected by widespread mtDNA introgression is neutral or adaptive. We characterized and compared 16 complete mitogenomes of European (Capreolus capreolus) and Siberian (C. pygargus) roe deer, including four of Siberian origin introgressed into European species. The average sequence divergence of species‐specific lineages was estimated at 2.8% and varied across gene classes. Only 21 of 315 fixed differences identified in protein‐coding genes represented nonsynonymous changes. Only three of them were determined to have arisen in the C. pygargus lineage since the time to the most recent common ancestor (TMRCA) of both Capreolus species, reflecting a decelerated evolutionary ratio. The almost four‐fold higher d_N/d_S ratio described for the European roe deer lineage is constrained by overall purifying selection, especially pronounced in the ND4 and ND5 genes. We suggest that the highly divergent C. capreolus lineage could have maintained a capability for genomic incorporation of the well‐preserved and almost ancestral type of mtDNA present in C. pygargus. Our analyses did not indicate any signs of positive selection for Siberian roe deer mtDNA, suggesting that the present widespread introgression is evolutionarily neutral.     

Invertebrate Species with Nonpelagic Larvae Have Elevated Levels of Nonsynonymous Substitutions and Reduced Nucleotide Diversities

Under a nearly neutral model in which most amino acid substitutions are slightly deleterious, variation in demography, population structure, and other ecological factors among closely related species can potentially modify the effective population size or the selective regime, leading to differences in the rate of nonsynonymous substitution. Ratios of nonsynonymous to synonymous substitutions (d_N/d_S) between species were analyzed in a sea star genus (Patiriella) and a molluscan genus (Littorina), each with diverse modes of reproduction, including multiple lineages with pelagic and nonpelagic larvae. In both genera, lineages with nonpelagic larvae had significantly higher d_N/d_S ratios than lineages with pelagic larvae. The hypothesis that the elevated d_N/d_S ratios in species with nonpelagic larvae was due to reduced effective population size was tested by comparing nucleotide diversities in three genera of gastropod mollusks (Littorina, Crepidula, and Hydrobia), each with several modes of reproduction. Overall, there was a significant (p < 0.05) reduction in nucleotide diversity in species with nonpelagic larvae compared to species with pelagic larvae.     

Sequence analysis of a polymorphic Mhc class II gene in Pacific salmon

Polymorphism of the nucleotide sequences encoding 149 amino acids of linked major histocompatibility complex (Mhc) class II 131 and 132 peptides, and of the intervening intron (548–773 base pairs), was examined within and among seven Pacific salmon (Oncorhynchus) species. Levels of nucleotide diversity were higher for theB1 sequence than forB2 or the intron in comparisons both within and between species. For the codons of the peptide binding region of the BI sequence, the level of nonsynonymous nucleotide substitution (d_N) exceeded the level of synonymous substitution (d_S) by a factor of ten for within-species comparisons, and by a factor of four for between-species comparisons. The excess of d_N indicates that balancing selection maintains diversity at this salmonidMhc class II locus, as is common forMhc loci in other vertebrates. Levels of nucleotide diversity for both the exon and intron sequences were greater among than within species, and there were numerous species-specific nucleotides present in both the coding and noncoding regions. Thus, neighbor-joining analysis of both the intron and exon regions provided phylogenies in which the sequences clustered strongly by species. There was little evidence of shared ancestral (trans-species) polymorphism in the exon phylogeny, and the intron phylogeny depicted standard relationships among the Pacific salmon species. The lack of shared allelicB1 lineages in these closely related species may result from severe bottlenecks that occurred during speciation or during the ice ages that glaciated the rim of the north Pacific Ocean approximately every 100 000 years in the Pleistocene.The nucleotide sequence data reported in this paper have been submitted to the GenBank nucleotide sequence database and have been assigned the accession numbers U34692-U34720     

Evolutionary rate variation among genes involved in galactomannan biosynthesis in Coffea canephora

The endosperm cell walls of mature coffee seeds accumulate large amounts of mannan storage polysaccharides, which serve as nutrient reserve for embryo and contribute to beverage quality. Our study investigated the evolutionary patterns of key galactomannan (GM) biosynthesis genes using d_N/d_S ratio, synteny, and phylogenetic analysis and detected heterogeneity in rate of evolution among gene copies. Selection ratio index revealed evidence of positive selection in the branch editing gene Coffea canephora alpha (α) galactosidase (Cc‐alpha Gal) at Cc11_g15950 copy (ω = 1.12), whereas strong purifying selection on deleterious mutations was observed in the Coffea canephora uridine diphosphate (UDP)‐glucose 4′‐epimerase (Cc‐UG4E) and Coffea canephora mannose‐1P guanylytransferase (Cc‐MGT) genes controlling the crucial nucleotide carbon sugar building blocks flux in the pathway. Relatively low sequence diversity and strong syntenic linkages were detected in all GM pathway genes except in Cc‐alpha Gal, which suggests a correlation between selection pressure and nucleotide diversity or synteny analysis. In addition, phylogenetic analysis revealed independent evolution or expansion of GM pathway genes in different plant species, with no obvious inferable clustering patterns according to either gene family or congruent with evolutionary plants lineages tested due to high dynamic nature and specific biochemical cell wall modification requirements. Altogether, our study shows a significant high rate of evolutionary variation among GM pathway genes in the diploid C. canephora and demonstrates the inherent variation in evolution of gene copies and their potential role in understanding selection rates in a homogenously connected metabolic pathway.     

The Genomic Selfing Syndrome Accompanies the Evolutionary Breakdown of Heterostyly

The evolutionary transition from outcrossing to selfing can have important genomic consequences. Decreased effective population size and the reduced efficacy of selection are predicted to play an important role in the molecular evolution of the genomes of selfing species. We investigated evidence for molecular signatures of the genomic selfing syndrome using 66 species of Primula including distylous (outcrossing) and derived homostylous (selfing) taxa. We complemented our comparative analysis with a microevolutionary study of P. chungensis, which is polymorphic for mating system and consists of both distylous and homostylous populations. We generated chloroplast and nuclear genomic data sets for distylous, homostylous, and distylous–homostylous species and identified patterns of nonsynonymous to synonymous divergence (d_N/d_S) and polymorphism (π_N/π_S) in species or lineages with contrasting mating systems. Our analysis of coding sequence divergence and polymorphism detected strongly reduced genetic diversity and heterozygosity, decreased efficacy of purifying selection, purging of large-effect deleterious mutations, and lower rates of adaptive evolution in samples from homostylous compared with distylous populations, consistent with theoretical expectations of the genomic selfing syndrome. Our results demonstrate that self-fertilization is a major driver of molecular evolutionary processes with genomic signatures of selfing evident in both old and relatively young homostylous populations.     

Evolutionary dynamics of Rh2 opsins in birds demonstrate an episode of accelerated evolution in the New World warblers (Setophaga)

Low rates of sequence evolution associated with purifying selection can be interrupted by episodic changes in selective regimes. Visual pigments are a unique system in which we can investigate the functional consequences of genetic changes, therefore connecting genotype to phenotype in the context of natural and sexual selection pressures. We study the RH2 and RH1 visual pigments (opsins) across 22 bird species belonging to two ecologically convergent clades, the New World warblers (Parulidae) and Old World warblers (Phylloscopidae) and evaluate rates of evolution in these clades along with data from 21 additional species. We demonstrate generally slow evolution of these opsins: both Rh1 and Rh2 are highly conserved across Old World and New World warblers. However, Rh2 underwent a burst of evolution within the New World genus Setophaga, where it accumulated substitutions at 6 amino acid sites across the species we studied. Evolutionary analyses revealed a significant increase in d_N/d_S in Setophaga, implying relatively strong selective pressures to overcome long‐standing purifying selection. We studied the effects of each substitution on spectral tuning and found they do not cause large spectral shifts. Thus, substitutions may reflect other aspects of opsin function, such as those affecting photosensitivity and/or dark–light adaptation. Although it is unclear what these alterations mean for colour perception, we suggest that rapid evolution is linked to sexual selection, given the exceptional plumage colour diversification in Setophaga.     

Molecular evolution of the clustered <Emphasis Type="Italic">MIC</Emphasis>-<Emphasis Type="Italic">3</Emphasis> multigene family of <Emphasis Type="Italic">Gossypium</Emphasis> species

The Gossypium MIC-3 (Meloidogyne Induced Cotton-3) gene family is of great interest for molecular evolutionary studies because of its uniqueness to Gossypium species, multi-gene content, clustered localization, and root-knot nematode resistance-associated features. Molecular evolution of the MIC-3 gene family was studied in 15 tetraploid and diploid Gossypium genotypes that collectively represent seven phylogenetically distinct genomes. Synonymous (d_S) and non-synonymous (d_N) nucleotide substitution rates suggest that the second of the two exons of the MIC-3 genes has been under strong positive selection pressure, while the first exon has been under strong purifying selection to preserve function. Based on nucleotide substitution rates, we conclude that MIC-3 genes are evolving by a birth-and-death process and that a ‘gene amplification’ mechanism has helped to retain all duplicate copies, which best fits with the “bait and switch” model of R-gene evolution. The data indicate MIC-3 gene duplication events occurred at various rates, once per 1 million years (MY) in the allotetraploids, once per ~2 MY in the A/F genome clade, and once per ~8 MY in the D-genome clade. Variations in the MIC-3 gene family seem to reflect evolutionary selection for increased functional stability, while also expanding the capacity to develop novel “switch” pockets for responding to diverse pests and pathogens. Such evolutionary roles are congruent with the hypothesis that members of this unique resistance gene family provide fitness advantages in Gossypium.     

Draft genomes and genomic divergence of two Lepidurus tadpole shrimp species (Crustacea,Branchiopoda, Notostraca)

Crustaceans of the order Notostraca (Branchiopoda) are distributed worldwide and are known for the remarkable morphological stasis between their extant and Permian fossil species. Moreover, these crustaceans show relevant ecological traits and a wide range of reproductive strategies. However, genomic studies on notostracans are fairly limited. Here, we present the genome sequences of two notostracan taxa, Lepidurus arcticus and Lepidurus apus lubbocki. Taking advantage of the small genome sizes (~0.11 pg) of these taxa, genomes were sequenced for one individual per species with one run on the Illumina HiSeq X platform. We finally assembled 73.2 Mbp (L. arcticus) and 90.3 Mbp (L. apus lubbocki) long genomes. Assemblies cover up to 84% of the estimated genome size, with a gene completeness >97% for both genomes. In total, 13%–16% of the assembled genomes consist of repeats, and based on read mapping, L. apus lubbocki shows a significantly lower transposable element content than L. arcticus. The analysis of 2,376 orthologous genes indicates an ~7% divergence between the two Lepidurus taxa, with a nucleotide substitution rate significantly lower than that of Daphnia taxa. K_a/K_s analysis suggests purifying selection in both branchiopod lineages, raising the question of whether the low substitution rate of Lepidurus is correlated with morphological conservation or is linked to specific biological traits. Our analysis demonstrates that, in these organisms, it is possible to obtain high‐quality draft genomes from single individuals with a relatively low sequencing effort. This result makes Lepidurus and Notostraca interesting models for genomic studies at taxonomic, ecological and evolutionary levels.     

Increased positive selection pressure within the complementarity determining regions of the T‐cell receptor β gene in New World monkeys

Because of the long‐term co‐evolution of TCR and MHC molecules, numerous nucleotide substitutions have accumulated within the domains of TCRβ genes. We previously found that nonsynonymous nucleotide substitutions occurred more frequently in complementarity determining region (CDR)β than in CDRα, even though only a limited number of common marmoset (Callithrix jacchus) and human T‐cell receptor β variable (TRBV) sequences were compared. This interesting finding raised the question of whether the increased selective pressure within CDRβ was species‐specific. In this study, we identified 21 TRBV region sequences from the common marmoset and performed comparative sequence analyses of the T‐cell receptor α variable (TRAV) and TRBV regions from human, chimpanzee, rhesus monkey, cotton‐top tamarin, Ma's night monkey, and common marmoset. The ratios of the number of nonsynonymous nucleotide substitutions per site (d_N) to the d_S values (d_N/d_S) were less than 1 within the framework regions (FRs) of TRAV and TRBV region sequences, suggesting that purifying selection is largely dominant within the FRs. In contrast, the d_N values were statistically significantly greater for CDRβ than for CDRα only in New World monkeys. Also, increased d_N/d_S ratios (d_N/d_S>1) were observed within CDRβ between humans and New World monkeys and, interestingly, between New World monkeys, which share a relatively recent common ancestor. Moreover, phylogenetic analysis by maximum likelihood analysis provided firm evidence to support that positive selection occurred within CDRβ along New World monkey lineages. These results suggest that increased positive selection pressure within CDRβ is common in New World monkeys rather than being species‐specific. This study provides an intriguing insight into the co‐evolution of TCR and MHC molecules within primates. Am. J. Primatol. 73:1082–1092, 2011. © 2011 Wiley‐Liss, Inc.     

Likelihood Analysis of the Chalcone Synthase Genes Suggests the Role of Positive Selection in Morning Glories (<Emphasis Type="Italic">Ipomoea</Emphasis>)

Chalcone synthase (CHS) is a key enzyme in the biosynthesis of flavonoides, which are important for the pigmentation of flowers and act as attractants to pollinators. Genes encoding CHS constitute a multigene family in which the copy number varies among plant species and functional divergence appears to have occurred repeatedly. In morning glories (Ipomoea), five functional CHS genes (A–E) have been described. Phylogenetic analysis of the Ipomoea CHS gene family revealed that CHS A, B, and C experienced accelerated rates of amino acid substitution relative to CHS D and E. To examine whether the CHS genes of the morning glories underwent adaptive evolution, maximum-likelihood models of codon substitution were used to analyze the functional sequences in the Ipomoea CHS gene family. These models used the nonsynonymous/synonymous rate ratio ( = d_N/d_S) as an indicator of selective pressure and allowed the ratio to vary among lineages or sites. Likelihood ratio test suggested significant variation in selection pressure among amino acid sites, with a small proportion of them detected to be under positive selection along the branches ancestral to CHS A, B, and C. Positive Darwinian selection appears to have promoted the divergence of subfamily ABC and subfamily DE and is at least partially responsible for a rate increase following gene duplication.     

Correlated evolution of nucleotide substitution rates and allelic variation in Mhc-DRB lineages of primates

Background  

The major histocompatibility complex (MHC) is a key model of genetic polymorphism. Selection pressure by pathogens or other microevolutionary forces may result in a high rate of non-synonymous substitutions at the codons specifying the contact residues of the antigen binding sites (ABS), and the maintenance of extreme MHC allelic variation at the population/species level. Therefore, selection forces favouring MHC variability for any reason should cause a correlated evolution between substitution rates and allelic polymorphism. To investigate this prediction, we characterised nucleotide substitution rates and allelic polymorphism (i.e. the number of alleles detected in relation to the number of animals screened) of several Mhc class II DRB lineages in 46 primate species, and tested for a correlation between them.     

Post‐glacial colonization routes coincide with a life‐history breakpoint along a latitudinal gradient

Although adaptive divergence along environmental gradients has repeatedly been demonstrated, the role of post‐glacial colonization routes in determining phenotypic variation along gradients has received little attention. Here, we used a hierarchical Q_ST–F_ST approach to separate the roles of adaptive and neutral processes in shaping phenotypic variation in moor frog (Rana arvalis) larval life histories along a 1,700 km latitudinal gradient across northern Europe. This species has colonized Scandinavia via two routes with a contact zone in northern Sweden. By using neutral SNP and common garden phenotypic data from 13 populations at two temperatures, we showed that most of the variation along the gradient occurred between the two colonizing lineages. We found little phenotypic divergence within the lineages; however, all phenotypic traits were strongly diverged between the southern and northern colonization routes, with higher growth and development rates and larger body size in the north. The Q_ST estimates between the colonization routes were four times higher than F_ST, indicating a prominent role for natural selection. Q_ST within the colonization routes did not generally differ from F_ST, but we found temperature‐dependent adaptive divergence close to the contact zone. These results indicate that lineage‐specific variation can account for much of the adaptive divergence along a latitudinal gradient.     

Complex dynamics underlie the evolution of imperfect wing pattern convergence in butterflies

Adaptive radiation is characterized by rapid diversification that is strongly associated with ecological specialization. However, understanding the evolutionary mechanisms fueling adaptive diversification requires a detailed knowledge of how natural selection acts at multiple life‐history stages. Butterflies within the genus Adelpha represent one of the largest and most diverse butterfly lineages in the Neotropics. Although Adelpha species feed on an extraordinary diversity of larval hosts, convergent evolution is widespread in this group, suggesting that selection for mimicry may contribute to adaptive divergence among species. To investigate this hypothesis, we conducted predation studies in Costa Rica using artificial butterfly facsimiles. Specifically, we predicted that nontoxic, palatable Adelpha species that do not feed on host plants in the family Rubiaceae would benefit from sharing a locally convergent wing pattern with the presumably toxic Rubiaceae‐feeding species via reduced predation. Contrary to expectations, we found that the presumed mimic was attacked significantly more than its locally convergent model at a frequency paralleling attack rates on both novel and palatable prey. Although these data reveal the first evidence for protection from avian predators by the supposed toxic, Rubiaceae‐feeding Adelpha species, we conclude that imprecise mimetic patterns have high costs for Batesian mimics in the tropics.     

A revised evolutionary history of hepatitis B virus (HBV)

Previous studies of the evolutionary history of hepatitis B virus (HBV) have been compromised by intergenotype recombination and complex patterns of nucleotide substitution, perhaps caused by differential selection pressures. We examined the phylogenetic distribution of recombination events among human HBV genotypes and found that genotypes A plus D, and genotypes B plus C, had distinct patterns of recombination suggesting differing epidemiological relationships among them. By analyzing the nonoverlapping regions of the viral genome we found strong bootstrap support for some intergenotypic groupings, with evidence of a division between human genotypes A–E from the viruses sampled from apes and human genotype F. However, the earliest events in the divergence of HBV remain uncertain. These uncertainties could not be explained by differential selection pressures, as the ratio of nonsynonymous-to-synonymous substitutions (d _N/d _S) did not vary extensively among lineages and there is no strong evidence for positive selection across the whole tree. Finally, we provide a new estimate of the mean substitution rate in HBV, 4.2 × 10⁻⁵, which suggests that divergence of HBV in humans and apes has occurred only in the last 6000 years.     

Whole Genome Characterization,Phylogenetic and Genome Signature Analysis of Human Pandemic H1N1 Virus in Thailand, 2009–2012

Background

Three waves of human pandemic influenza occurred in Thailand in 2009–2012. The genome signature features and evolution of pH1N1 need to be characterized to elucidate the aspects responsible for the multiple waves of pandemic.

Methodology/Findings

Forty whole genome sequences and 584 partial sequences of pH1N1 circulating in Thailand, divided into 1^st, 2^nd and 3^rd wave and post-pandemic were characterized and 77 genome signatures were analyzed. Phylogenetic trees of concatenated whole genome and HA gene sequences were constructed calculating substitution rate and d_N/d_S of each gene. Phylogenetic analysis showed a distinct pattern of pH1N1 circulation in Thailand, with the first two isolates from May, 2009 belonging to clade 5 while clades 5, 6 and 7 co-circulated during the first wave of pH1N1 pandemic in Thailand. Clade 8 predominated during the second wave and different proportions of the pH1N1 viruses circulating during the third wave and post pandemic period belonged to clades 8, 11.1 and 11.2. The mutation analysis of pH1N1 revealed many adaptive mutations which have become the signature of each clade and may be responsible for the multiple pandemic waves in Thailand, especially with regard to clades 11.1 and 11.2 as evidenced with V731I, G154D of PB1 gene, PA I330V, HA A214T S160G and S202T. The substitution rate of pH1N1 in Thailand ranged from 2.53×10⁻³±0.02 (M2 genes) to 5.27×10⁻³±0.03 per site per year (NA gene).

Conclusions

All results suggested that this virus is still adaptive, maybe to evade the host''s immune response and tends to remain in the human host although the d_N/d_S were under purifying selection in all 8 genes. Due to the gradual evolution of pH1N1 in Thailand, continuous monitoring is essential for evaluation and surveillance to be prepared for and able to control future influenza activities.     

The evolution of mitochondrial genomes in modern frogs (Neobatrachia): nonadaptive evolution of mitochondrial genome reorganization

Background

Although mitochondrial (mt) gene order is highly conserved among vertebrates, widespread gene rearrangements occur in anurans, especially in neobatrachians. Protein coding genes in the mitogenome experience adaptive or purifying selection, yet the role that selection plays on genomic reorganization remains unclear. We sequence the mitogenomes of three species of Glandirana and hot spots of gene rearrangements of 20 frog species to investigate the diversity of mitogenomic reorganization in the Neobatrachia. By combing these data with other mitogenomes in GenBank, we evaluate if selective pressures or functional constraints act on mitogenomic reorganization in the Neobatrachia. We also look for correlations between tRNA positions and codon usage.

Results

Gene organization in Glandirana was typical of neobatrachian mitogenomes except for the presence of pseudogene trnS (AGY). Surveyed ranids largely exhibited gene arrangements typical of neobatrachian mtDNA although some gene rearrangements occurred. The correlation between codon usage and tRNA positions in neobatrachians was weak, and did not increase after identifying recurrent rearrangements as revealed by basal neobatrachians. Codon usage and tRNA positions were not significantly correlated when considering tRNA gene duplications or losses. Change in number of tRNA gene copies, which was driven by genomic reorganization, did not influence codon usage bias. Nucleotide substitution rates and d_N/d_S ratios were higher in neobatrachian mitogenomes than in archaeobatrachians, but the rates of mitogenomic reorganization and mt nucleotide diversity were not significantly correlated.

Conclusions

No evidence suggests that adaptive selection drove the reorganization of neobatrachian mitogenomes. In contrast, protein-coding genes that function in metabolism showed evidence for purifying selection, and some functional constraints appear to act on the organization of rRNA and tRNA genes. As important nonadaptive forces, genetic drift and mutation pressure may drive the fixation and evolution of mitogenomic reorganizations.

Electronic supplementary material

The online version of this article (doi:10.1186/1471-2164-15-691) contains supplementary material, which is available to authorized users.     

Molecular evolution of <Emphasis Type="Italic">PKD2</Emphasis> gene family in mammals

PKD2 gene encodes a critical cation channel protein that plays important roles in various developmental processes and is usually evolutionarily conserved. In the present study, we analyzed the evolutionary patterns of PKD2 and its homologous genes (PKD2L1, PKD2L2) from nine mammalian species. In this study, we demonstrated the orthologs of PKD2 gene family evolved under a dominant purifying selection force. Our results in combination with the reported evidences from functional researches suggested the entire PKD2 gene family are conserved and perform essential biological roles during mammalian evolution. In rodents, PKD2 gene family members appeared to have evolved more rapidly than other mammalian lineages, probably resulting from relaxation of purifying selection. However, positive selection imposed on synonymous sites also potentially contributed to this case. For the paralogs, our results implied that PKD2L2 genes evolved under a weaker purifying selection constraint than PKD2 and PKD2L1 genes. Interestingly, some loop regions of transmembrane domain of PKD2L2 exhibited higher P _N/P _S ratios than expected, suggesting these regions are more functional divergent in organisms and worthy of special attention. Electronic supplementary material  The online version of this article (doi:) contains supplementary material, which is available to authorized users. Chun Ye, Huan Sun have contributed equally to this work.