Protein sites with more coevolutionary connections tend to evolve slower, while more variable protein families acquire higher coevolutionary connections

Background: Amino acid exchanges within proteins sometimes compensate for one another and could therefore be co-evolved. It is essential to investigate the intricate relationship between the extent of coevolution and the evolutionary variability exerted at individual protein sites, as well as the whole protein. Methods: In this study, we have used a reliable set of coevolutionary connections (sites within 10Å spatial distance) and investigated their correlation with the evolutionary diversity within the respective protein sites. Results: Based on our observations, we propose an interesting hypothesis that higher numbers of coevolutionary connections are associated with lesser evolutionary variable protein sites, while higher numbers of the coevolutionary connections can be observed for a protein family that has higher evolutionary variability. Our findings also indicate that highly coevolved sites located in a solvent accessible state tend to be less evolutionary variable. This relationship reverts at the whole protein level where cytoplasmic and extracellular proteins show moderately higher anti-correlation between the number of coevolutionary connections and the average evolutionary conservation of the whole protein. Conclusions: Observations and hypothesis presented in this study provide intriguing insights towards understanding the critical relationship between coevolutionary and evolutionary changes observed within proteins. Our observations encourage further investigation to find out the reasons behind subtle variations in the relationship between coevolutionary connectivity and evolutionary diversity for proteins located at various cellular localizations and/or involved in different molecular-biological functions.

: In this study, we have used a reliable set of coevolutionary Methods connections (sites within 10Å spatial distance) and investigated their correlation with the evolutionary diversity within the respective protein sites.
: Based on our observations, we propose an interesting hypothesis that Results higher numbers of coevolutionary connections are associated with lesser evolutionary variable protein sites, while higher numbers of the coevolutionary connections can be observed for a protein family that has higher evolutionary variability. Our findings also indicate that highly coevolved sites located in a solvent accessible state tend to be less evolutionary variable. This relationship reverts at the whole protein level where cytoplasmic and extracellular proteins show moderately higher anti-correlation between the number of coevolutionary connections and the average evolutionary conservation of the whole protein.
: Observations and hypothesis presented in this study provide Conclusions intriguing insights towards understanding the critical relationship between coevolutionary and evolutionary changes observed within proteins. Our observations encourage further investigation to find out the reasons behind subtle variations in the relationship between coevolutionary connectivity and evolutionary diversity for proteins located at various cellular localizations and/or involved in different molecular-biological functions.

Referee Status:
Invited Referees According to the neutral theory of evolution, the functionality of a protein with a disadvantageous mutation can be restored by another mutation that compensates for the first to sustain the function 1 . Such compensating mutations, together with other factors arising due to common functional, structural and folding constraints, lead to correlations between different positions in a protein or protein family. Coordinated changes of amino acid residues are typically acquired by examining covariation between two aligned positions. A large number of computational methods have been proposed [2][3][4][5][6][7][8][9][10][11] to quantify the covariation between two protein sites in a given multiple sequence alignment (MSA). Most methods are based on variation of mutual information 12-17 , maximum likelihood approximations 18 , Bayesian probabilities 19 , and phylogenetic approaches 20,21 . Newer methods successfully implement direct coupling analysis 22 , Protein Sparse Inverse COVariance: PSICOV 23 and Matrix Match Maker 24 algorithms to identify coevolving sites. These previous studies demonstrate that sequence covariation is powerful in detecting protein-protein interactions, ligand receptor binding, and the folding structure of a protein. In addition to direct physical interactions, distantly located coevolving amino acid residues are reported to be energetically coupled 25 or subject to similar functional constraints 26 . Compensated amino acid substitutions have been described in previous works in terms of their locations in structure and their physico-chemical properties 3,20,21 . Coevolutionary signals coming from residue charge compensating mutations have been found to be stronger compared to size compensating mutations 3,21,27 . Despite the fact that coevolution has been found to be rather weak in many cases, correlated mutations have had comparative success in predicting protein secondary and tertiary structures, and in some cases protein interaction partners [28][29][30] .
Coevolution is difficult to detect due to various reasons, such as the variable nature of compensatory mutations, the strong dependence of covariations on evolutionary distances, and the number of sequences in the alignment. Hence, it is crucial to understand how coevolutionary processes are related to evolutionary diversity within protein families. Despite significant efforts in this field, the relationship between evolutionary conservation and the extent of coevolution is not well understood. For example, it is not clear whether families with higher evolutionary diversity would exhibit more coevolutionary connections or not. Similarly, at the residue level, this relationship needs to be thoroughly examined. An earlier study by Fodor and Aldrich 31 observed a lack of agreement between correlated mutation methods, and the resultant differences might have been caused by differing sensitivities to background conservation. In a previous study, it was also indicated that residues, which form many coevolutionary connections with other residues, are more evolutionary conserved and are involved in specific functionally important interactions and conformational changes 32 .
A complete understanding of protein evolution and coevolution will require a large scale analysis of important factors that determine the selective forces acting on different residues of a protein to be coevolved. Here, we present a study that undertook a detailed analysis to investigate the relationship between evolutionary conservation and the extent of coevolution within a protein. This relationship could be dependent on the reliability of the predicted coevolved sites as there are no direct ways to validate the coevolutionary connectivity. Therefore, it is a good idea to use multiple coevolution extracting algorithms and filter out a reliable set within protein sites. Similarly, spatial proximity between the coevolved sites might provide additional reliability about the predicted coevolved sites 33 . We examined the evolutionary conservation using the popular AL2CO  Our study suggests the hypothesis that a higher number of coevolutionary connection is likely to be observed for a particular site that is less evolutionary variable, while a higher number of coevolutionary connections can be observed for a protein family that has higher evolutionary variability. We found that the sites with a higher number of coevolutionary connections have a much higher tendency to be conserved compared to the sites with a smaller number of connections. These sites might act as 'hub points', and therefore changes in these sites would affect many other connected sites. We further investigated the impact of important structural properties, like secondary structures, solvent accessibilities and hydrogen bonding of the coevolved sites, to understand the reasons behind the observed correlation between coevolution and evolutionary diversity. Our findings indicate that coevolved sites are generally preferred at a solvent accessible/hydrogen bonded/helical state compared to a solvent buried/non-hydrogen bonded/β strand state. However, discernable differences in evolutionary conservation between the higher and lesser coevolved sites were observed only for sites located at solvent accessible states compared to buried states. We also examined whether the observed negative (anti) correlation between coevolution and evolutionary conservation for a protein family can be under the

Amendments from Version 1
Version 2 of this manuscript includes the following changes: 1: Change of a line in the Abstract from "Correlated mutation or coevolution of positions in a protein is tightly linked with the protein's respective evolutionary rate", to "Amino acid exchanges within proteins sometimes compensate for one another and could therefore be co-evolved." 2: Included detail of coevolutionary programs' range of values.
3: Incorporated result of the MISTIC server to compare with our results.

4:
We have included Figure S7 as a MISTIC server Circos representation for CD01291 and CD00164 families. Table S2 which contains the MISTIC server's web link results for a few randomly selected CDD families.

5: We have included
6: Included mention of project GAP 362 in the Grant information.

REVISED
influence of its cellular localization or the type of functions with which it is involved. Coevolution analysis for the whole protein suggests that the cytoplasmic and extracellular proteins possess moderately stronger negative (anti) correlation between the number of coevolutionary connections and their average evolutionary conservation.

Dataset
We collected 753 protein domain alignments from the Conserved Domain Database (CDD; https://www.ncbi.nlm.nih.gov/Structure/ cdd/cdd.shtml 35 ) version 2.13, for which at least one 3D structure entry and more than 50 protein sequences are available. An alignment length threshold (>=100) was also applied to exclude smaller proteins. A complete list of protein families is provided in Table S1.
Identification of coevolved sites Mutual information 6 (Suppl. Mat. Ref.) is widely used measure to estimate the covariation between sites in protein families. In this analysis, we used a mutual information based method to estimate coevolutionary connection between two sites of a protein family. This method (MIp) is based on information theory that accurately estimates the expected levels of background coming from random and phylogenetic signals. Removal of the phylogenetic and random background allows identifying substantially more coevolving positions in protein families. Altogether we identified 19,736 (out of total 36,616) coevolved site pairs located within 10Å spatial distance from the 753 family alignments, with a MIp Z-score cutoff of 4.0 or higher.
McBASC 27 (http://fodorwebsite.appspot.com//covariance1_1.zip) was used to calculate the simple inter-position coevolution for the 753 protein family alignments. McBASC provides high score for non-conserved and co-varying positions from a multiple sequence alignment. The calculation of McBASC was performed as described in Fodor and Aldrich 2004, using the software provided by the authors (http://www.afodor.net/). McBASC does not use any structural or phylogenetic information in the calculation of coevolution. We identified 35,514 (out of total 95,866) coevolved site pairs located within 10Å spatial distance from the 753 family alignments with McBASC Z-score cutoff of 4.0 or higher.
DCA 22 (Direct Coupling Analysis) aims at predicting coevolving residues based on the maximum entropy principle. DCA is also used in predicting inter and intra domain contacts. This method is used in separating direct and indirect correlation between residues. DCA analysis was implemented with MATLAB code kindly provided to us by Domenico L. Gatti (Supplementary File 1). We identified 50,217 (out of total 1,61,332) coevolved site pairs located within 10Å spatial distance from the 753 family alignments with DCA Z-score cutoff of 4.0 or higher.
PSICOV 23 (Protein Sparse Inverse COVariance) method is developed with the specific goal of separating direct from indirect coupling between residues. PSICOV takes into account the global correlations between pairs. Modified MATLAB code (without the default minimum requirement of 500 sequences), which was kindly provided to us by Domenico L. Gatti (Supplementary File 1), was used in this study. We identified 56,879 (out of total 162,336) coevolved site pairs located within 10Å spatial distance from the 753 family alignments with PSICOV Z-score cutoff of 4.0 or higher.

Random selection of non-coevolved sites and pairs
Site pairs other than those involved in coevolutionary connections were considered as non-coevolutionary sites. We randomly selected non-coevolved sites from each protein family (Supplementary File 2). For each randomly selected non-coevolved site (i), neighboring non-coevolved sites were selected based on the structural distance (<10Å) and sequence distance filters (>i±6 positions). Similar numbers of non-coevolved site pairs were selected randomly 10 times. We performed similar correlation analysis between the numbers of spatial neighbors and evolutionary conservation of non-coevolved sites.

Calculation of amino acid conservation
Analysis of positional conservation in a sequence alignment can aid in the detection of functionally and/or structurally important residues. The AL2CO 34 program performs conservation analysis in a comprehensive and systematic way. It was used to calculate the conservation index for each position for a given multiple sequence alignment. Twelve different strategies of conservation index calculation have been implemented in the AL2CO program (http:// prodata.swmed.edu/al2co/al2co.php). For this analysis, we used independent count (sequence weighting scheme) and matrix based sum-of-pair 36 (conservation calculation method) measure scoring scheme to calculate evolutionary conservation of each coevolved sites or column in the alignment. A higher AL2CO score indicates higher conservation index.

Calculation of spatial distances and structural properties
Representative three-dimensional (3D) structures were collected for each family from the Protein Data Bank (PDB; http://www. rcsb.org/pdb/home/home.do) 37 . Spatial distances were calculated using atom coordinates supplied in the individual PDB file. Structural properties, such as solvent accessibility, secondary structures, and hydrogen bonds, were computed from the protein structure using the JOY package 38 (http://mizuguchilab.org/joy/) Solvent accessibility was measured using the PSA program from the JOY package, and residues that had an accessible surface area <7% were treated as solvent buried or inaccessible. Similarly, secondary structures (helix, strand and coil) and hydrogen bonding patterns were estimated using the SSTRUC and HBOND programs from the JOY package 39 , respectively.

Collection of Gene Ontology information
The Gene Ontology (http://www.geneontology.org/) 40 covers three classes/domains: cellular localization, molecular function and biological process. Functional information of each CDD family was collected from Gene Ontology database using the UNIPROT 39 ID of the representative protein structure as a query. We mapped 517, 720, and 634 protein domain families into cellular localization, molecular function and biological process, respectively.
Mapping conservation and coevolutionary connection onto 3D structure Mapping of evolutionary conservation and coevolutionary information onto the 3D structure was done using in-house Perl scripts (Supplementary File 3). B-Factor column in PDB file was substituted with evolutionary conservation score and colored according to B-Factor ranging blue (low conservation) to red (high conservation). Lines connecting C-alpha atoms of residues represent coevolutionary connection between those residues.

Results and discussion
Coevolution versus evolutionary diversity at the site level Coevolutionary connections between protein sites were identified from multiple sequence alignments of 753 protein domain families by algorithms employing differing approaches, such as mutual information (MIp program 6 ), McLachlan amino acid similarity matrix based techniques (McBASC program 27 ), Direct Coupling Analysis (DCA program 22 ), and Protein Sparse Inverse COVariance method (PSICOV program 23 ). Minimal overlaps were observed for coevolved sites predicted by these programs ( Figure S1), supporting previous interpretations that differences in the preferred level of background conservation may exist within each program to identify coevolved residue pairs 6 .
The pattern of evolutionary diversity within the coevolved sites was examined using evolutionary conservation scoring approaches (e.g., AL2CO). Figure 1 plots the average conservation scores of sites having higher or lower coevolutionary connections (A: MIp; B: McBASC; C: DCA; D: PSICOV programs, respectively). Figure 1 suggests that highly coevolved sites possess higher average AL2CO scores, depicting higher evolutionary conservation. Coevolutionary connections, even though selected based on a strong statistically significant threshold (Z-score >4), might contain background noise resulting in an unreliable relationship between coevolution and evolutionary conservation. To disprove this, we performed similar analysis using noncoevolutionary random sites and found that there is a smaller correlation between non-coevolved sites having higher or lower structural distance based neighbors (<10Å) and their evolutionary conservation ( Figure S2).
Observation of strong positive correlation between coevolutionary connections and evolutionary conservation within the coevolved sites selected based on structural proximity suggests that highly coevolved protein sites tend to evolve slower.
Influence of structural environment. The structural environment of a protein site is a critical factor that can influence its evolutionary diversity pattern 41,42 . To understand the reasons behind the observed phenomenon where higher coevolutionary connections are found for sites that are less diversified, we investigated the roles of structural environments, such as solvent accessibility state, and secondary structural content of the coevolved sites.
We have observed more coevolutionary connections for sites that are solvent accessible compared to that observed within buried sites ( Figure 2A). Interestingly, solvent accessible sites that possess lower numbers (<3) of coevolutionary connections (LCC) are consistently less conserved compared to the sites that have relatively higher number (>3) of coevolutionary connections (HCC) (Figure 2A). Although the similar trend is also observed within the solvent buried sites, the differences of conversation indices between the HCC and LCC are more prominent within solvent accessible state compared to that observed at buried state ( Figure 2B).
Higher abundance of coevolutionary connections is also observed for sites that are involved in hydrogen bonding compared to those are not involved in hydrogen bonding. However, no discernable differences in evolutionary conservation were observed between the higher and lesser coevolved sites involved in hydrogen bonding compared to those that do not have hydrogen bonding ( Figure S3).
Slightly higher abundance of coevolutionary connections was observed for sites that were located in helix compared to those forming strands. No discernable differences in evolutionary conservation were observed between the higher and lesser coevolved sites located at helical environments compared to those that form strands ( Figure S4).

Influence of functional involvement.
We also investigated the relationship between coevolutionary connection and evolutionary conservation for protein sites with respect to their functional involvement. However, functional sites (e.g., active sites, protein or ligand binding sites) do not show significantly higher positive correlation between coevolutionary connection and evolutionary conservation, and no discernable differences were observed among the correlation coefficients between coevolutionary connection and evolutionary conservation observed for various types of functional sites (data not visualised).
Coevolution versus evolutionary diversity at the protein/ family level It is important to know how the evolutionary conservation profile of the whole protein or family influences the coevolutionary connections within its sites. Figure 3 and Figure S5 plot the average conservation scores of protein families (considering all gapless columns of the family alignment) with respect to the total number of coevolved sites observed within those families. Our results suggest strong negative correlation between the number of coevolved sites found within a protein family and its average conservation score. This finding indicates that, in general, more conserved proteins/families tend to possess lower coevolutionary connections, whereas proteins/families with less stringent evolutionary pressure might engineer more intra-coevolutionary connections.
We further investigated the influence of cellular localization and biological-molecular functions of the proteins that displayed correlation between the coevolutionary connections and evolutionary conservation. We categorized the representative proteins from 517, 720, 634 families into cellular localization, molecular function and biological processes, respectively, using their Gene Ontology annotations. For example, 54%, 15% and 12% of the 517 families, having at least one pair of coevolved sites, reside within cytoplasm, nucleus and membrane, respectively ( Figure S6). Similarly, 55%, 17% and 10% coevolved protein families are involved in catalysis (enzyme), nucleic acid and ion binding functions, respectively. Coevolved proteins were also found to be abundant in various metabolic functions ( Figure S6). Table 1 provides the R 2 and slope (m) values between the coevolutionary  connection and evolutionary conservation for proteins categorized in certain cellular localization. Cytoplasmic and extracellular proteins show slightly stronger anti-correlation between the number of their coevolutionary connections and evolutionary conservation. Similarly, proteins involved in catalysis and nucleic acid binding type of molecular functions show moderately stronger negative correlation, whereas proteins involved in miscellaneous metabolic processes, which mostly include generic carbohydrate and glutamine metabolisms and nitrogen fixation processes, exhibit stronger negative correlation between coevolutionary connections within the protein and its average conservation ( Table 1). Figure 4A provides an example case where coevolutionary connections are overlaid with evolutionary conservation scores onto the 3D structure of a representative protein (PDB code: 1DJ0) from the pseudouridine synthase domain family (CDD code: CD01291). 8, 30, 20 and 46 coevolutionary connections were predicted by MIp 6 , McBASC 27 , DCA 22 and PSICOV23 methods, respectively. Interestingly, in this family, the average conservation score (AL2CO score: 0.65) for all sites are quite low (as shown by color coding), despite having higher coevolutionary connections. Hence, observations in this family support the hypothesis that a higher number of coevolutionary connections can be expected for a protein family that has higher evolutionary variability or lower evolutionary conservation. Similarly, Figure 4B provides a case where coevolutionary connections are projected onto the 3D structure of a representative protein (PDB code: 1SRO) from the ribosomal protein S1 domain (CDD code: CD00164). It is evident from Figure 4B that the number of coevolutionary connections is relatively low in this family, while the overall  Files of coevolutionary sites predicted by four programs with conservation score predicted by AL2CO program with 10Å filter.

Example cases
In order to compare our findings, we have performed a similar analysis using the MISTIC 43 (Mutual Information Server To Infer Coevolution) server, taking 20 randomly selected protein families from our dataset as case study. Interestingly, the observed coevolutionary network predicted for CD01291 and CD00164 families (discussed above) are similar to our study ( Figure S7). MISTIC results also show that the CD01291 family has higher coevolutionary network connections for fewer variables sites whereas the CD00164 family has less coevolutionary connections and is overall less conserved. The MISTIC server's web link results for other protein families, are available in Table S2.

Conclusions
Over the years, it has become apparent that intra protein coevolution is an important evolutionary phenomenon to maintain proteins' functional flexibility. However, the signs of coevolution are subtle, and as a consequence, hard to detect. The majority of sites in a protein coevolve to some degree, in that they contribute more or less to structural integrity and, thus, function of the protein. However, some sites will more directly influence each other. By definition, coevolution is closely connected to the evolutionary variability of a protein. Hence, it is essential to investigate the intricate relationship between the extent of coevolution and the evolutionary variability exerted at individual protein sites, as well as the whole protein. However, it is also relevant to check the reliability of the predicted coevolved sites before deriving any hypothesis between coevolution and evolutionary conservation. Therefore, we employed multiple algorithms for the detection of coevolutionary connection and used a structural proximity based filtration system to validate the coevolutionary connections within protein sites.
In this study we have not checked/compared the difference between the two concepts of covariation and coevolution. We have used different programs (MIp, McBASC, DCA and PSICOV) which calculate covariation among protein sites in tree-independent manner. In this study, it was assumed that observed patterns of covariation are caused by molecular coevolution and they were treated synonymously. To the best of our knowledge, this is the first time where such a detailed analysis is performed to investigate any existing correlation between the coevolution and evolutionary conservation. Based on our observations, we propose an interesting hypothesis that a higher number of coevolutionary connection is associated for a protein site that is less evolutionary variable, while a higher number of the coevolutionary connections can be observed for a protein family that has higher evolutionary variability. The obvious question is why such apparently contrasting relationship exists. One probable explanation could be that these highly coevolved sites might act as 'coevolutionary hubs', and therefore changes at these sites would affect many other connected sites. On the contrary, the evolutionary selection pressure needs to be lower at the whole protein for more sites to be involved in covariation.
Probably, sites that are critical to maintain structural integrity and functional flexibility are co-varying with many other sites, but the extent of variation is limited. Hence, the critical balance between covariation and evolutionary conservation is maintained via these 'coevolutionary hub' sites. However, to be rich in a coevolutionary connection, a protein requires evolutionary flexibility so that correlated or compensatory mutations can be arranged with response to an initial change. Hence, higher coevolutionary connection is observed for families that are more evolutionary variable than others.

Dataset 1: Predicted data for coevolution and conservation.
Files of coevolutionary sites predicted by four programs with conservation score predicted by AL2CO program with 10Å filter. doi, 10.5256/f1000research.11251.d157108 44 Author contributions SC designed the work, analyzed and interpreted data; SM conceived the experiments, analyzed and interpreted data. SC and SM wrote the manuscript.

Competing interests
No competing interests were disclosed.

Grant information
The work is supported by the CSIR network project fund (BSC 0121) and GAP 362.
The funders had no role in study design, data collection and analysis, decision to publish, or preparation of the manuscript.  Supplementary File 1: Program for coevolutionary connection prediction. MATLAB code for DCA and PSICOV coevolutionary connection prediction program. This code is provided by Domenico L. Gatti with permission. Click here to access the data.

Supplementary File 2: Program to extract non-coevolved sites.
Click here to access the data.

Supplementary File 3: Program and data-files to map conservation and coevolutionary information onto PDB.
Click here to access the data.

Figure S1: Number of coevolved pairs predicted (in blue, pink, green and yellow) by different programs (MIp, MCBASC, DCA and PSICOV, respectively) and common pairs (in black) between them.
Click here to access the data. Panels provide correlation data for the non-coevolved sites (i) that are located within or equal to 10Å (and at sequence position between i and >i±6). The coefficient of determination (R 2 ) indicates how well the data points fit to the linear regression model between coevolutionary connection and evolutionary conservation. Click here to access the data.    (1) Representative proteins from 517 CDD families were assigned to cellular localization, whereas the same from 720 and 624 families could be assigned to at least one (2) molecular function or (3) biological process (3), respectively. Details can be found in Methods. Click here to access the data. The evolutionary conservation of a large number of predicted co-evolving residue pairs has been investigated for possible correlation between the extent of conservation and the strength of co-evolving residue networks. Co-evolution has been predicted by four different popular algorithms. Residues with a high networking of co-evolved residues are found to be more evolutionarily conserved. However, the same trend is not true at the entire protein domain family level and evolutionarily conserved protein families appear to exhibit less co-evolved network of residues. Likewise, solvent-accessible residues were predicted to retain more co-evolutionary connections in comparison to solvent-buried residues. These are interesting observations, but the connections between these individual observations and possible implications/applications will be good to include in the paper. Queries: The first sentence in the Abstract could be changed to "Amino acid exchanges within proteins sometimes compensate for one another and could therefore be co-evolved." since this fact of tight linking is not well-known and forms one of the questions in this study.
well-known and forms one of the questions in this study.
Page 4: It will be nice to explain how the conservation score (within the AL2CO program) is calculated.
Was there no check for consensus in predicting the co-evolved residues? For instance, to see whichever are predicted by three or more methods … It will be interesting to examine the results for subset of such highly predicted co-evolved sites.
Page 4: "Representative three-dimensional (3D) structures were collected for each family from the Protein Data Bank" -to provide details as to how they were selected?
Page 5: This statement "Observation of strong positive correlation between coevolutionary connections and evolutionary conservation within the coevolved sites selected based on structural proximity suggests that highly coevolved protein sites tend to evolve slower." seems to be apparently counter-intuitive. How can highly conserved sites be co-evolving also? Highly conserved sites usually imply high degree of identity (self-amino acid preservation). If so, how a co-evolutionary index can be set up for two spatially proximate residues which remain identical? Please explain for the benefit of the readers.

If applicable, is the statistical analysis and its interpretation appropriate? Yes
Are all the source data underlying the results available to ensure full reproducibility? Yes The evolutionary conservation of a large number of predicted co-evolving residue pairs has been investigated for possible correlation between the extent of conservation and the strength of co-evolving residue networks. Co-evolution has been predicted by four different popular algorithms. Residues with a high networking of co-evolved residues are strength of co-evolving residue networks. Co-evolution has been predicted by four different popular algorithms. Residues with a high networking of co-evolved residues are found to be more evolutionarily conserved. However, the same trend is not true at the entire protein domain family level and evolutionarily conserved protein families appear to exhibit less co-evolved network of residues. Likewise, solvent-accessible residues were predicted to retain more co-evolutionary connections in comparison to solvent-buried residues. These are interesting observations, but the connections between these individual observations and possible implications/applications will be good to include in the paper.
We thank the reviewer for the positive comments. Current study does not provide practical application in its current form, but does offer insight into the underlying properties of covariation/coevoluution methods and the relationship of these methods with evolutionary rate. However, the knowledge regarding the intricate relationship between evolutionary variability and coevolutionary connection is very important to gain insight about the dynamics and pattern of evolutionary history of protein families. The variable nature of this intricate balance is perhaps crucial in determining the overall conservation and/or flexibility of functionally important sites within certain protein families.
The first sentence in the Abstract could be changed to "Amino acid exchanges within proteins sometimes compensate for one another and could therefore be co-evolved." since this fact of tight linking is not well-known and forms one of the questions in this study.
We thank the reviewer for her thoughtful comment. We agree to change the line in the abstract.
Page 4: It will be nice to explain how the conservation score (within the AL2CO program) is calculated.
Information on conservation score calculation is provided in manuscript subsection "Calculation of amino acid conservation". The AL2CO (manuscript reference: 34) program performs conservation analysis in a comprehensive and systematic way. We used independent count (sequence weighting scheme) and matrix based sum-of-pair (conservation calculation method) measure scoring scheme of AL2CO program to calculate evolutionary conservation of each coevolved sites or column in the alignment. These scoring functions use the summation of the products of frequencies for the column for every combination of amino acid and and multiplies the products a b by the corresponding BLOSUM62 matrix amino acid substitution frequencies.
Was there no check for consensus in predicting the co-evolved residues? For instance, to see whichever are predicted by three or more methods … It will be interesting to examine the results for subset of such highly predicted co-evolved sites.
Page 4: "Representative three-dimensional (3D) structures were collected for each family from the Protein Data Bank" -to provide details as to how they were selected?
We have selected 3D structure of each family from the conserved domain database (CDD) alignment (represented as first sequence in alignment file).
Page 5: This statement "Observation of strong positive correlation between coevolutionary connections and evolutionary conservation within the coevolved sites selected based on structural proximity suggests that highly coevolved protein sites tend to evolve slower." seems to be apparently counter-intuitive. How can highly conserved sites be co-evolving also? Highly conserved sites usually imply high degree of identity (self-amino acid preservation). If so, how a co-evolutionary index can be set up for two spatially proximate residues which remain identical? Please explain for the benefit of the readers.
We thank the reviewer for these useful comments. We agree with the apparent counter intuitiveness of the sentence. However, it is the fact that we observed. The obvious question is why such apparently contrasting relationship exists. Perhaps, both coevolutionary and evolutionary changes are dynamic processes and for a given protein site at a certain point of evolutionary conservation status, highest coevolutionary connections are observed. This evolutionary conservation status of the site is perhaps selected and maintained. One probable explanation could be that these highly coevolved sites might act as 'coevolutionary hub' and therefore changes at these sites would affect many other connected sites. However, we must mention that in this study the higher conservation is with respect to other coevolving sites and not necessarily meant completely conserved sites.  (Ref 32) and shows that the number of inter-residue coevolutionary relationships can be correlated with the evolutionary conservation of a protein site and protein family. The authors applied four different algorithms to calculate the coevolutionary relationships between sites and an overall trend observed in this study is confirmed by different methods. Interestingly, the absolute scale of site conservation for sites with the same number of coevolutionary relationships can differ drastically between methods (for example McBASC and MIp on Figure 1). I wonder if the sets of pairwise correlated sites overlap between different methods. I would also suggest using MISTIC server which can provide information on conservation, coevolution and structure mapping. The relationships between coevolution and diversity of protein families is interesting and intriguing, can it be related to the quality of alignments, one of the major factors defining the accuracy of coevolutionary detection quality of alignments, one of the major factors defining the accuracy of coevolutionary detection algorithms? It is important also to discuss the difference between covariation and coevolution, the latter is not necessarily the cause, see some recent studies: PMID:25944916).

Is the work clearly and accurately presented and does it cite the current literature? Yes
Is the study design appropriate and is the work technically sound? Yes

If applicable, is the statistical analysis and its interpretation appropriate? Partly
Are all the source data underlying the results available to ensure full reproducibility? Yes

Are the conclusions drawn adequately supported by the results? Partly
No competing interests were disclosed.

Competing Interests:
I have read this submission. I believe that I have an appropriate level of expertise to confirm that it is of an acceptable scientific standard, however I have significant reservations, as outlined above.
Author Response 14 Jun 2017 , Indian Institute of Chemical Biology, India Saikat Chakrabarti This paper expands on a previous study (Ref 32) and shows that the number of inter-residue coevolutionary relationships can be correlated with the evolutionary conservation of a protein site and protein family. The authors applied four different algorithms to calculate the coevolutionary relationships between sites and an overall trend observed in this study is confirmed by different methods. Interestingly, the absolute scale of site conservation for sites with the same number of coevolutionary relationships can differ drastically between methods (for example McBASC and MIp on Figure 1).
We thank the reviewer for the positive comments. We agree with reviewer's comment, the observed scale of values obtained from multiple coevolutionary programs varies a lot. The probable reason for such observation can be the algorithm used by individual programs for calculation of covariation/coevolution. We are providing our point-by-point response in the following.

I wonder if the sets of pairwise correlated sites overlap between different methods.
Figure S1 in supplementary file 1 shows number of coevolved pairs predicted (in blue, pink, green Figure S1 in supplementary file 1 shows number of coevolved pairs predicted (in blue, pink, green and yellow) by different programs (MIp, MCBASC, DCA and PSICOV, respectively) and common pairs (in black) between them. I would also suggest using MISTIC server which can provide information on conservation, coevolution and structure mapping.
We thank the reviewer for the comment. MISTIC ( ) is mutual information server to infer coevolution an online server, hence running for large number of protein families is not feasible. However, as case studies, we have performed the analysis on MISTIC server for 20 protein families (including CD01291 and CD00164 families provided in the paper). Result of MISTIC server for CD01291 is available at and for CD00164 at http://mistic.leloir.org.ar/results.php?jobid=201705252335594338 . http://mistic.leloir.org.ar/results.php?jobid=201705269510226 Circos representation of result for both the families is as follows: CD01291 (Pseudouridine synthases family): CD00164 (Ribosomal protein S1-like RNA-binding domain): MI Circo is a sequential circular representation of the MSA and the information it contains. Coloured square boxes of the second circle indicate the MSA position conservation (highly conserved positions are in red, while less conserved ones are in blue). Lines connect pairs of positions with MI greater than 6.5 (Marino Buslje , 2009). edges et al Red represent the top 5%, ones are between 70% and 95%, and edges account for the black gray remaining 70%.
Interestingly observed coevolutionary network predicted for both the families are similar to our study. Where CD01291 family has higher coevolutionary network connections for fewer variables sites whereas CD00164 family has less coevolutionary connections and overall less conserved.
No competing interests were disclosed. Competing Interests: