extfieldsofstudy: Medicine Geography

text:
Inter- and intracontinental migrations and local differentiation have shaped the contemporary epidemiological landscape of canine parvovirus in South America

Abstract Canine parvovirus (CPV) is a fast-evolving single-stranded DNA virus that causes one of the most significant infectious diseases of dogs. Although the virus dispersed over long distances in the past, current populations are considered to be spatially confined and with only a few instances of migration between specific localities. It is unclear whether these dynamics occur in South America where global studies have not been performed. The aim of this study is to analyze the patterns of genetic variability in South American CPV populations and explore their evolutionary relationships with global strains. Genomic sequences of sixty-three strains from South America and Europe were generated and analyzed using a phylodynamic approach. All the obtained strains belong to the CPV-2a lineage and associate with global strains in four monophyletic groups or clades. European and South American strains from all the countries here analyzed are representative of a widely distributed clade (Eur-I) that emerged in Southern Europe during 1990–98 to later spread to South America in the early 2000s. The emergence and spread of the Eur-I clade were correlated with a significant rise in the CPV effective population size in Europe and South America. The Asia-I clade includes strains from Asia and Uruguay. This clade originated in Asia during the late 1980s and evolved locally before spreading to South America during 2009–10. The third clade (Eur-II) comprises strains from Italy, Brazil, and Ecuador. This clade appears in South America as a consequence of an early introduction from Italy to Ecuador in the middle 1980s and has experienced extensive local genetic differentiation. Some strains from Argentina, Uruguay, and Brazil constitute an exclusive South American clade (SA-I) that emerged in Argentina in the 1990s. These results indicate that the current epidemiological scenario is a consequence of inter- and intracontinental migrations of strains with different geographic and temporal origins that set the conditions for competition and local differentiation of CPV populations. The coexistence and interaction of highly divergent strains are the main responsible for the drastic epidemiological changes observed in South America in the last two decades. This highlights the threat of invasion from external sources and the importance of whole-genome resolution to robustly infer the origin and spread of new CPV variants. From a taxonomic standpoint, the findings herein show that the classification system that uses a single amino acid to identify variants (2a, 2b, and 2c) within the CPV-2a lineage does not reflect phylogenetic relationships and is not suitable to analyze CPV evolution. In this regard, the identification of clades or sublineages within circulating CPV strains is the first step towards a genetic and evolutionary classification of the virus.

Introduction
The parvoviruses (Parvoviridae family) are small (25 nm), nonenveloped viruses with icosahedral symmetry that infect a wide range of animal species. The subfamily Densovirinae infects invertebrates, while the subfamily Parvovirinae has diverse vertebrate hosts, including humans, pigs, rodents, cattle, and several domestic and wild carnivore species (Allison et al. 2013). Canine parvovirus (CPV) belongs to the species Carnivore protoparvovirus 1 (genus Protoparvovirus), a highly genetically related cluster of viruses that also comprises raccoon parvovirus, mink enteritis virus (MEV), and feline panleukopenia virus (FPV) (Cotmore et al. 2013).
CPV is a well-recognized cause of acute and often fatal hemorrhagic gastroenteritis in puppies and needs to be controlled by early vaccination (Decaro and Buonavoglia 2012). The virus has a linear single-stranded DNA genome (5.2 kb) with two open reading frames (ORFs) (Reed et al. 1988). The left ORF encodes nonstructural proteins 1 and 2 (NS1 and NS2), which are involved in DNA replication, capsid assembly, and viral cellular transport (Wang et al. 1998). The N-terminal region of NS1 and NS2 are identical in sequence, whereas the C-terminal region of NS2 is generated from differential splicing and is translated in a different reading frame. The right ORF encodes the viral capsid proteins 1 and 2 (VP1 and VP2), which are the main antigens that induce protective antibodies (Tsao et al. 1991;Agbandje et al.1995;Nelson et al. 2007). VP1 and VP2 are splice variants and are identical in sequence, except for a 143-amino-acid Nterminal region that is unique to VP1 and essential for viral cell infection (Tattersall and Ward 1976). In infectious virions, VP2 may be proteolytically cleaved near the amino terminus to form VP3 (Weichert et al. 1998). At both ends of the genome, there are non-translated regions with hairpin structures that prime viral DNA replication (Tattersall and Ward 1976).
CPV emerged in the late 1970s as a consequence of an interspecies jump from other carnivores to dogs, a process that involved a small number of amino acid changes in the VPs (Parrish 1999). The newly emerged CPV, referred as CPV type-2 (CPV-2), rapidly spread worldwide causing a pandemic in the new canine host (Decaro and Buonavoglia 2012). In 1979 appeared CPV-2a, a new genetic lineage that differs at five residues in VP2 and regained the ability to infect cats and other carnivores (Parrish et al. 1988;Truyen et al. 1996;Decaro et al. 2010Decaro et al. , 2011. This CPV-2a lineage replaced the original CPV-2 strains within two years and became the ancestor of all the variants that currently circulate worldwide. This early period of CPV emergence was characterized by global spreading and strain replacement (Hoelzer et al. 2008).
The CPV-2a lineage underwent further evolution, showing an intrinsically high substitution rate that resembled fastevolving RNA viruses (Shackelton et al. 2005). Changes in the residue 426 of VP2 have received considerable interest because they alter the antigenic properties of the capsid and have a geographic and temporal pattern. According to the amino acid in the 426 position, variants are denoted 2a (Asn), 2b (Asp), and 2c (Glu). The CPV-2a lineage was comprised of only 2a variants until the 2b variant was detected in the United States in 1984 (Parrish et al. 1991). The newly emerged 2b variant reached high prevalence in some countries and is currently circulating worldwide, but with different frequencies (Decaro and Buonavoglia 2012;Miranda and Thompson 2016). The 2c variant was identified around 2000 in Italy and is now detected mainly in Europe and America (Buonavoglia et al. 2001;Miranda and Thompson 2016).
Although the virus dispersed over long distances in the past, several current CPV populations have acquired particular nucleotide substitutions, which is evidence of local differentiation. This fact leads to the hypothesis that current populations are spatially confined and show only a few instances of migration between specific localities (Doki et al. 2006;Hoelzer et al. 2008;Clegg et al. 2011). However, it is unclear whether these dynamics persist in other scenarios due to the difficulty of comparing strains using only the 2a/2b/2c classification and partial VP2 sequences. To infer the emergence and spread of the genetic variants, it is critical to extend the analysis to other genomic regions and focus on virus populations from different places and periods.
In South America (Argentina, Brazil, Colombia, Ecuador, and Uruguay), CPV has shown singular characteristics that supports a significant temporal and spatial genetic variability (Pereira et al. 2007;Streck et al. 2009;Gallo Calderó n et al. 2012Pé rez et al. 2012;Aldaz et al. 2013;Duque-García et al. 2017). Events of local diversification were detected in Brazil during 1980-2000 by analyzing positively selected residues in VP1 (Pereira et al. 2007). In the early 2000s, the emergence and spread of 2c strains seemingly changed the epidemiological scenario of CPV in South America (Pé rez et al. 2007). More recently, a divergent 2a strain emerged in the Uruguayan 2c homogenous population and reached a frequency of 85% within 2 years (Maya et al. 2013).
Questions remain about the origin and spread of South American variants as most strains have been only partially sequenced; complete genomes available in South America are exclusively from Uruguay (Pé rez et al. 2014).
In this study, the complete coding genome of global CPV populations was analyzed to infer the origin and patterns of genetic variability in South America.

Genome dataset
To generate a comprehensive genome dataset, most full-length sequences available in the GenBank database were retrieved (Supplementary Table S2). The only genomes available from South America corresponded to forty-one Uruguayan strains that were previously reported by our research group (Pé rez et al. 2014).

DNA extraction, PCR, and sequencing
Viral DNA was extracted from dog feces using the QIAamp DNA Mini kit (Qiagen, Valencia, CA, USA). Full-length genome amplification was performed with a Long PCR Enzyme Mix (Thermo Scientific, Waltham, MA, USA), and PCR products were directly sequenced using internal primers according to Pé rez et al.

Sequence assembly and alignment
Sequence assembly and consensus sequences were obtained using the Lasergene Genomics Suite (DNASTAR). DNA sequences were aligned using the MAFFT algorithm (Katoh et al. 2002).

Recombination analysis
Identification of potential recombinant, parental sequences and localization of possible recombinant breakpoints were performed using the recombination detection program software, which implements seven distinct algorithms for characterization of recombinant sequences (Martin et al. 2015). Only potential recombination events detected by two or more algorithms, coupled with phylogenetic evidence of recombination, were considered significant using the highest acceptable P-value cutoff of 0.05.

Evolutionary rates and population dynamics
The evolutionary rates (measured as nucleotide substitutions per site per year, s/s/y), the ages of the most recent common ancestor (tMRCA) and internal nodes, and the geographic movement patterns of the CPV population were jointly estimated using the Bayesian Markov Chain Monte Carlo approach implemented in BEAST v1.7.5 package (Drummond and Rambaut 2007;Drummond et al. 2012). In order to minimize sampling bias, some countries were pooled together: China, India, and Japan were pooled as Asia, Albania, and Italy as Southern Europe and France and Germany as Western Europe. A single sample from Paraguay was pooled with Argentina. The Bayesian skyline plot method was used as the coalescent tree prior, under a GTRþIþG4 substitution rate model and a relaxed uncorrelated lognormal molecular clock model (Drummond et al. 2006). Migratory events were inferred from a Continuous Time Markov Chain model for the reconstruction of discrete ancestral characters. For comparison of reversible and nonreversible matrices of geographic movement, path, and stepping-stone samplings were used to determine Bayes factors (BFs). The Bayesian stochastic search variable selection approach was used to identify significant migratory events over non-zero rates. Markov Chain Monte Carlo was ran for 100 million generations and convergence of all parameters was checked by observing Effective Sample Size > 200 with TRACER v1.5 (Drummond and Rambaut 2007). The maximum clade credibility tree was obtained from the posterior distribution of trees using TreeAnnotator v1.8 (part of the BEAST package) after discarding 10% of steps as burn-in. For each node, posterior probability (PP) and posterior state probability (PSP) was collected. Phylogenetic uncertainty for each parameter was informed as 95% highest probability density (HPD) values. Maps, migration routes, trees and population effective number through time were plotted using in-house R scripts.

Informative amino acid residues and selective pressures
Identification of variable residues in the protein alignments was performed using DIVEIN web server (Deng et al. 2010). To evaluate selective pressures on individuals sites of the codon alignments, the synonymous (dS) and non-synonymous (dN) substitution rates per site were estimated using different methods available in DATAMONKEY web server: Single-Likelihood Ancestor Counting (SLAC), Fixed Effects Likelihood (FEL), Fast Unconstrained Bayesian AppRoximation (FUBAR), Internal Fixed Effects Likelihood (IFEL) and Mixed Effects Model of Evolution (MEME) (Pond and Frost 2005).

Strain sequencing
Sequences of sixty-three CPV genomes (4629 nt) were obtained from South American and European strains (Supplementary Table S1). The sixty-three newly sequenced genomes were added to the genomes retrieved from the GenBank (n ¼ 110) (Supplementary Table S2), to compile a dataset of 173 full-length sequences of Carnivore protoparvovirus 1. Sequences came from strains isolated mainly from dogs but also from raccoons (n ¼ 9), felines (n ¼ 11), and minks (n ¼ 7).

Recombination analysis
As recombination may lead to confounding the signal of vertical evolution when reconstructing phylogenetic trees, the first approach was to analyze the strains to detect recombination events. Thirteen recombinant strains were detected and characterized (Supplementary Tables S2 and S3).

Phylogenetic analysis using full-length genomes
An initial phylogenetic analysis was performed with the nonrecombinant dataset that included 160 genomes of Carnivore protoparvovirus 1. This analysis shows two distinct and highly supported genetic groups, which are referred as the FPV-and the CPV-like viruses (Hoelzer and Parrish 2010) (Fig. 1). The FPVlike group contains the viruses from felines and minks going back to the first isolates that were collected in the early 1960s. The CPV-like group splits into two branches: the original pandemic CPV-2 strains that were collected in the USA during 1978-79 and the CPV-2a lineage that includes the variants collected from 1980 onward (Fig. 1).

Phylodynamic analysis and identification of amino acid signatures
The dataset for the phylodynamic analysis (n ¼ 144) included only the viruses from the CPV-like group (CPV-2 and CPV-2a lineage). Most of the strains of the CPV-2a lineage come from domestic dogs but there are also some raccoon strains that appear grouped in the phylogeny. The tMRCA of the 2a lineage was estimated to be 36 years, calculated from the most recent analyzed strain (2014) with a 95% HPD ranging from 34 to 38 years. South American CPVs associate with global strains in four monophyletic groups or clades (Europe I, Asia I, Europe II, and South America I) with high PP ( Fig. 2A, Supplementary Figs S1 and S2).
The Europe I clade (Eur-I) includes European (Italy, France, and Albania) and South American (Argentina, Brazil, Ecuador, Paraguay, and Uruguay) strains ( Fig. 2A). This clade has Southern Europe as its most probable ancestral geographic location (PSP ¼ 0.57) ( Supplementary Fig. S3); and has strains from the United States and Western Europe (Germany) as sister branches ( Fig. 2A). The tMRCA of this clade was estimated to be 20 years with a 95% HPD interval of 16-24 years. The strains from this clade are relatively homogenous (p-distance: 2.7 Â 10 À3 and 1.6 Â 10 À3 for nucleotides and amino acids, respectively). This clade has a characteristic combination of residues that includes the 426E, which identifies strains originally denoted as 2c, but there is not a single amino acid in the NS and VP sequences that is unique (i.e. exclusively found) for the Eur-I clade ( Table 1). The strains of the Eur-I clade group with no geographic pattern, with the exception of Ecuadorian strains that cluster in a single branch and associate with Southern European strains ( Fig. 2A).
The Europe II clade (Eur-II) comprises strains from Europe (Italy) and South American strains from Ecuador and Northeast Brazil. The Eur-II clade has Southern Europe as the geographic origin (PSP ¼ 0.98) ( Supplementary Fig. S3) and a tMRCA of 30 years with a 95% HPD interval of 25-31 years ( Fig. 2A). This clade has high variability (5.6 Â 10 À3 and 4.4 Â 10 À3 for nucleotide and amino acids, respectively). The NS1 residue 544F is shared by all the Italian and South American strains (Table 1). Brazilian strains have characteristic residues in NS1 (369K, 583K, and 616N) and NS2 (92K, 125K); they also show a one nucleotide (G) deletion at position 2343 of the intronic region of VP. Ecuadorian strains have a 668E substitution in NS1, and a T440S change in VP2; some Ecuadorian strains have 139I and 426E residues in VP2 (Supplementary Table S4).

Effective population size through time
The global full-length sequences were also analyzed using the Bayesian skyline plot method to determine the population size fluctuations that occurred during the 1978-2014 period of the CPV epidemic. The Bayesian skyline plot reveals that CPV has undergone two relatively recent events of population size expansion. The first event occurred in the 1990s and the second one in the 2000s (Fig. 2B).

Significant dispersal pathways involving South American strains
The analysis herein dated the emergence of the CPV-2a lineage in the USA during 1978 (95% HPD: 1976-80). From then, the virus began expanding across continents, including Asia (BF ¼ 62) and Southern Europe (BF ¼ 56) (Fig. 3A and B). Southern Europe was a secondary hub from where the virus disseminated to Western Europe and South America (Fig. 3A). There are significant inter-and intracontinental migration events that involved South American strains (Fig. 3C). Highly significant intercontinental movements include the migration from Southern Europe to Ecuador (BF ¼ 139) and Argentina (BF ¼ 12), and the migration of the Asia-I clade from Asia to Uruguay (BF ¼ 91). The most significant regional movements (i.e. across countries borders) were from Argentina to Uruguay (BF ¼ 91) and Brazil (BF ¼ 49) (Fig. 3B and C).

Discussion
The original pandemic CPV (type 2) was described from dog samples collected simultaneously in 1978 from the USA, Europe and Australia (Appel et al. 1979). By 1980, the disease had spread worldwide because of the transportation of dogs or contaminated fomites (Hoelzer and Parrish 2010). CPV-2 was replaced 2 years later by the CPV-2a lineage (Parrish et al. 1988). The CPV-2a lineage has been incredibly successful and originated all the strains that are actually circulating in the field (Hoelzer and Parrish 2010). Data here obtained indicate that the CPV-2a lineage emerged in the USA (PSP ¼ 1) (Supplementary Fig. S3) during the late 1970 s and spread to Eurasia; this constitutes the first report of the geographic origin and the dispersal routes that this lineage underwent during its intercontinental spreading (Fig. 3A). Europe served as a second dispersion hub from where certain CPV strains arrived in South America (Figs 2A and 3C). Early records of CPV in South America dated back to 1980 in Brazil and Chile (Abalos et al. 1982;Pereira et al. 2007). Brazilian CPV strains belonged to the 2a lineage and were extremely similar to strains collected in other continents during the same period, representing a clear example of the global dispersion occurred during the early period of the CPV-2a lineage emergence (Pereira et al. 2007;Hoelzer et al. 2008). The present pattern of genetic variability in the continent is remarkably different; the heterogeneity in the viral population and the phylogenetic relationship with foreign strains suggests that long migrations and local differentiation have been driving forces for CPV evolution in South America in the recent past. Viral migrations seem to be acquiring a more predominant role also in other continents, as evidenced by the increase in the number of reports of strains that share high homology with variants from distant geographic regions (Zhao et al. 2015;Zhou et al. 2017).

Eur-I clade: origin and spread in South America
The findings here reported indicate that some European and South American strains from all the countries here analyzed are representative of a relevant clade (Eur-I) within the 2a lineage. The Eur-I clade was circulating in Europe during 1990-98 to later spread to South American territories. Southern Europe appears as the ancestral location of this clade. The first representative strains of this clade were reported in Italy in 2001 (Buonavoglia et al. 2001) and were denoted as '2c strains' based on a single amino acid position (426E). Notable, the Eur-I clade is not related with Italian strains collected during the late 1990s (1997-2000) ( Fig. 2A), suggesting that Eur-I clade would have emerged in another European country and later introduced in Italy around 2001. Strains with the 426E residue were circulating in Germany in 1996 and in Poland in 1995 but genome sequences are not available to confirm their precise origin (Decaro et al. 2007;Majer-Dziedzic et al. 2011). Soon after emerging, strains of the Eur-I clade successfully spread in several European dog populations in few years, replacing previous variants (Martella et al. 2005;Decaro et al. 2007Decaro et al. , 2009Decaro and Buonavoglia 2012). Interestingly, there was a significant rise in the CPV effective population size just before the emergence of this clade in Europe (Fig. 2B). It is possible that the ancestor of the Eur-I clade had a better fitness than previously circulating variants and infect dogs more effectively, leading to virus immune escape, vaccine failure and the subsequent population increase (Decaro et al. 2008, Decaro andBuonavoglia 2012). To our knowledge, there has been no previous phylodynamic evidence of a rapid expansion growth in the global CPV population associated with the emergence of a CPV strain with particular amino acid combinations (Fig. 2B, Table 1). Similar expansion events have been observed in other viruses with worldwide distribution, such as HIV and influenza, associated with the emergence of genetic variants (Ekong et al. 2017;Li et al. 2017).
The Eur-I strains invaded and spread in all South American countries analyzed so far (Argentina, Brazil, Ecuador, Paraguay, and Uruguay) (Pé rez et al. 2007(Pé rez et al. , 2012Streck et al. 2009;Calderó n et al. 2011;Pinto et al. 2012;Aldaz et al. 2013). Some regions of Uruguay, Argentina and Brazil have seen Eur-I strains becoming the only variant circulating in some period during the last decade. This spreading and replacement behavior was even more prominent than in Europe, where this clade has been circulating in some countries for almost 20 years but remains scarce or absent in others (Decaro and Buonavoglia 2012). The introduction in South America of the Eur-I clade was soon after its emergence in Europe and occurred independently in Argentina and Ecuador (Fig. 3C). From Argentina, the strains migrated to Uruguay and Brazil, revealing their successful spread and invasion capability. The Argentinean origin is coincident with epidemiological data, as strains from this clade were reported in Argentina during 2003(Calderó n et al. 2009), before the first description in Uruguay (2006) andBrazil (2008) (Pé rez et al. 2007;Streck et al. 2009). The second migration wave to Ecuador of the Eur-I clade is evidenced by an Ecuadorian subclade that is closely related with some strains of the ancestral South European CPV population.
The invasion and replacement dynamics of the Eur-I clade in South America are unique and had been only observed for the CPV-2a lineage that primary replaced CPV-2 and, in a lesser extent, for the 2b antigenic type that reached high frequencies in the USA during the 1980s (Parrish et al. 1991). Notable, the South American spread of the Eur-I clade correlates with a second CPV population expansion during the early 2000s (Fig. 2B). The occurrence of the population expansion inferred from genomic data is supported by an apparent increase in the number of cases in South America (Gallo Calderó n et al. 2012).
The Glu 426 residue, which was used to describe the variant 2c, is a retained ancestral character for the Eur-I clade but also occurs in unrelated strains of the Eur-II clade from Ecuador (Table 1) and therefore cannot be used as a diagnostic character for the group. The group identity is based on non-coding changes and particular residue combinations (molecular signature) in the NS and VP viral proteins (Table 1).

The introduction of the Asia-I clade in South America
Some Uruguayan strains are representative of a large Asiatic clade (Asia-I) within the CPV-2a lineage (Fig. 2A). This clade originated in Asia during the late 1980s and evolved in this continent by the acquisition of particular nucleotide and amino acid changes. The amino acid signature of the Asia-I clade is very characteristic and includes unique residues in the NS and VP (Table 1).
After emerging and evolving in Asia, the Asia-I clade was introduced to South America around 2009 (Fig. 3C). This scenario is completely concordant with previous studies; Asiatic strains from this clade were already reported in 2004 (Ying et al. 2009) but in Uruguay, they were first detected in 2010 (Pé rez et al. 2012). Interestingly, the Uruguayan strains reached a frequency of 85% in 2 years, replacing the previous variant of the Eur-I clade that was prevalent in the territory (Pé rez et al. 2012). Although the spreading was fast, this strain has not been described in the neighboring country of Argentina (Gallo Calderó n et al. 2015). This could be indicating that migration and replacement events are restricted to certain strains with particular characteristics and that not all strains have the same fitness in different regions. CPV dynamics seem to follow the changing connection between the countries as the introduction of the Asiatic clade in South America was associated with an increase in the cultural and economic connections between China and Uruguay.
The Uruguayan strains of the Asia-I clade are genetically homogenous, most likely because its presence in South America is the consequence of a single and recent introduction event. During the co-circulation of the Eur-I and Asia-I clades in South America, one recombinant strain was detected in Uruguay (Pé rez et al. 2014). This strain did not increase its frequency in the population to acquire a more significant role in CPV evolution. Similar sporadic events of recombination have been identified within CPV strains or between CPV and FPV or MEV (Mochizuki et al. 2008;Ohshima and Mochizuki 2009;Wang et al. 2012) (Supplementary Table S3). Remarkably, nine of the recombinant strains here detected are closely related and should be considered as circulating recombinant forms (Supplementary Table S3). These strains were collected in China from 2009 to 2014 and originated from a single recombinant event between an Asiatic strain of the Asia-I clade and another undetermined strain of the CPV-2a lineage (Wang et al. 2016). This circulating recombinant form is a clear evidence that recombination could strongly influence CPV evolution by creating new lineages.

Early introduction of European strains to South America
Some strains from Ecuador, Italy, and Brazil form a third clade (Eur-II) within the CPV-2a lineage. This clade appears in South America as a consequence of an early introduction from Europe to Ecuador in the middle 1980s ( Fig. 3B and C). The strains of the Eur-II clade have significant levels of variability (5.6 Â 10 À3 for nucleotide and 4.4 Â 10 À3 for amino acids) and signs of local differentiation (Table 1). The existence of relevant amino acid changes indicates that some viruses were undergoing diversification that may be related to selection pressures, a fact that is supported by the presence of positively selected residues in the molecular signature (Table 1) (Truyen et al. 1995;Decaro et al. 2009;Yoon et al. 2009). One remarkable change is the occurrence of the 426E residue in some Ecuadorian strains of the Eur-II clade. This residue emerged de novo in Ecuador and does not have the same ancestral origin that the 426E change of the Eur-I clade ( Fig. 2A, Table 1). Notably, Brazilian strains also show a conserved deletion in the VP1 intron that occurs in most isolates after 1990 but was not present in strains from other regions of the world (Pereira et al. 2007).

An exclusive South American genetic clade
Some South American strains from Argentina, Brazil and Uruguay constitute a clade (SA-I) that has evolved exclusively in South America and harbors a particular amino acid signature ( Table 1). The strains of the SA-I clade can be considered as excellent representatives of the evolution of CPV strains in the continent. The strains of the SA-I clade acquired changes by local differentiation and may have been circulating before the Eur-I and Asia-I introductions to South America. A similar situation occur in other countries that show spatial clustering and a limited international spread in different continents (Doki et al. 2006;Hoelzer et al. 2008;Clegg et al. 2011).

CPV classification and phylodynamics
The classification system that uses a single amino acid to identify variants (2a, 2b, and 2c) within the CPV-2a lineage does not reflect phylogenetic relationships of the strains and it is not suitable to analyze CPV evolution. In this regard, the identification of clades or sublineages within circulating CPV strains is the first step towards a genetic and evolutionary classification of the virus. The new classification system should take into account the phylogenetic clustering using the coding genome, the presence of unique amino acid signatures, and the ancestral geographic location of the genetic clade. The clade nomenclature can be based on the continent or territory of origin to refer to them in a consistent and standardized way. The wellcharacterized complete genome sequences of CPV here analyzed will serve as the reference or representative sequences for future studies that deal with newly identified CPV strains.

Conclusion
CPV evolution in South America is the consequence of the coexistence of different strains from multiple origins that sets the condition for competition and local genetic differentiation. Evolutionary dynamics involve locally differentiated populations that spread within or between continents and began to interact with CPV strains of the receptive population. One potential advantage of this type of evolution is the coexistence of highly divergent strains in the same territory that can compete and even recombine to allow the virus to evolve into more adapted forms. This could eventually produce impactful evolutionary events that can drastically change the epidemiological scenario in a short period. From a taxonomic standpoint, the results here obtained support that the amino acid changes that have been traditionally used for CPV typing (i.e. 426 for 2a/2b/2c strains) occurred independently in the populations and are therefore not suitable markers to study CPV evolution. Viral epidemiology requires a new classification system that takes into account phylodynamic data and the inclusion of recombinant circulating forms.

Supplementary data
Supplementary data are available at Virus Evolution online.
Conflict of interest: None declared.


Study of Different Qualitative Traits in Dolichos Bean (Dolichos lablab L., Var. typicus Prain) Germplasm

1 College of Horticulture, Sri Konda Laxman Telangana State Horticultural University, Rajendranagar, Hyderabad-500030, Telangana, India 2 Department of Horticulture, College of Agriculture, Professor Jayashankar Telangana State Agricultural University, Rajendranagar, Hyderabad-500030, Telangana, India 3 Department of Genetics and Plant Breeding, College of Horticulture, SKLTSHU, Rajendranagar, Hyderabad-500030, Telangana, India 4 Economic Botany, NBPGR Regional Station, Rajendranagar, Hyderabad-500030, Telangana, India

Introduction
Dolichos bean (Dolichos lablab L. var. typicus Prain) (2n=2x=22) belongs to family Leguminaceae sub-family Faboideae, tribe Phaseoleae and sub-tribe Phaseolinae which is commonly known as Indian bean, Hyacinth bean, sem, Egyptian kidney bean, bonavist bean, avarai and avari chikkudu is one of the most popular traditional vegetables extensively cultivated in India. The pods are naturally rich in carbohydrates, proteins, fat and fibers, as well as minerals which include Ca, P and Fe (Naeem et al., 2009).
The protein content of pods and seeds ranges from 10-19% and 15-25% respectively. Among the legumes, dolichos bean  March, 2017. Quality traits revealed that there is a considerable variability in dolichos bean germplasm for most of the traits like stem colour (green, dark green and purple), leaf vein colour (light green and green, purple), leaf density (sparse and intermediate), flower colour (purple, white and dark purple), growth habit (pole and bush), pod shape (straight, curved and intermediate), pod colour dark purple in IC-446584, pod beak (long, short and medium), pod suture (green, dark green and purple), pod curvature (straight and curved), pod surface (smooth and wrinkled), seed colour (black, brown and cream) and seed shape (flat, round and oblong). Whereas there was no variability observed in traits like leaf vein colour, pod suture and pod beak colour. These characters can be utilized as morphological markers in selections and breeding programmes.
constituents an important source of therapeutic agents used in the modern as well as traditional system of medicine (Morris, 2009). In fact it is considered as a multipurpose crop since it is used for forage, soil improvement, soil protection and weed control (Mass, 2005).

Materials and Methods
Thirty five dolichos bean germplasm lines along with the three checks were grown in Randomized Block Design with two replications. The data on qualitative characters were recorded as per minimal descriptors of NBPGR (Srivastava et al., 2001). Qualitative data on 13 traits were recorded on ten randomly selected plants in each genotype and the details of trait, classification and stage of scoring are presented in Table 1.

Results and Discussion
Thirty five genotypes of dolichos bean germplasm under present investigation were characterized based on 13 qualitative traits (Table 2) and (Plate 1.1 and 1.2). The present observations are similar with the findings of Chattopadhyay and Dutta (2010) and Chaitanya et al., (2014). Stem colour revealed that one genotype exhibited dark green stem colour, four genotypes showed purple stem colour and 30 genotypes showed green stem colour. Leaf vein colour of 28 genotypes exhibited light green in colour, five genotypes showed purple colour and 2 genotypes are green in colour. Leaf density of 30 genotypes was found to be intermediate and the remaining five genotypes showed sparse leaf density. Plant growth habit revealed that 33 genotypes were pole type plant growth habit and remaining two genotypes were found to be bush type plant growth habit.
Flower colour revealed that 25 genotypes exhibited white flower colour, nine genotypes showed purple flower colour and one genotype showed dark purple flower colour. Pod colour of 25 genotypes was found to be green in colour, three genotypes showed dark green in colour, five genotypes showed light green in colour, one genotype showed purple colour and remaining one genotype showed dark purple colour. Pod shape of six genotypes was straight, two genotypes were curved and remaining 27 genotypes were intermediate with respect to it.
Pod surface of 28 genotypes had smooth pod surface and remaining seven genotypes had wrinkled pod surfaces. Pod curvature of thirty genotypes showed curved pod curvature and remaining five genotypes showed straight pod curvature. Pod beak revealed that 19 genotypes were with short pod beak, 13 genotypes had medium pod beaks and remaining three genotypes with long pod beaks. Pod suture colour of thirty genotypes showed green colour pod suture and remaining five genotypes showed cream colour pod suture. Seed colour of sixteen genotypes seeds were found to be brown in colour, eight genotypes showed black colour seeds and remaining 11 genotypes seeds were cream in colour. Seed shape of 25 genotypes was found oblong, five genotypes were round in shape and remaining five genotypes were flat in shape.


Impact of G-CSF on expressions of Egr-1 and VEGF in acute ischemic cerebral injury

The aim of the present study was to investigate the protective effect of granulocyte colony-stimulating factor (G-CSF) on acute ischemic cerebral injury, and its mechanism through the impact of G-CSF on early growth response-1 (Egr-1) and vascular endothelial growth factor (VEGF) expressions. Male Sprague-Dawley (SD) rats were divided them into three groups, i.e., the sham, model and G-CSF groups to measure the effect of G-CSF on the volume of cerebral infarction and level of lactate dehydrogenase (LDH) in rats. Hematoxylin and eosin (H&E) staining method was performed for histopathological examination. Reverse transcription-polymerase chain reaction (RT-PCR) and western blot analysis were used to detect the mRNA and protein expressions of Egr-1 and VEGF in different groups. Furthermore, Statistical Product and Service Solutions (SPSS) 17.0 software was applied to detect the differences in the expression of Egr-1 and VEGF between the two groups. Compared with the sham group, we found that the volume of cerebral infarction and LDH content in the model group were significantly elevated. By contrast, in the model group, those indicators in the G-CSF group were obviously decreased. H&E staining results also showed that G-CSF could decrease the necrotic area in cerebral infarction and the incidence of inflammation, and sustain the integrity of the molecular structure. Immunofluorescence staining results revealed that the protein expressions of Egr-1 and VEGF in the model group were all significantly decreased, while those in the G-CSF group were remarkably elevated. RT-PCR and western blot analysis revealed that the mRNA and protein expressions of Egr-1 and VEGF in the model group were decreased obviously, but those in the G-CSF group were elevated significantly, and the differences between the two groups showed statistical significance (P<0.05). G-CSF manifests a significant protective effect on the acute ischemic cerebral injury, which may be realized through its effect on the expressions of Egr-1 and VEGF.

Introduction
In recent years, acute ischemic cerebral injury has become a kind of disease with an extremely high incidence rate and poor prognosis, and for some severe conditions, it can even be life-threatening (1). Generally, it is originated from the cerebral thrombosis and embolism, and, after onset of ischemia, timely recovery of blood supply is critical to the reduction of apoptosis in neurons and repair of neurological functions. Ischemic angiogenesis is closely related to the progression and outcome of cerebrovascular diseases. However, current treatment method for acute ischemic cerebral injury requires further research and development (2)(3)(4).
Granulocyte colony-stimulating factor (G-CSF), as a member of hematopoietic growth factor family, can promote the proliferation, differentiation and survival of hematopoietic stem cells with the obvious protective effect on neurons in peripheral and central nerve system (5). In the present study, we aimed to establish the models of acute ischemic cerebral injury in rats to investigate the changes in early growth response-1 (Egr-1) and vascular endothelial growth factor (VEGF) in cerebral tissues after intervention with G-CSF, thereby exploring the potential mechanism of G-CSF in protecting nerves, and providing new pathways for clinical treatment of acute ischemic cerebral injury.

Materials
Experiment animals and grouping. After Sprague-Dawley (SD) rats (180-220 g) were acclimatized to the environment, they were randomly divided into the sham group (n=10), the model group (n=10) and the G-CSF group (100 µg/kg, n=10). G-CSF was given through neck subcutaneous injection once per day for 7 consecutive days. Subsequently, access of food to rats was forbidden for 12 h but they had free access to water. The rats were kept in cage with controlled temperature and light cycles (24̊C and 12/12 light cycles). Longa's method was applied to establish the middle cerebral artery occlusion models in rats. The humidity was 60±10%. The study was approved by the Ethics Committee of Zhongshan Torch Development Zone Hospital (Zhongshan, China).

Experiment methods Measurement of the volume of cerebral infarction and relevant indicators.
After rats in each group were anesthetized, they were immediately sacrificed and brain tissue was removed. Samples were then stained for 30 min at 37˚C in the dark in 1% triphenyltetrazolium chloride (TTC), and then fixed in 10% formaldehyde. The infarction area in each layer was then calculated, and LDH content in brain tissues in the three groups was assayed in accordance with the LDH kit.
Histopathological examination. Brain tissues in each group were embedded in paraffin to prepare the paraffin samples which were later sliced into sections of 5 µm. After hematoxylin and eosin (H&E) staining, the sections were placed under a microscope (x100) (Roche Diagnostics, Basel, Switzerland) for histopathological observation.
Immunofluorescence analysis. Paraffin sections of brain tissues were collected from each group, and, after dewaxing with xylene and dehydration using alcohol of gradient concentration, the antigen retrieval was carried out. Sections were then rinsed with 0.01 M phosphate-buffered saline (PBS, pH 7.4) 3 times (5 min/time), and blocked in a wet box containing 10% bull serum albumin (BSA) for 30 min. Subsequently, on sections, fluorescence-labeled antibodies that were diluted appropriately (1:70) were placed in the wet box and the sections were then incubated at 4˚C overnight. After sections were rinsed with PBS (pH 7.4) 3 times, fluorescent secondary antibodies (1:100) were added onto the sections in the dark, and placed in the wet box for incubation for 2 h at 37˚C. The sections were then observed and photographed under the fluorescent microscopes (Olympus Corporation, Tokyo, Japan).
RT-PCR analysis. Brain tissues that were collected from each group were rapidly transferred to the Eppendorf (EP) tubes containing RNAiso Plus extraction solution and placed at room temperature for 5 min for sufficient lysis. After centrifugation at 12,000 x g and 4˚C for 5 min, the supernatant was well mixed with 0.2 ml chloroform and placed at room temperature for 5 min. The solution was centrifuged again at 12,000 x g and 4˚C for 15 min, and, isopropanol of the same volume was added and well mixed to the supernatant. After being placed at room temperature for 10 min, it was centrifuged at 12,000 x g and 4˚C for 10 min, and sediment was preserved while the supernatant was discarded. In the sediment, 1 ml 75% alcohol was added and mixed well followed by centrifugation at 12,000 x g and 4˚C for 5 min, which was repeated once after the supernatant was then removed. After RNA sediment was fully removed through rinsing, the liquid was replaced with RNase-free water. Some of the total RNA solution was extracted and diluted with RNase-free water to a concentration of 1 µg/µl. In accordance with the requirement of the PrimeScript ® RT reagent kit with gDNA Eraser, reaction solution for reverse transcription was prepared, and cDNA was obtained through reverse transcription with the RNA samples and preserved at -20˚C. Thereafter, according to the instruction of SYBR ® Premix Ex Taq™ II (Tli RNaseH Plus), we detected the level of mRNA expression. The relevant RNA primer sequences used in this procedure are shown in Table I.
Western blotting assay. Brain tissues collected from all the groups were rinsed with icy normal saline followed by measurement of protein concentration with BCA kit and preservation at -80˚C for later use. According to the instructions of total protein extraction kit, the IP lysis buffer supplemented with phenylmethanesulfonyl fluoride (PMSF) and protease inhibitor was added in the tissues which were later ground on ice. Tissue samples were then centrifuged at 4˚C and 12,000 x g for 10 min of homogenization, and the supernatant was taken for centrifugation at 4˚C and 12,000 x g for 20 min. With the supernatant, protein samples containing the same content of total protein, after protein quantification according to the relevant kit instruction, were loaded in the sampling holes for electrophoresis under a constant voltage of 220 V. Electrophoresis was stopped until the bromophenol blue reached the bottom of gel. Based on the molecular weight of targeting proteins, gel was trimmed and placed in the transfer buffer. Polyvinylidene fluoride (PVDF) membrane was also trimmed in accordance with the size of the gel, and then soaked in methanol for 10 sec. The PVDF membrane and filter paper were placed in the transfer buffer, and later, those materials were all placed in the transfer machine in the following sequence: positive electrode, filter paper, PVDF membrane, gel, filter paper and negative electrode. Membrane transfer was then carried out under a constant voltage of 110 V for appro-

Gene name
Primer sequence priate time until the protein on the gel was transferred on the PVDF membrane. Thereafter, PVDF membrane with protein was blocked in 5% skimmed milk on a shaker for 3 h, and the rabbit anti-rat Egr-1, and GAPDH monoclonal antibodies (1:1,000; cat. nos. 4154 and 2118), rabbit anti-rat VEGF monoclonal antibody (1:1,000; cat. no. 19003-1-AP) was added on the membrane for incubation at 4˚C overnight. The next day, membrane was sufficiently washed by Tween-20 + Tris-buffer saline (TTBS) 3 times (10 min/time), and then incubated with goat anti-rabbit secondary polyclonal antibody (1:2,000; cat. no. 7074) for 1 h at room temperature. Enhanced chemiluminescence (ECL) reagent was added for color development and photographing after the membrane was washed with TTBS 3 times (10 min/time).

Statistical analysis. Statistical Product and Service Solutions
(SPSS) 17.0 software (SPSS, Inc., Chicago, IL, USA) was applied in statistical analysis, and experiment data were presented as mean ± standard deviation (SD). Data analysis was carried out through analysis of variance and the post hoc test was LSD test for multiple comparisons or t-test for intergroup comparisons. P<0.05 suggested that the difference had statistical significance.

Effect of G-CSF on the volume of cerebral infarction and
LDH content in brain tissues in rats. As shown in Fig. 1 in comparison with the sham group, the volume of cerebral infarction in rats of the model group was significantly enlarged.
By contrast, compared with the model group, a significant decrease was identified in the volume of cerebral infarction in rats of the G-CSF group. In addition, G-CSF could reduce the level of LDH in brain tissues significantly.
Effect of G-CSF on the histopathological examination of rats with acute ischemic cerebral injury. The H&E staining results of histopathological examination of brain tissues are shown in Fig. 2. In the sham group, the brain structure of rats appeared to be integral, while rats in the model group manifested massive necrotic cells, and destructed molecular structures. Compared with the model group, we found that the necrotic area was reduced in the G-CSF group, suggesting the cerebral injury was ameliorated.
Immunofluorescent staining results of Egr-1 and VEGF. The immunofluorescence analysis revealed that compared with the sham group, the protein expressions of Egr-1 and VEGF in brain tissues of rats in the model group were reduced significantly, while those in the G-CSF group were increased evidently (P<0.01) (Fig. 3). Fig. 4, the RT-PCR results revealed that compared with the sham group, the mRNA expressions of Egr-1 and VEGF in brain tissues of rats in the model group were reduced significantly. By contrast, in comparison with the model group, the mRNA expression of Egr-1 and VEGF was increased significantly in brain tissues of rats in the G-CSF group.  Western blot analysis results of the effect of G-CSF on the protein expression of Egr-1 and VEGF in MCAO rats. In Fig. 5, the western blot analysis revealed that in comparison with the sham group, protein expression of Egr-1 and VEGF in the brain tissues of rats in the model group was significantly reduced. Compared with the model group, the protein expression of Egr-1 and VEGF in brain tissues of rats in the G-CSF group was significantly elevated. This result confirmed that G-CSF can ameliorate the acute ischemic cerebral injury by regulating the expressions of Egr-1 and VEGF.

Discussion
Acute ischemic cerebral injury has gradually become a global killer with characteristics including high incidence rate, high mortality rate and high morbidity rate (6,7). Cerebral ischemia can damage the brain cells, and once the blood perfusion is recovered, ischemic injury is further aggravated, known as cerebral ischemia-reperfusion injury (8)(9)(10). Pathogenesis leading to acute ischemic cerebral injury remains elusive, which makes the search for effective, reliable and feasible treatment methods for this disease more urgent. In recent years, G-CSF has gained more and more attention of researchers for its special effect (11). G-CSF, as a kind of cytokine in hematopoietic stem cells, can promote the survival, proliferation and maturity of neutrophilic cell lines. G-CSF exists mostly in the bone marrow stem cell of patients (12)(13)(14). Previous findings have shown that G-CSF has a promising protective effect on nerves, and it can promote angiogenesis and exert an anti-inflammatory effect by regulating varying molecular mechanisms, such as mobilizing the generation of hematopoietic stem cells, inhibiting their apoptosis, facilitating the differentiation of neurons and antagonizing the inflammatory responses (15)(16)(17)(18). Previous findings have shown that Egr-1 and VEGF have key roles in the regulation of growth and development of vessels, and are also critical to ischemic diseases. Additionally, they can facilitate the migration of hematopoietic stem cells with protective effects on nerves (19,20).
In this study, male SD rats used as subjects, were divided the sham, model and G-CSF groups to measure the effect of G-CSF on the volume of cerebral infarction and level of LDH in rats. The H&E staining method was performed for histopathological examination. Through RT-PCR and western blot analysis, we detected the mRNA and protein expressions of Egr-1 and VEGF in different groups. Furthermore, SPSS 17.0 software was applied to detect the differences in expression of Egr-1 and VEGF between the two groups. Results showed that compared with the sham group, the volume of cerebral infarction and LDH content in the model group were significantly elevated, whereas in the model group, those indicators in the G-CSF group were obviously decreased. H&E staining results also showed that G-CSF could decrease the necrotic area in cerebral infarction and the incidence of inflammation, and sustain the integrity of molecular structure. Immunofluorescence staining results revealed that the protein expression of Egr-1 and VEGF in the model group were all significantly decreased, while those in the G-CSF group were significantly elevated. RT-PCR and western blot analysis revealed that the mRNA and protein expression of Egr-1 and VEGF in the model group was decreased obviously, but those in the G-CSF group were elevated significantly.
Taken together, G-CSF manifests a significant protective effect on acute ischemic cerebral injury, which is generated through its impact on the expression of Egr-1 and VEGF. The results of this study are expected to provide new therapeutic procedures for treatment of acute ischemic cerebral injury, and other effects of G-CSF on acute ischemic cerebral injury require more in-depth studies.


s2fieldsofstudy: Biology Medicine

text:
Prediction of novel target genes and pathways involved in tall cell variant papillary thyroid carcinoma

Supplemental Digital Content is available in the text

Introduction
Papillary thyroid carcinoma (PTC) is the most common endocrine neoplasm and the most frequent malignant thyroid tumor. PTC consists of several histological variants, the most common of which are classical/conventional variant PTC (cPTC), follicular-variant PTC (FVPTC), and tall cell variant PTC (TCPTC), which collectively account for the vast majority of PTCs. [1,2] Some rare PTC variants, such as Hobnail, Columnar, Solid and diffuse sclerosing variant PTC, have also been reported. [3,4] Most patients with PTC have excellent prognoses, but some histological variants of PTC, such as TCPTC, are characterized by more aggressive phenotypes. [2,5] TCPTC consists predominantly of tall cells (TC), which have a height of more than at least double their width, as well as an eosinophilic cytoplasm and hyperchromatic basilar nuclei. [1,6] The incidence of TCPTC ranges from 2% to 19% among PTC patients described in the literatures. [5,7,8] The clinical characteristics and prognosis of TCPTC have been reported by a number of studies. The consensus is that TCPTC is related to aggressive clinicopathological parameters, including tumor stage, tumor multifocality, extrathyroidal extension, lymphovascular invasion, initial lymph node metastasis, preablation lung metastasis, and BRAF mutations. [1,5,6,[9][10][11] However, studies on the outcomes of patients with TCPTC have revealed controversial results. [5,7,11,12] Although numerous studies that have investigated the biological characteristics of TCPTC, [13,14] the molecular mechanisms underlying TCPTC remain poorly understood. Therefore, recognizing molecular alterations and identifying novel diagnostic biomarkers for TCPTC is important. A comprehensive and systematic analysis of gene expression profiles in TCPTC is also lacking. The Cancer Genome Atlas (TCGA) database catalogs genetic information, including microRNA and mRNA expression profiles, somatic mutations and copy number variations, and covers 33 types of cancers, thereby facilitating studies on gene expression signatures and tumorigenesis mechanisms. [15] In this study, the gene mutation and mRNA expression profiles of patients with TCPTC from TCGA database were obtained by RNA sequencing. The clinical features and prognoses of TCPTC patients were examined. Differentially expressed genes (DEGs) were identified, and functional enrichment and pathway analyses were performed. Subsequently, protein-protein interaction (PPI) networks were constructed to identify hub genes. We also performed a receiver operating characteristic (ROC) analysis to evaluate the diagnostic values of the candidate genes. Our study may provide a meaningful contribution towards exploring the tumorigenesis mechanism of TCPTC and assisting in the discovery of novel candidate diagnostic markers for TCPTC.

TCGA data access
The standardized level 3 RNA sequencing data of PTC patients and the corresponding clinical records in TCGA were obtained from the cBioPortal for Cancer Genomics (http://www.cbiopor tal.org/index.do) and FireBrowse (http://firebrowse.org). Among the total 507 cases in TCGA thyroid carcinoma (TCGA-THCA) dataset, there were 359 cPTCs, 37 TCPTCs (≥50% tall cell feature) and 59 normal tissue samples. A total of 357 cPTC and 35 TCPTC samples with complete mRNA sequencing and clinical data and 59 normal tissue samples were included in the subsequent analyses (the case IDs are shown in Supplementary Table S1, http://links.lww.com/MD/C715). The RNAseq by Expectation-Maximization (RSEM) values were utilized to quantify the mRNA expression levels. Histological images of all samples are openly available on the Cancer Slide Digital Archive website (CSDA, http://cancer.digitalslidearchive.net).

Identification of DEGs
Differentially expressed mRNAs between the PTC (including cPTC and TCPTC) and normal tissues, and between the TCPTC and cPTC tissues were analyzed with the Limma package of R. The linear model fitting (lmFit) function in the Limma package was used to calculate fold changes, and empirical Bayes statistics (eBayes) were used to estimate standard errors. Only genes with jlog2FCj > 1.0 and P < .05 were considered to be significantly DEGs. The R script used in Limma package is provided in Supplementary File 1, http://links.lww.com/MD/C715. Hierarchical analysis of the identified DEGs was achieved by Cluster 3.0 (calculating average linkage based on Pearson correlation coefficient) and then visualized via the Java TreeView 1.16r4 (http://www.treeview.net).

Functional enrichment analysis
The Database for Annotation, Visualization, and Integrated Discovery (DAVID, https://david.ncifcrf.gov/) online system was utilized for candidate DEGs functions and pathway enrichment analyses. Gene Ontology analysis (GO, classification including biological process (BP), cellular component (CC) and molecular function (MF)) and Kyoto Encyclopedia of Genes and Genomes (KEGG) pathway enrichment analyses were performed. The cutoff was set at a P value < .05.

Protein-protein interaction network analysis
The PPI network was constructed using the Search Tool for the Retrieval of Interacting Genes/Proteins (STRING, https://stringdb.org, the minimum required interaction score was set at 0.4) database, and the results were visualized with Cytoscape software 3.5.1. Subnetwork models were selected using the plugin molecular complex detection (MCODE) application in Cytoscape. The important coexpression modules of DEGs were detected from the PPI network, and the criteria for defining a module were an MCODE score ≥5 and a number of nodes ≥4. [16] Functional enrichment of the genes in the subnetworks was performed with STRING. Hub genes were explored using the CytoHubba plugin application with the topological analysis of Maximal Clique Centrality (MCC, top 10 nodes ranked by MCC). [17] Cytohubba is a tool that provides 11 topological analysis and 6 centralities based on shortest paths to detect the hub genes (key genes). [18]

Statistical analysis
Statistical analysis was carried out using SPSS 21.0 and R studio. The Chi-squared or Fisher exact test was used to analyze the associations between different histological diagnoses and clinicopathological parameters. Kaplan-Meier and log-rank tests were applied to illustrate the risk of disease-free and overall survival rates between patients with different pathologic types. A P value < .05 was considered statistically significant. ROC curves and the area under the curve (AUC) were determined to test the specificity and sensitivity of each candidate gene for diagnostic predictions. ROC analysis was used to find a cut-off point of each hub gene for distinguishing different pathologic types.

Somatic mutations in TCPTC and cPTC
The gene mutation frequencies in patients with TCPTC and cPTC were obtained by selecting relevant case IDs through the cBioPortal website. In total, altered genes (mutations in more than 2 cases) occurred in 29 of the 37 TCPTC cases. The top 10 gene mutation profiles are shown in Figure 1. Alterations in BRAF (28 BRAF V600E mutations and 1 deletion mutation) were particularly prevalent in TCPTC with a mutation frequency of 78%. Patients with the other 9 gene mutations also harbored BRAF mutations. None of the total 37 patients harbored RAS gene mutations. BRAF mutation was also the predominant mutation type in cPTC (occurring in 56% of cPTC patients, Supplementary Fig. S1, http://links.lww.com/MD/C715), fol- lowed by RET fusion (with a relatively low occurrence rate of 7%) and RAS mutations (including NRAS, HRAS and KRAS mutations with frequencies of 4%, 1.7%, and 0.71%, respectively).

TCPTC was associated with aggressive clinicopathological parameters
To address the clinical differences between TCPTC and cPTC, the correlations of the two histological diagnoses with clinicopathological parameters were examined. The results showed that TCPTC was significantly correlated with a patient age >45 years, tumor multifocality, extrathyroidal extension, a higher T stage, advanced AJCC TNM stages, and BRAF V600E mutation. No significance was observed for patient gender, N stage or M stage (Table 1). Additionally, Kaplan-Meier analysis revealed that patients with TCPTC had a poorer disease-free survival rate, whereas no difference in the overall survival rate was evident between the TCPTC and cPTC patients (Fig. 2). These results suggested that TCPTC patients had more aggressive clinicopathological features and poorer outcomes than cPTC patients.  normal tissues, except for 7 genes) were much higher in the TCPTC tissues than those in the cPTC tissues. The expression levels of nearly all putative tumor suppressors (genes downregulated both in TCPTC and cPTC compared with normal tissues, except for 1 gene) were much lower in the TCPTC tissues than those in the cPTC tissues.

GO and KEGG pathway enrichment analysis of TCPTC-related and TCPTC-specific DEGs
To gain a better understanding of the gene functions and signaling pathways of TCPTC-related and TCPTC-specific DEGs, online GO and KEGG pathway enrichment analyses were conducted using DAVID. Due to gene quantitative restriction, we enriched the upregulated genes and the downregulated TCPTC-related genes separately. For the TCPTC-related upregulated DEGs, the GO analysis results showed that the DEGs were significantly enriched in cell adhesion at the BP level, in integral component of membrane at the CC level and in calcium ion binding at the MF level. The top 10 enriched terms are shown in Figure 5A. Furthermore, the enriched KEGG pathways of the TCPTC-related upregulated DEGs included pathways in cancer, cytokine-cytokine receptor interactions, HTLV-I infection, PI3K-Akt signaling, and cell adhesion (Fig. 5B). For the TCPTC-related downregulated DEGs, the GO analysis results showed that the DEGs were significantly enriched in chemical synaptic transmission at the BP level, plasma membrane at the CC level, and calcium ion binding at the MF level (Fig. 5C). Furthermore, the enriched KEGG pathways of the TCPTC-related downregulated DEGs included pathways in cancer, neuroactive ligand-receptor interactions, PI3K-Akt signaling, MAPK signaling, and cAMP signaling (Fig. 5D). For the TCPTC-specific DEGs, the top enriched GO terms for BP, CC and MF were extracellular matrix organization, extracellular space, and extracellular matrix structural constituent, respectively (Fig. 5E). The enriched KEGG pathways of the TCPTC-specific DEGs included pathways in cytokine-cytokine receptor interactions, ECM-receptor interaction, neuroactive ligand-receptor interaction, cell adhesion molecules, and focal adhesions (Fig. 5F).

PPI network construction and hub gene selection
To identify TCPTC-specific biomarkers, a PPI network was constructed with only 301 TCPTC-specific DEGs by the STRING database and then visualized with Cytoscape software. The PPI network consisted of 166 nodes and 321 edges (Fig. 6A). A total of 3 subnetworks, comprising 25 nodes and 76 edges, were found using the defined criteria (Fig. 6B). GO enrichment revealed that these genes were annotated in extracellular matrix organization at the BP level, extracellular space at the CC level, and extracellular matrix structural constituent at the MF level. The enriched KEGG pathways included pathways in cytokinecytokine receptor interactions, ECM-receptor interactions, chemokine signaling, focal adhesion, and PI3K-Akt signaling. All enriched GO and KEGG terms are shown in Supplementary

Diagnostic value of hub genes in TCPTC
ROC analyses were performed to evaluate the specificities and sensitivities of the hub genes for diagnostic predictions. The diagnostic performances were carried out in 3 independent groups (namely, TCPTC vs. normal tissues, cPTC vs. normal tissues, and TCPTC vs. cPTC tissues). To determine the discriminatory abilities of these candidate genes in each group, an AUC > 0.7 was fixed as the threshold. The ROC curve analysis demonstrated that all 8 hub genes yielded high diagnostic accuracy in discriminating TCPTC from normal tissues (Fig. 7A). Similarly, except for CCR8, the aforementioned hub genes could distinguish cPTC from normal tissues (Fig. 7B). Furthermore, all 7 upregulated hub genes could robustly distinguish TCPTC from cPTC tissues (Fig. 7C). Collectively, 6 hub genes (COL5A1, COL1A1, COL10A1, COL11A1, CCL20, and CXCL5) were identified in all 3 independent groups. The AUCs and cut-off values of these 6 genes for distinguishing the 3 independent groups are presented in Table 2. Our results suggest that a combination of these hub genes can be used not only for the differential diagnosis of PTC, including TCPTC and cPTC, from normal samples, but also for the differential diagnosis of TCPTC from cPTC samples.

Discussion
In this study, we studied the gene mutation patterns and clinical characteristics of TCPTC by utilizing data obtained from TCGA-THCA dataset. Multi-group analyses were performed to identify aberrantly expressed TCPTC-related and TCPTC-specific genes, and these analyses were followed by functional and pathway enrichment analyses. Importantly, we highlighted the diagnostic values of several hub genes not only for distinguishing TCPTC and cPTC from normal tissues, but also for distinguishing TCPTC from cPTC tissues. By analyzing the TCGA-THCA data, Stokowy et al revealed molecular differences, including micro-RNA and mRNA differences, between TCPTC and FVPTC. [19] To the best of our knowledge, this investigation is the first global study of gene mutations, screening for TCPTC-specific DEGs, and diagnostic approaches for candidate biomarkers of TCPTC by utilizing TCGA-THCA data. We found that BRAF V600E was the predominant gene mutation type with a high incidence rate (78%), whereas other mutation types were rare in patients with TCPTC. BRAF V600E was also the predominant gene mutation type, and other gene mutation types, including RET fusion and RAS mutations, were also common in patients with cPTC. Consistent with other reports, BRAF V600E-positive PTC was often the conventional or tall cell variant. [20] BRAF mutation was more frequent in TCPTC than that in cPTC, and more frequent in tall cell variant papillary thyroid microcarcinoma (PTMC) than that in classic PTMC. [21,22] Some rare gene mutations, such as the Titin (TTN) mutation, have been observed in TCPTC. A truncating TTN mutation has been identified as a cause of congenital myopathy that is reported as centronuclear myopathy, [23] but no evidence has shown that the TTN mutation is associated with cancers. Some rare gene mutation types that were not identified in our study, including the RET/PTC rearrangement and TERT promoter mutation were also reported to be prevalent in TCPTC. [24,25] We disclosed that TCPTC was associated with more aggressive clinicopathological parameters and poorer outcomes than cPTC. This finding is consistent with the results of the other studies. The overall survival analysis showed a better survival trend of patients with TCPTC compared to those with cPTC. This may be related to the number of TCPTC cases and the duration of follow-up time. A large multicenter study launched by Shi et al demonstrated the differential prognostic risk of the 3 major PTC variants and established a unique risk order of TCPTC > cPTC >> FVPTC. [1] A recent meta-analysis also revealed that TCPTC was associated with more aggressive clinicopathological characteristics and poorer prognoses. [10] The percentage of TC necessary to diagnose TCPTC has been debated. Beninato et al revealed that the aggressive features conferred by the presence of TC in PTCs occurred with as little as 10% TC and were maintained with increasing percentages. [26] Oh et al revealed that PTC with 10% to 50% TC showed a similar BRAF mutation rate and clinicopathological features to TCPTC. [27] Ganly et al suggested that consideration should be given to using a 30% TC threshold to diagnose TCPTC, [11] whereas Ito et al suggested that diagnosing TCPTC with a TC threshold ≥50% is appropriate. [28] The TC percentage was ≥50% in all of the TCPTC patients included in this study. These Matriptase-2 (TMPRSS6) is a member of the type 2 transmembrane serine protease family. TMPRSS6 has been implicated in the progression of cancers, including breast, prostate, and colorectal cancers. [29] Strong TMRRSS6 expression has been detected in thyroid tissues, [30] which may indicate that TMRRSS6 plays a specialized role in PTC. COL11A1 dysregulation across cancer types is particularly striking, including colorectal, ovarian, breast, head and neck, lung, and brain cancers. [31] This result indicates that COL11A1 may serve as a remarkable biomarker for various types of cancers and as a target for cancer therapy. Park et al revealed that the COL11A1 gene might be associated with PTC and, in particular, that the T allele of rs1763347 and rs2229783 might contribute to reducing the risk of PTC. [32] Carbonic anhydrase 4 (CA4) is a member of the carbonic anhydrases family. Davidov conducted a tissue microarray study that comprised 26 follicular thyroid carcinomas and 53 follicular adenomas from patients with indeterminate thyroid nodules. By analyzing the staining results for 17 immunohistochemical biomarkers, the authors revealed that loss of immunoreactivity of CA4 was associated with malignancy in indeterminate thyroid specimens. [33] These genes are the topmost dysregulated genes among the TCPTC-related and TCPTCspecific DEGs. However, the functions and specific mechanisms of these genes need elucidation. GO and KEGG pathway analyses of the DEGs were performed to understand the gene functions and interactions of the TCPTC-related and TCPTC-specific DEGs. The enriched KEGG pathways of the TCPTC-related DEGs included typical pathways involved in thyroid cancer, such as the MAPK and PI3K-Akt signaling pathways. [34][35][36][37] These data suggest that TCPTC shares some same signaling pathway abnormalities with cPTC. Intriguingly, the results of our study showed that the dysregulated levels of both putative oncogenes and tumor suppressors in TCPTC were higher than those in cPTC. Furthermore, the higher degree of abnormality of the enriched functions and signaling pathways associated with the TCPTC-specific DEGs may be one reason for the higher degree of malignancy of TCPTC.
TCPTC is usually associated with a poor outcome, and the correct pre-surgery recognition of this variant is important for clinical management. Evranos et al conducted a retrospective study to compare the malignancy rates determined by fine needle aspiration (FNA) and histopathological results. Pre-surgical cytology showed malignancy or suspicion of malignancy in 92% of the aggressive (including TCPTC) variants of histopathologically confirmed PTC. [38] These results suggest that the Bethesda classification is a reliable indicator of the malignancy of nodules with aggressive variant PTC. However, other studies have revealed that pathologists often face the dilemma of a proper diagnosis of TCPTC, based on cytology as well as histology. [13] A   Figure 7. ROC curves for the diagnostic performances of the hub genes. A. Diagnostic performance of the hub genes for differentiating TCPTC from normal tissues. All 8 hub genes showed high diagnostic accuracy. B. Diagnostic performance of the hub genes for differentiating cPTC from normal tissues. Seven hub genes could be used to distinguish cPTC from normal tissues. C. Diagnostic performance of the hub genes for differentiating TCPTC from cPTC tissues. All 7 upregulated hub genes showed high diagnostic accuracy. cPTC = classical/conventional variant Papillary Thyroid Carcinoma, ROC = receiver operating characteristic, TCPTC = tall cell variant papillary thyroid carcinoma. Table 2 Cut-off values and AUCs of 6 genes for distinguishing the 3 independent groups.  few studies have described the identification of biomarkers to improve the differential diagnosis of TCPTC from cPTC lesions. [39,40] For example, Nardone et al have revealed that TCPTC specimens expressed significantly greater levels of c-Met than other forms of PTC and benign thyroid disease. [40] In this study, we revealed for the first time that several biomarkers showed excellent predictive values not only for the differential diagnosis of TCPTC from normal tissues but also for the differential diagnosis of TCPTC from cPTC tissues. Future studies with more cost-effective approaches through FNA or postoperative histopathologies, such as RT-PCR and immunohistochemical staining, should be conducted instead of RNA sequencing to validate the clinical utility of the individual or combined candidate biomarkers identified in our study.
In summary, we have identified notable alterations of TCPTC, including gene mutations, DEGs and related pathways, and evaluated the diagnostic performances of hub genes for TCPTC. We disclosed that TCPTC was associated with more aggressive clinicopathological parameters and poorer outcomes than cPTC. This could be partly attributed to the high incidence rate of BRAF mutation. Our work may provide additional diagnostic targets for TCPTC. The use of molecular profiles for differential diagnoses of TCPTC from cPTC should be evaluated by studies with larger patient cohorts. More work is needed to elucidate the functional mechanisms of these genes and to validate the clinical utility of these biomarkers with more cost-effective approaches.

Author contributions
Xia F and Li X. conception and design; All authors. methodology and development; Xia F, Jiang B, Chen Y, and Li X. data acquisition and analysis; Xia F and Li X. manuscript preparation and editing.