BACKGROUND
proteins carry out their biological function within intricate networks of interacting molecules.
high throughput techniques such as whole genome sequencing have led to accurate representation of gene networks  <cit> .
key insights regarding the influence of natural selection on genes can be obtained by taking into account the network topology  <cit> .
the evolutionary rate of a protein , represented by the ratio of the rate of its nonsynonymous substitutions  to the rate of its synonymous substitutions , is used as an indicator of the selective constraints acting on proteins.
studies have indicated that many factors affect the evolutionary rate of genes within gene pathways and networks.
for example, gene position  <cit> , protein length  <cit> , the number of protein-protein interactions  <cit>  and codon bias  <cit>  have influence on gene evolution.
signaling pathways mediate the sensing and processing of both extracellular and intracellular stimuli.
they rely on receptors that recognize a signaling molecule and trigger a series of events leading to the transmission of signals to the downstream region of the pathway.
within the drosophila toll and imd signaling pathways, downstream genes were more conserved, indicating a relatively stronger evolutionary constraints than upstream genes  <cit> .
this observation is consistent with trends reported for the yeast hog-signaling pathway  <cit>  and the caenorhabditis elegans and drosophila insulin/tor pathway  <cit> .
in contrast, evolutionary patterns of the insulin/foxo signaling pathway across metazoan species genomes indicated that components within the middle of the pathway were rather under stronger purifying selection  <cit> .
moreover, within the drosophila ras signaling pathway, upstream genes have lower evolutionary rate, suggesting more constraints, than downstream genes  <cit> .
obviously, different patterns of selective constraints exist for different signaling pathways.
this is expected, as signaling inputs and network architecture vary among the pathways  <cit> .
the tlr signaling pathway represents the best characterized component of the innate immune system  <cit> .
tlr signaling involves the dimerization of receptors in response to pathogenic microbial products.
this is followed by the recruitment of various adaptor molecules.
the recruitment of these adaptor molecules leads to the recruitment of downstream signaling molecules that activate transcription factors leading to the release of inflammatory cytokines and type i interferons.
the tlr signaling is therefore a sequential process.
thus, the effects of pathway parameters such as gene position on the evolutionary rate of genes warrant further analyses.
the tlr signaling pathway has been selectively constrained over time, indicating its essential role for host survival  <cit> .
in humans, population genetic studies have indicated that purifying selection is stronger within tlr adaptors relative to receptors  <cit> .
a recent study on the tlr signaling pathway-related genes from eight vertebrate genomes showed that the selective constraints of genes was negatively correlated with gene position along the tlr signaling pathway  <cit> .
thus, different components of the tlr signaling pathway appear to differ in their evolutionary rates.
moreover, taxon specific differences among pathway components have been demonstrated for the tlr pathway  <cit> .
new insights relating to the influence of network parameters on the evolutionary rate of genes within the tlr signaling pathway can be gained through studies on closely related species not yet investigated in detail.
members of the family suidae have successfully inhabited broad geographic locations including eurasia, mainland southeast asia, island southeast asia  <cit>  and the african continent  <cit> .
given the numerous pathogenic challenges across these diverse environments, the tlr signaling pathway may have been critical in the survival of suids.
suid species are increasingly becoming important biomedical models given their anatomical, physiological and immunogenetics similarities with humans.
comparison of protein sequences of the domestic pig, a member of the suidae, to human orthologues revealed  <dig> positions where the domestic pig protein has the same amino acid that is implicated in a human disease  <cit> .
thus, findings regarding suid tlr signaling genes evolutionary rates can be relevant to their hominid orthologues.
the availability of a complete porcine genome and resequencing libraries of species within the family suidae provide an opportunity to investigate the evolution of the tlr signaling pathway genes within the context of gene position and other network parameters.
the suid species have evolved over a relative short time span of 1– <dig> million years  <cit> .
the low divergence among these species indicate that synonymous substitutions are not likely to be saturated, making the estimation and comparison of selective pressure  among genes less prone to bias.
this low divergence among the species also ensures reliability in aligning tlr sequences for subsequent analysis.
in the current study, our goal was to investigate the evolutionary constraints of genes within the tlr signaling pathway of the family suidae.
we sought to investigate the extent to which patterns of protein evolution observed in other organisms extend to members of the family suidae.
we hypothesized that there is a polarity of purifying selection along the tlr signaling pathway within the family suidae.
we tested our hypothesis by answering the following questions 1) is there a relationship between the strength of purifying selection and gene position in the tlr signaling pathway and; 2) are there any network parameters that contribute to the polarity in the strength of purifying selection?
methods
ethics statement
blood samples were obtained from captive-held suid species by qualified individuals at each institution, according to their respective animal welfare oversight governance and institutional permits .
all animal procedures were conducted in accordance with the guidelines approved by the american society of mammalogists  <cit> , and approved by the university of illinois’ institutional animal care and use committee .
sources of animals are provided in additional file 1: table s <dig>
dna extraction and sequencing
dna was extracted from whole blood using the qiaamp dna blood spin kit  and quantity and quality parameters were performed using the qubit  <dig>  fluorometer  and visualized using a 1 % agarose gel.
library construction and re-sequencing of individual members of the family suidae were done with 1–3 μg of genomic dna according to the illumina library prepping protocols.
the library insert size was 300–500 bp and sequencing was performed using a  <dig> paired-end sequencing kit.
all dna samples were sequenced to approximately 8x depth.
quality trimmed reads  were aligned to the sus scrofa reference genome build  <dig>  using the unique alignment option of mosaik aligner .
the aligned reads from each of the animals  together with the sus scrofa reference genome were stored as bam files for each individual animal.
ortholog identification
genes involved in the tlr signaling pathway  were examined for this study .
in addition tlr <dig>  which is not found in the kegg database tlr signaling pathway, was included since it cooperates with tlr <dig>  senses lipopeptides  <cit>  and activates the tlr signaling pathway through association with myd <dig>  <cit> .
the genomic coordinates of sus scrofa genes within the tlr signaling pathway  were obtained from the ensembl database .
based on these genomic coordinates, sequences of gene orthologs were then retrieved from aligned bam files  for sus scrofa , sus verrucosus, sus celebensis, sus scebifrons, sus barbatus, babyrousa babyrussa, potamochoerus larvatus, potamochoerus porcus and phacochoerus africanus to identify gene orthologs.
the resulting sequences for each species were then screened against the genome of sus scrofa using blast to ensure similarity with the sus scrofa genes.
to obtain the coding sequence of genes for each species, exonic regions were retrieved based on genomic coordinates of exons from the sus scrofa gene transcripts  and concatenated.
if a gene was found to have more than one transcript, the longest transcript was chosen for analysis.
coding sequences were aligned using clustalw  <dig>   <cit> .
aligned sequences are provided in additional file  <dig> fig.
<dig> the tlr signaling pathway genes.
redrawn following  <cit> , with modification to include tlr <dig>  the direction of signal transduction is indicated by the arrows.
tlrs  <dig>   <dig>   <dig>  and  <dig> are included in the first pathway position category.
the numbers on the left side represent the position of the pathway genes.
genes used in this study are shown in bold
some of the coding sequences had missing nucleotides.
the following criteria were used to discard such sequences:  sequences with more than 50 % of missing amino acid residues; and  aligned sequences absent in any one of the species involved in this study.
consequently, the final data set was composed of  <dig> genes .
impact of natural selection
the impact of natural selection on genes within the tlr signaling pathway was determined by estimating nonsynonymous substitution rate , synonymous substitution rate  and their ratio  using the m <dig> model implemented in the program codeml from the paml package version  <dig>  <cit> .
to determine whether specific codon positions were under positive selection within gene orthologs, codeml site model m1a, a nearly neutral evolution model where sites are assumed to be evolving under either purifying selection  or neutral evolution  was compared to model m2a that allows positive selection  among sites.
m <dig>  which allows sites to evolve under either purifying selection or neutrally, was compared to model m <dig>  which allows for positively selected sites  .
models were compared using a likelihood ratio test .
the f3x <dig> model of codon frequencies was used for the analyses.
models were run in duplicates with ω of  <dig>  and  <dig>  to increase the probability of convergence of model parameters.
multiple testing for positive selection on genes was corrected for by conducting a false discovery rate  test  <cit>  at a q value of  <dig> .
the phylogenetic relationship among the species of suidae  inferred from near complete genome sequences of each species  <cit>  was used for the analysis involving the codeml program.fig.
<dig> phylogenetic relationships among species of the family suidae obtained from near complete genome data of each species.
the posterior probability at each node is 1
network level analysis
to determine whether the evolution of molecules within each gene was affected by network structure, network parameters were computed for each gene and were correlated with each evolutionary parameters estimated by model m <dig>  the number of protein-protein interactions  for proteins encoded by each gene was determined from the sus scrofa interaction network in the string database .
string is a database of known and predicted interactions which include direct  and indirect  associations from various sources including high throughput experiments and genomic contexts  <cit> .
the magnitude for codon bias for individual genes was measured by the effective number of codons .
the codon usage bias of each orthologous gene group was measured as the mean of enc of each species gene.
enc values of each gene were obtained using dnasp software  <cit> .
bivariate correlation analysis between network parameters ) and length of 3′utr ) and evolutionary parameters  were conducted using spearman’s rank correlations.
a false discovery rate  test  <cit>  was performed to correct for multiple testing for correlations controlling for q value at  <dig> .
bivariate correlations and corrections for multiple testing were performed using sas software .
figures illustrating bivariate correlations were obtained with the package r .
multivariate analyses
the two multivariate analyses methods  were used to determine whether the observed bivariate correlations were due to direct or indirect influences.
partial correlation analysis measures the strength of relationship between two variables, while holding one or more variables constant.
path analysis estimates direct and indirect relationships under a user defined causal model.
for the path analysis, pathway position, ppi, protein length, enc and l3′utr were considered as exogenous variables whereas ω and dn were considered as endogenous variables.
prior to performing the path analysis, data were log-transformed to improve normality.
all statistics related to path analysis were calculated by the lavaan package in r.
RESULTS
orthologs of  <dig> sus scrofa tlr signaling pathway genes were identified for  <dig> members of the family suidae.
a total of  <dig> sequences of the tlr signaling pathway genes ranging from  <dig> to  <dig> nucleotides were therefore used in analysis.
the orthologous sequences were aligned for a series of analyses including test for selection pressures acting on genes, bivariate correlations and multivariate analysis.
analysis of protein sequence evolution
in order to estimate the selective pressure acting on genes within the suid tlr signaling pathway, the m <dig> model, which provides a single estimate of ω across all codons and lineages, was utilized.
the ω values ranged from  <dig>   to  <dig>   .
the mean ω was  <dig> .
to test for positive selection, the models m1a  vs m2a ; m <dig>  vs m <dig>  were used.
the two tests were in agreement for all genes  and demonstrated that a proportion of sites within tlr <dig>  tlr <dig>  tlr <dig> and irak <dig>  were under positive selection.
however, tlr <dig> was the only gene that had sites under positive selection after correcting for multiple testing at q =  <dig> .
thus, results here indicated that genes within the tlr signaling pathway have evolved under strong functional constraint.
relationship between evolutionary rates and pathway position
the relationship between the evolutionary parameter ω of genes and their position within the tlr signaling pathway was determined using a spearman’s rank correlation test between the two variables.
the ω values were negatively correlated with pathway position  indicating that downstream genes are under stronger purifying selection than upstream genes.
the robustness of the correlation between ω and pathway position was tested by removing the highest ω value  and re-estimating the correlation.
this was done to ensure that results were not influenced by high ω values.
the correlation remained significant .
to test for the evolutionary parameter accounting for the ω differences among genes, the correlation between pathway position and dn, pathway position and ds was also determined.
the dn was negatively correlated with the pathway position .
there was no significant correlation between ds and pathway position .
these results indicated that the decrease in ω  along the tlr signaling pathway  is attributable to a decrease in dn .fig.
<dig> ω and dn versus pathway position, number of protein-protein interactions  and protein length.
all relationships are significant.
continuous lines represent regression lines
relationship between evolutionary rates and network parameters
network parameters other than pathway position may influence evolutionary parameters  and account for inferred polarity in ω and dn.
in order to determine such relationships, spearman’s rank correlation test was performed between ppi of proteins  encoded by genes within the tlr signaling pathway, protein length, codon bias measured as effective number of codons , l3′utr region and the evolutionary parameters .
the ω values were negatively correlated with ppi  and positively correlated with protein length .
similarly, dn values were negatively correlated with ppi  and positively correlated with protein length  indicating that apart from pathway position, the number of protein-protein interactions a given protein is involved in and its protein length influence it evolution.
in addition, ds was also negatively correlated with enc  indicating stronger selection based on codon usage in genes with high codon bias than genes with low codon bias  <cit> .
since pathway position, ppi and protein length are intercorrelated and are each correlated with evolutionary parameters  , observed associations of these network parameters with the evolutionary parameters could be indirect .
thus, to distinguish between direct and indirect effects, multivariate analysis  were utilized.
partial correlation analysis revealed that when controlling for ppi, the correlation between ω and pathway position remained significant .
when controlling for protein length, the correlation between ω and pathway position remained significant .
thus, gene position within the tlr signaling pathway directly influenced the evolution of genes in the tlr signaling pathway.
the correlation between ω and ppi while controlling for pathway position remained significant .
it also remained significant  when controlling for protein length.
this result indicated that ppi is an important factor affecting the evolution of genes within the tlr signaling pathway.
the correlation between ω and protein length was not significant  when pathway position was held constant indicating that the correlation between ω and protein length was rather mediated by pathway position.
nonsynonymous  substitution rates were not significantly  correlated with pathway position when controlling for protein length.
in contrast, the correlation between dn and ppi remained significant  when controlling for protein length.
if pathway position was held constant, the correlation between dn and ppi remained significant .
correlation between dn and protein length was still significant  after controlling for pathway position and remained significant  if controlling for ppi.
results for the relationship between dn values and network parameters indicated that the number of protein-protein interactions and protein length have direct effects on the rate of dn evolution.
the relationship between evolutionary parameters  and the network parameters  were further analyzed using path analysis.
as depicted in fig.
<dig>  path analysis indicated that dn values were affected by gene position within the tlr signaling pathway .
pathway position had the largest direct effect on dn.
the dn values were negatively associated with the number of protein-protein interactions  and positively associated with protein length .
ω values were only associated with enc .
thus pathway position, the number of protein interactions and protein length were factors that influenced tlr signaling pathway dn substitution rates in the suidae.fig.
<dig> the causal model used to analyze the relationship among evolutionary and network parameters.
the network parameter pathway position, number of protein-protein interactions , protein length, codon bias ) and length of the 3′utr region are considered as exogenous variables.
continuous and dashed lines represent significant and nonsignificant relationships, respectively.
single-headed arrows indicate causal relationship between variables.
double-headed arrows indicate correlations between exogenous variables.
numbers on the arrows represent the standardized regression weights
discussion
our results support the hypothesis that there is polarity of purifying selection along the tlr signaling pathway  within the family suidae with purifying selective pressure increasing along the pathway.
protein-protein interactions and protein length accounted for this polarity in the strength of purifying selection, where number of interactors of tlr signaling protein molecules increases and protein length decreasing along the tlr signaling pathway.
the evolutionary parameter influenced by these network parameters was dn , which is actually the metric of selective pressure  <cit> .
fewer amino acid substitutions were tolerated along the signaling pathway reflective of stronger purifying selection.
in a previous study on the evolution of tlr signaling pathway-related genes from eight vertebrate genomes   <cit> , selective constraint was negatively correlated with gene position in the signaling pathway and dn was also influenced by position and protein length in a similar manner as in the present study.
these similarities indicate that some patterns revealed across the eight vertebrate genomes are also seen in the closely related suid species involved in this study.
several other studies conducted on various biological pathways have revealed a polarity in the strength of purifying selection.
upstream genes revealed higher constrains than downstream genes in the carotenoid biosynthesis pathway  <cit> , in the plant terpenoid biosynthesis pathway  <cit>  and in the dopamine catabolic pathway across mammals  <cit> .
furthermore, in the drosophila ras signaling pathway, upstream genes evolve under higher levels of selective constrains relative to their downstream counterparts  <cit>  the relatively slow rate of evolution of upstream genes reflects these genes being required for a wider range of end products and therefore being more pleiotropic  <cit>  and in the case of biosynthesis pathways, upstream genes exert a higher control over the flux of metabolites throughout the pathway  <cit> .
a negative correlation between strength of purifying selection and gene position within pathway, as inferred in this study, have also been reported for the drosophila toll and imd signaling pathways  <cit> , caenorhabditis insulin/tor signal transduction pathway  <cit>  and the yeast hog pathway  <cit> .
this pattern of selective constraint might be attributed to purifying selection acting to maintain the function of downstream signal transduction elements and a concentration of adaptive changes in the upstream genes due to their interaction with external environment  <cit> .
in this study, tlr <dig> was the only gene inferred to be under positive selection after fdr correction.
therefore, adaptive changes were not prominent within upstream genes of the tlr signaling pathway to result in the inferred pattern of selective constraints.
with the inferred significantly negative correlation between pathway position and the evolutionary parameters ω and dn, a negative correlation between pathway position and number of protein-protein interactions for a protein, as indicated in this study, implies a role for the number of protein-protein interactions in the polarity of selective constraints.
this negative relationship is indicative of increasing selective constraints for proteins along the tlr signaling pathway, as a result of increase in the number of interactors along the signaling pathway.
evidence for higher connectivity  of downstream genes has also been demonstrated for the human signal transduction network  <cit> .
within metabolic pathways, connectivity plays a major role in constraining evolutionary rates with proteins interacting with more proteins being subject to stronger selective constraints  <cit> .
the increasing number of interactors along the tlr signaling pathway indicates the essential role of protein molecules downstream of tlrs in propagating and amplifying signals initiated by pathogen binding of tlrs and the role of downstream protein molecules in other functions apart from tlr signaling.
this places stronger selective constraint on the intracellular signaling protein molecules.
the finding of an effect of connectivity on dn in this study is in contrast to the study of song et al.
<cit>  that did not show a significant influence of connectivity on dn in a multivariate analysis of factors influencing evolution of genes in the tlr signaling pathway across  <dig> vertebrate genomes.
a possible explanation for this difference may be that the use of suid species that are more closely related in this study than species involved in the study of song et al.
<cit>  led to greater precision in the estimation of dn and thus the relationship between connectivity and dn was more evident.
the finding of a significant negative correlation between protein length and pathway position in this study corroborates studies in the insulin/tor signaling of drosophila  <cit>  and indicated that downstream genes encode shorter and more actively translated proteins  <cit> .
given that partial correlation and pathway analysis were both in agreement as to a positive relationship between protein length and dn, protein length appears to be a factor responsible for the polarity in dn along tlr signaling pathway.
thus, along the tlr signaling pathway, the decrease in protein length is associated with decreases in amino acid substitutions.
a similar pattern has been reported for the woody perennial plant populus tremula, where protein length is the main factor affecting selective constraints, with purifying selection weaker in genes with longer coding regions  <cit> .
a possible explanation for the relationship between purifying selection and protein length is that selection at more than one site should cause an overall reduction in the effectiveness of selection   <cit> .
in that case, for longer proteins, which may have many sites under selection simultaneously  <cit> , there will be a reduced efficiency in natural selection.
though path analysis indicated a negative association between ω and enc, pathway position and enc were not significantly correlated, therefore enc cannot explain correlation between selective constraints and enc.
in this study, length of the gene 3utr region did not have a relationship with any evolutionary parameter in disagreement with studies in the toll/imd pathway in the drosophila toll and imd signaling pathways  <cit>  and human and mouse mirna target prediction data  <cit> .
genes within the signaling pathway had an average ω value <  <dig>  indicating that the tlr signaling pathway is selectively constrained across members of the suidae.
this is in keeping with the essential role of the pathway in innate immunity and host survival.
selectively constrained regions within the genome are likely to be functionally important  <cit> .
the strongest selective constraints inferred for the map kinases  with dn substitution rates of zero for these genes suggest their involvement in critical roles during tlr signaling.
the map kinases cascade components play important role in the production of proinflammatory mediators  <cit> .
recently, it has been demonstrated that map kinases play a role in agonist dependent regulation of cognate tlr mrna levels  <cit> .
map kinases are required for the regulation of cellular development and differentiation processes  <cit> .
the conservation of the map kinase genes in this study could therefore be attributed to their involvement in many processes.
the stronger purifying selection detected in downstream genes of tlrs relative to upstream genes suggests that there is a greater need to protect the integrity of proteins downstream  <cit> .
thus, genes downstream might be essential for survival of suid species involved in this study.
tlr signaling evolution is partly dictated by the need to conserve the capability to recognize specific pathogen signature  <cit> .
for the homogenous group of suid species involved in this study, the evolutionary dynamics of the tlr signaling pathway as revealed in this study may be reflective of the range of related pathogens that adversely affected them during their evolution.
our study has provided insight into factors that have shaped the present day tlr signaling pathway genes of the suid species, a convenient number of which are still present in the world  <cit> .
as the wild suid species are among animals at the greatest risk of declines due to cross species transmission of disease with their domesticated counterparts   <cit> , evolutionary analysis of the tlr signaling pathway as presented here can aid in identifying genes that are essential in mediating the immune response to disease causing bacteria and viruses.
for example, a protein that has a higher connectivity and selectively constrained within the tlr signaling pathway would be expected to be of functional relevance in terms of disease prevention.
CONCLUSIONS
the findings here provide evidence of the role of network topology in the polarity of purifying selection along the tlr signaling pathway of the suidae.
our work indicates that the suidae are an attractive system of closely related but well characterized species to study the evolutionary dynamics of tlr signaling genes.
furthermore, our work opens other avenues of research including the use of the suidae in understanding patterns of molecular evolution in other important pathways such as metabolic and gene-regulation pathways.
availability of supporting data
all supporting data are included as additional files in the form of additional files  <dig>   <dig>   <dig>   <dig>   <dig> and  <dig>
additional files
additional file 1: table s <dig>  sample information.
additional file 2: table s <dig>  genes of the sus scrofa tlr signaling pathway used in querying genomes of other species of the suidae.
additional file 3:
multiple sequence alignments of the tlr signaling pathway orthologs across the suid species.
asteriks indicate same nucleotide and “-“ represents missing nucleotide.
additional file 4: table s <dig>  summary statistics for genes.
additional file 5: table s <dig>  results of codon based test for positive selection.
additional file 6: table s <dig>  bivariate correlations among variables.
competing interests
the authors declare that they have no competing interests.
authors’ contributions
conceived and designed the experiments: kado mamg lbs.
performed the experiments: kado.
analyzed the data: kado hjm alr mamg lbs.
contributed reagents/materials/analysis tools: hjm alr mamg lbs.
wrote the paper: kado lbs.
all authors read and approved the final manuscript.

 proteins carry out their biological function within intricate networks of interacting molecules. high throughput techniques such as whole genome sequencing have led to accurate representation of gene networks. key insights regarding the influence of natural selection on genes can be obtained by taking into account the network topology. the evolutionary rate of a protein, represented by the ratio of the rate of its nonsynonymous substitutions to the rate of its synonymous substitutions, is used as an indicator of the selective constraints acting on proteins. studies have indicated that many factors affect the evolutionary rate of genes within gene pathways and networks. for example, gene position, protein length, the number of protein-protein interactions and codon bias have influence on gene evolution. signaling pathways mediate the sensing and processing of both extracellular and intracellular stimuli. they rely on receptors that recognize a signaling molecule and trigger a series of events leading to the transmission of signals to the downstream region of the pathway. within the drosophila toll and imd signaling pathways, downstream genes were more conserved, indicating a relatively stronger evolutionary constraints than upstream genes. this observation is consistent with trends reported for the yeast hog-signaling pathway and the caenorhabditis elegans and drosophila insulin/tor pathway. in contrast, evolutionary patterns of the insulin/foxo signaling pathway across metazoan species genomes indicated that components within the middle of the pathway were rather under stronger purifying selection. moreover, within the drosophila ras signaling pathway, upstream genes have lower evolutionary rate, suggesting more constraints, than downstream genes. obviously, different patterns of selective constraints exist for different signaling pathways. this is expected, as signaling inputs and network architecture vary among the pathways. the tlr signaling pathway represents the best characterized component of the innate immune system. tlr signaling involves the dimerization of receptors in response to pathogenic microbial products. this is followed by the recruitment of various adaptor molecules. the recruitment of these adaptor molecules leads to the recruitment of downstream signaling molecules that activate transcription factors leading to the release of inflammatory cytokines and type i interferons. the tlr signaling is therefore a sequential process. thus, the effects of pathway parameters such as gene position on the evolutionary rate of genes warrant further analyses. the tlr signaling pathway has been selectively constrained over time, indicating its essential role for host survival. in humans, population genetic studies have indicated that purifying selection is stronger within tlr adaptors relative to receptors. a recent study on the tlr signaling pathway-related genes from eight vertebrate genomes showed that the selective constraints of genes was negatively correlated with gene position along the tlr signaling pathway. thus, different components of the tlr signaling pathway appear to differ in their evolutionary rates. moreover, taxon specific differences among pathway components have been demonstrated for the tlr pathway. new insights relating to the influence of network parameters on the evolutionary rate of genes within the tlr signaling pathway can be gained through studies on closely related species not yet investigated in detail. members of the family suidae have successfully inhabited broad geographic locations including eurasia, mainland southeast asia, island southeast asia and the african continent. given the numerous pathogenic challenges across these diverse environments, the tlr signaling pathway may have been critical in the survival of suids. suid species are increasingly becoming important biomedical models given their anatomical, physiological and immunogenetics similarities with humans. comparison of protein sequences of the domestic pig, a member of the suidae, to human orthologues revealed positions where the domestic pig protein has the same amino acid that is implicated in a human disease. thus, findings regarding suid tlr signaling genes evolutionary rates can be relevant to their hominid orthologues. the availability of a complete porcine genome and resequencing libraries of species within the family suidae provide an opportunity to investigate the evolution of the tlr signaling pathway genes within the context of gene position and other network parameters. the suid species have evolved over a relative short time span of 1– million years. the low divergence among these species indicate that synonymous substitutions are not likely to be saturated, making the estimation and comparison of selective pressure among genes less prone to bias. this low divergence among the species also ensures reliability in aligning tlr sequences for subsequent analysis. in the current study, our goal was to investigate the evolutionary constraints of genes within the tlr signaling pathway of the family suidae. we sought to investigate the extent to which patterns of protein evolution observed in other organisms extend to members of the family suidae. we hypothesized that there is a polarity of purifying selection along the tlr signaling pathway within the family suidae. we tested our hypothesis by answering the following questions 1) is there a relationship between the strength of purifying selection and gene position in the tlr signaling pathway and; 2) are there any network parameters that contribute to the polarity in the strength of purifying selection? methods ethics statement blood samples were obtained from captive-held suid species by qualified individuals at each institution, according to their respective animal welfare oversight governance and institutional permits. all animal procedures were conducted in accordance with the guidelines approved by the american society of mammalogists, and approved by the university of illinois’ institutional animal care and use committee. sources of animals are provided in additional file 1: table s dna extraction and sequencing dna was extracted from whole blood using the qiaamp dna blood spin kit and quantity and quality parameters were performed using the qubit fluorometer and visualized using a 1 % agarose gel. library construction and re-sequencing of individual members of the family suidae were done with 1–3 μg of genomic dna according to the illumina library prepping protocols. the library insert size was 300–500 bp and sequencing was performed using a paired-end sequencing kit. all dna samples were sequenced to approximately 8x depth. quality trimmed reads were aligned to the sus scrofa reference genome build using the unique alignment option of mosaik aligner. the aligned reads from each of the animals together with the sus scrofa reference genome were stored as bam files for each individual animal. ortholog identification genes involved in the tlr signaling pathway were examined for this study. in addition tlr which is not found in the kegg database tlr signaling pathway, was included since it cooperates with tlr senses lipopeptides and activates the tlr signaling pathway through association with myd. the genomic coordinates of sus scrofa genes within the tlr signaling pathway were obtained from the ensembl database. based on these genomic coordinates, sequences of gene orthologs were then retrieved from aligned bam files for sus scrofa, sus verrucosus, sus celebensis, sus scebifrons, sus barbatus, babyrousa babyrussa, potamochoerus larvatus, potamochoerus porcus and phacochoerus africanus to identify gene orthologs. the resulting sequences for each species were then screened against the genome of sus scrofa using blast to ensure similarity with the sus scrofa genes. to obtain the coding sequence of genes for each species, exonic regions were retrieved based on genomic coordinates of exons from the sus scrofa gene transcripts and concatenated. if a gene was found to have more than one transcript, the longest transcript was chosen for analysis. coding sequences were aligned using clustalw. aligned sequences are provided in additional file fig. the tlr signaling pathway genes. redrawn following, with modification to include tlr the direction of signal transduction is indicated by the arrows. tlrs and are included in the first pathway position category. the numbers on the left side represent the position of the pathway genes. genes used in this study are shown in bold some of the coding sequences had missing nucleotides. the following criteria were used to discard such sequences: sequences with more than 50 % of missing amino acid residues; and aligned sequences absent in any one of the species involved in this study. consequently, the final data set was composed of genes. impact of natural selection the impact of natural selection on genes within the tlr signaling pathway was determined by estimating nonsynonymous substitution rate, synonymous substitution rate and their ratio using the m model implemented in the program codeml from the paml package version. to determine whether specific codon positions were under positive selection within gene orthologs, codeml site model m1a, a nearly neutral evolution model where sites are assumed to be evolving under either purifying selection or neutral evolution was compared to model m2a that allows positive selection among sites. m which allows sites to evolve under either purifying selection or neutrally, was compared to model m which allows for positively selected sites. models were compared using a likelihood ratio test. the f3x model of codon frequencies was used for the analyses. models were run in duplicates with ω of and to increase the probability of convergence of model parameters. multiple testing for positive selection on genes was corrected for by conducting a false discovery rate test at a q value of. the phylogenetic relationship among the species of suidae inferred from near complete genome sequences of each species was used for the analysis involving the codeml program.fig. phylogenetic relationships among species of the family suidae obtained from near complete genome data of each species. the posterior probability at each node is 1 network level analysis to determine whether the evolution of molecules within each gene was affected by network structure, network parameters were computed for each gene and were correlated with each evolutionary parameters estimated by model m the number of protein-protein interactions for proteins encoded by each gene was determined from the sus scrofa interaction network in the string database. string is a database of known and predicted interactions which include direct and indirect associations from various sources including high throughput experiments and genomic contexts. the magnitude for codon bias for individual genes was measured by the effective number of codons. the codon usage bias of each orthologous gene group was measured as the mean of enc of each species gene. enc values of each gene were obtained using dnasp software. bivariate correlation analysis between network parameters ) and length of 3′utr ) and evolutionary parameters were conducted using spearman’s rank correlations. a false discovery rate test was performed to correct for multiple testing for correlations controlling for q value at. bivariate correlations and corrections for multiple testing were performed using sas software. figures illustrating bivariate correlations were obtained with the package r. multivariate analyses the two multivariate analyses methods were used to determine whether the observed bivariate correlations were due to direct or indirect influences. partial correlation analysis measures the strength of relationship between two variables, while holding one or more variables constant. path analysis estimates direct and indirect relationships under a user defined causal model. for the path analysis, pathway position, ppi, protein length, enc and l3′utr were considered as exogenous variables whereas ω and dn were considered as endogenous variables. prior to performing the path analysis, data were log-transformed to improve normality. all statistics related to path analysis were calculated by the lavaan package in r.  orthologs of sus scrofa tlr signaling pathway genes were identified for members of the family suidae. a total of sequences of the tlr signaling pathway genes ranging from to nucleotides were therefore used in analysis. the orthologous sequences were aligned for a series of analyses including test for selection pressures acting on genes, bivariate correlations and multivariate analysis. analysis of protein sequence evolution in order to estimate the selective pressure acting on genes within the suid tlr signaling pathway, the m model, which provides a single estimate of ω across all codons and lineages, was utilized. the ω values ranged from to. the mean ω was. to test for positive selection, the models m1a vs m2a; m vs m were used. the two tests were in agreement for all genes and demonstrated that a proportion of sites within tlr tlr tlr and irak were under positive selection. however, tlr was the only gene that had sites under positive selection after correcting for multiple testing at q = . thus, results here indicated that genes within the tlr signaling pathway have evolved under strong functional constraint. relationship between evolutionary rates and pathway position the relationship between the evolutionary parameter ω of genes and their position within the tlr signaling pathway was determined using a spearman’s rank correlation test between the two variables. the ω values were negatively correlated with pathway position indicating that downstream genes are under stronger purifying selection than upstream genes. the robustness of the correlation between ω and pathway position was tested by removing the highest ω value and re-estimating the correlation. this was done to ensure that results were not influenced by high ω values. the correlation remained significant. to test for the evolutionary parameter accounting for the ω differences among genes, the correlation between pathway position and dn, pathway position and ds was also determined. the dn was negatively correlated with the pathway position. there was no significant correlation between ds and pathway position. these results indicated that the decrease in ω along the tlr signaling pathway is attributable to a decrease in dn.fig. ω and dn versus pathway position, number of protein-protein interactions and protein length. all relationships are significant. continuous lines represent regression lines relationship between evolutionary rates and network parameters network parameters other than pathway position may influence evolutionary parameters and account for inferred polarity in ω and dn. in order to determine such relationships, spearman’s rank correlation test was performed between ppi of proteins encoded by genes within the tlr signaling pathway, protein length, codon bias measured as effective number of codons, l3′utr region and the evolutionary parameters. the ω values were negatively correlated with ppi and positively correlated with protein length. similarly, dn values were negatively correlated with ppi and positively correlated with protein length indicating that apart from pathway position, the number of protein-protein interactions a given protein is involved in and its protein length influence it evolution. in addition, ds was also negatively correlated with enc indicating stronger selection based on codon usage in genes with high codon bias than genes with low codon bias. since pathway position, ppi and protein length are intercorrelated and are each correlated with evolutionary parameters, observed associations of these network parameters with the evolutionary parameters could be indirect. thus, to distinguish between direct and indirect effects, multivariate analysis were utilized. partial correlation analysis revealed that when controlling for ppi, the correlation between ω and pathway position remained significant. when controlling for protein length, the correlation between ω and pathway position remained significant. thus, gene position within the tlr signaling pathway directly influenced the evolution of genes in the tlr signaling pathway. the correlation between ω and ppi while controlling for pathway position remained significant. it also remained significant when controlling for protein length. this result indicated that ppi is an important factor affecting the evolution of genes within the tlr signaling pathway. the correlation between ω and protein length was not significant when pathway position was held constant indicating that the correlation between ω and protein length was rather mediated by pathway position. nonsynonymous substitution rates were not significantly correlated with pathway position when controlling for protein length. in contrast, the correlation between dn and ppi remained significant when controlling for protein length. if pathway position was held constant, the correlation between dn and ppi remained significant. correlation between dn and protein length was still significant after controlling for pathway position and remained significant if controlling for ppi. results for the relationship between dn values and network parameters indicated that the number of protein-protein interactions and protein length have direct effects on the rate of dn evolution. the relationship between evolutionary parameters and the network parameters were further analyzed using path analysis. as depicted in fig.  path analysis indicated that dn values were affected by gene position within the tlr signaling pathway. pathway position had the largest direct effect on dn. the dn values were negatively associated with the number of protein-protein interactions and positively associated with protein length. ω values were only associated with enc. thus pathway position, the number of protein interactions and protein length were factors that influenced tlr signaling pathway dn substitution rates in the suidae.fig. the causal model used to analyze the relationship among evolutionary and network parameters. the network parameter pathway position, number of protein-protein interactions, protein length, codon bias ) and length of the 3′utr region are considered as exogenous variables. continuous and dashed lines represent significant and nonsignificant relationships, respectively. single-headed arrows indicate causal relationship between variables. double-headed arrows indicate correlations between exogenous variables. numbers on the arrows represent the standardized regression weights discussion our results support the hypothesis that there is polarity of purifying selection along the tlr signaling pathway within the family suidae with purifying selective pressure increasing along the pathway. protein-protein interactions and protein length accounted for this polarity in the strength of purifying selection, where number of interactors of tlr signaling protein molecules increases and protein length decreasing along the tlr signaling pathway. the evolutionary parameter influenced by these network parameters was dn, which is actually the metric of selective pressure. fewer amino acid substitutions were tolerated along the signaling pathway reflective of stronger purifying selection. in a previous study on the evolution of tlr signaling pathway-related genes from eight vertebrate genomes, selective constraint was negatively correlated with gene position in the signaling pathway and dn was also influenced by position and protein length in a similar manner as in the present study. these similarities indicate that some patterns revealed across the eight vertebrate genomes are also seen in the closely related suid species involved in this study. several other studies conducted on various biological pathways have revealed a polarity in the strength of purifying selection. upstream genes revealed higher constrains than downstream genes in the carotenoid biosynthesis pathway, in the plant terpenoid biosynthesis pathway and in the dopamine catabolic pathway across mammals. furthermore, in the drosophila ras signaling pathway, upstream genes evolve under higher levels of selective constrains relative to their downstream counterparts the relatively slow rate of evolution of upstream genes reflects these genes being required for a wider range of end products and therefore being more pleiotropic and in the case of biosynthesis pathways, upstream genes exert a higher control over the flux of metabolites throughout the pathway. a negative correlation between strength of purifying selection and gene position within pathway, as inferred in this study, have also been reported for the drosophila toll and imd signaling pathways, caenorhabditis insulin/tor signal transduction pathway and the yeast hog pathway. this pattern of selective constraint might be attributed to purifying selection acting to maintain the function of downstream signal transduction elements and a concentration of adaptive changes in the upstream genes due to their interaction with external environment. in this study, tlr was the only gene inferred to be under positive selection after fdr correction. therefore, adaptive changes were not prominent within upstream genes of the tlr signaling pathway to result in the inferred pattern of selective constraints. with the inferred significantly negative correlation between pathway position and the evolutionary parameters ω and dn, a negative correlation between pathway position and number of protein-protein interactions for a protein, as indicated in this study, implies a role for the number of protein-protein interactions in the polarity of selective constraints. this negative relationship is indicative of increasing selective constraints for proteins along the tlr signaling pathway, as a result of increase in the number of interactors along the signaling pathway. evidence for higher connectivity of downstream genes has also been demonstrated for the human signal transduction network. within metabolic pathways, connectivity plays a major role in constraining evolutionary rates with proteins interacting with more proteins being subject to stronger selective constraints. the increasing number of interactors along the tlr signaling pathway indicates the essential role of protein molecules downstream of tlrs in propagating and amplifying signals initiated by pathogen binding of tlrs and the role of downstream protein molecules in other functions apart from tlr signaling. this places stronger selective constraint on the intracellular signaling protein molecules. the finding of an effect of connectivity on dn in this study is in contrast to the study of song et al. that did not show a significant influence of connectivity on dn in a multivariate analysis of factors influencing evolution of genes in the tlr signaling pathway across vertebrate genomes. a possible explanation for this difference may be that the use of suid species that are more closely related in this study than species involved in the study of song et al. led to greater precision in the estimation of dn and thus the relationship between connectivity and dn was more evident. the finding of a significant negative correlation between protein length and pathway position in this study corroborates studies in the insulin/tor signaling of drosophila and indicated that downstream genes encode shorter and more actively translated proteins. given that partial correlation and pathway analysis were both in agreement as to a positive relationship between protein length and dn, protein length appears to be a factor responsible for the polarity in dn along tlr signaling pathway. thus, along the tlr signaling pathway, the decrease in protein length is associated with decreases in amino acid substitutions. a similar pattern has been reported for the woody perennial plant populus tremula, where protein length is the main factor affecting selective constraints, with purifying selection weaker in genes with longer coding regions. a possible explanation for the relationship between purifying selection and protein length is that selection at more than one site should cause an overall reduction in the effectiveness of selection. in that case, for longer proteins, which may have many sites under selection simultaneously, there will be a reduced efficiency in natural selection. though path analysis indicated a negative association between ω and enc, pathway position and enc were not significantly correlated, therefore enc cannot explain correlation between selective constraints and enc. in this study, length of the gene 3utr region did not have a relationship with any evolutionary parameter in disagreement with studies in the toll/imd pathway in the drosophila toll and imd signaling pathways and human and mouse mirna target prediction data. genes within the signaling pathway had an average ω value <  indicating that the tlr signaling pathway is selectively constrained across members of the suidae. this is in keeping with the essential role of the pathway in innate immunity and host survival. selectively constrained regions within the genome are likely to be functionally important. the strongest selective constraints inferred for the map kinases with dn substitution rates of zero for these genes suggest their involvement in critical roles during tlr signaling. the map kinases cascade components play important role in the production of proinflammatory mediators. recently, it has been demonstrated that map kinases play a role in agonist dependent regulation of cognate tlr mrna levels. map kinases are required for the regulation of cellular development and differentiation processes. the conservation of the map kinase genes in this study could therefore be attributed to their involvement in many processes. the stronger purifying selection detected in downstream genes of tlrs relative to upstream genes suggests that there is a greater need to protect the integrity of proteins downstream. thus, genes downstream might be essential for survival of suid species involved in this study. tlr signaling evolution is partly dictated by the need to conserve the capability to recognize specific pathogen signature. for the homogenous group of suid species involved in this study, the evolutionary dynamics of the tlr signaling pathway as revealed in this study may be reflective of the range of related pathogens that adversely affected them during their evolution. our study has provided insight into factors that have shaped the present day tlr signaling pathway genes of the suid species, a convenient number of which are still present in the world. as the wild suid species are among animals at the greatest risk of declines due to cross species transmission of disease with their domesticated counterparts, evolutionary analysis of the tlr signaling pathway as presented here can aid in identifying genes that are essential in mediating the immune response to disease causing bacteria and viruses. for example, a protein that has a higher connectivity and selectively constrained within the tlr signaling pathway would be expected to be of functional relevance in terms of disease prevention.  the findings here provide evidence of the role of network topology in the polarity of purifying selection along the tlr signaling pathway of the suidae. our work indicates that the suidae are an attractive system of closely related but well characterized species to study the evolutionary dynamics of tlr signaling genes. furthermore, our work opens other avenues of research including the use of the suidae in understanding patterns of molecular evolution in other important pathways such as metabolic and gene-regulation pathways. availability of supporting data all supporting data are included as additional files in the form of additional files and additional files additional file 1: table s sample information. additional file 2: table s genes of the sus scrofa tlr signaling pathway used in querying genomes of other species of the suidae. additional file 3: multiple sequence alignments of the tlr signaling pathway orthologs across the suid species. asteriks indicate same nucleotide and “-“ represents missing nucleotide. additional file 4: table s summary statistics for genes. additional file 5: table s results of codon based test for positive selection. additional file 6: table s bivariate correlations among variables. competing interests the authors declare that they have no competing interests. authors’ contributions conceived and designed the experiments: kado mamg lbs. performed the experiments: kado. analyzed the data: kado hjm alr mamg lbs. contributed reagents/materials/analysis tools: hjm alr mamg lbs. wrote the paper: kado lbs. all authors read and approved the final manuscript. 
selective constraints could vary depending on the gene’s position within the pathway as tlr signaling is a sequential process and that genes downstream of the tlrs may be more selectively constrained to ensure efficient immune responses given the important role of downstream genes in the signaling process. a significant relationship was inferred for the network parameters gene position, number of protein-protein interactions, protein length and the evolutionary parameter dn such that downstream genes had lower nonsynonymous substitution rates, more interactors and shorter protein length than upstream genes. gene position was significantly correlated with the number of protein-protein interactions and protein length. in particular, downstream genes evolve more slowly as a result of being highly connected and having shorter protein lengths. results indicate that the level of selective constraints on genes within the tlr signaling pathway of the suidae is dependent on the gene’s position and network parameters. suidaeselective constraintsevolutionsignalingupstream and downstream genesus department of agriculture agriculture research service grant58-5438-2-307schook lawrence b. issue-copyright-statement© the author 2016 the members of the suidae examined included the european sus scrofa, asian sus scrofa, sus verrucosus, sus celebensis, sus scebifrons, sus barbatus, babyrousa babyrussa, potamochoerus larvatus, potamochoerus porcus and phacochoerus africanus. these findings highlight the critical role of gene network parameters in gene evolution. highly conserved proteins, indicative of their critical roles in host survival, characterize this pathway. the evolutionary parameter ω had an overall mean of across genes, indicating high functional conservation within the tlr signaling pathway. thus, we investigated whether gene position influenced protein evolution in the tlr signaling pathway of the suidae. a total of tlr signaling pathway genes in the suidae were retrieved from resequencing data. thus, the polarity in the selective constraint along the tlr signaling pathway was due to the number of molecules a protein interacted with and the protein’s length.
 the toll-like receptor signaling pathway constitutes an essential component of the innate immune system. highly conserved proteins, indicative of their critical roles in host survival, characterize this pathway. selective constraints could vary depending on the gene’s position within the pathway as tlr signaling is a sequential process and that genes downstream of the tlrs may be more selectively constrained to ensure efficient immune responses given the important role of downstream genes in the signaling process. thus, we investigated whether gene position influenced protein evolution in the tlr signaling pathway of the suidae. the members of the suidae examined included the european sus scrofa, asian sus scrofa, sus verrucosus, sus celebensis, sus scebifrons, sus barbatus, babyrousa babyrussa, potamochoerus larvatus, potamochoerus porcus and phacochoerus africanus.  a total of tlr signaling pathway genes in the suidae were retrieved from resequencing data. the evolutionary parameter ω had an overall mean of across genes, indicating high functional conservation within the tlr signaling pathway. a significant relationship was inferred for the network parameters gene position, number of protein-protein interactions, protein length and the evolutionary parameter dn such that downstream genes had lower nonsynonymous substitution rates, more interactors and shorter protein length than upstream genes. gene position was significantly correlated with the number of protein-protein interactions and protein length. thus, the polarity in the selective constraint along the tlr signaling pathway was due to the number of molecules a protein interacted with and the protein’s length.  results indicate that the level of selective constraints on genes within the tlr signaling pathway of the suidae is dependent on the gene’s position and network parameters. in particular, downstream genes evolve more slowly as a result of being highly connected and having shorter protein lengths. these findings highlight the critical role of gene network parameters in gene evolution. electronic supplementary material the online version of this article contains supplementary material, which is available to authorized users. keywords suidaeselective constraintsevolutionsignalingupstream and downstream genesus department of agriculture agriculture research service grant58-5438-2-307schook lawrence b. issue-copyright-statement© the author 2016 