we used the complete sets of annotated orthology data of 20 culicidae mosquito genomes ( supplementary table s1 , supplementary material online ) from http://cegg.unige.ch/orthodbmoz2 ( neafsey et al . we selected 2,008 single - copy genes which that contained sequences from all 20 species as our initial data set . we aligned all genes using the mafft software , version 7.182 ( katoh and toh 2008 ) based on their amino acid sequence , using the default settings ( automatic selection of the appropriate strategy , from l - ins - i , fft - ns - i , and fft - ns-2 , according to data size ) . 2006 ) to translate amino acid sequence alignments to codon sequence alignments , and the  automated  option of trimal ( capella - gutierrez et al . trimmed segments of the amino acid sequence alignments were deleted from their corresponding codon sequence alignments using custom perl scripts . following trimming , our data matrix consisted of 2,007 genes from 20 species and contained no missing data . to test how incongruence varied independent of gene alignment length , we also generated a data matrix that was comprised of only the first 999 bp of sequence from each gene . since there were 667 genes shorter than 999 bp , this data matrix contains 1,340 genes from 20 species , every one of which has a 999 bp long alignment and does not contain any missing data . for the codon sequence alignment of each gene , the un - rooted phylogenetic tree under the optimality criterion of ml was inferred using raxml , version 8.0.20 ( stamatakis 2014 ) , under the gtrgamma model and with the values of the nucleotide base frequencies fixed to  observed  and those of the substitution rate parameters estimated from the data ( raxmlhpc - pthreads - sse3 -t 8 -f a -x 12345 -p 12345 -n 100 -m gtrgamma -s alignment -n name ) . for the concatenation analysis , the un - rooted concatenation species phylogeny was inferred through a single ml search under the gtrgamma model in raxml , version 8.0.20 ( stamatakis 2014 ) , with the values of the nucleotide base frequencies fixed to  observed  and those of the substitution rate parameters estimated from the data ( raxmlhpc - pthreads - sse3 -t 8 -f a -x 12345 -p 12345 -n 100 -m gtrgamma -s alignment -n name ) . in all cases , note that the raxml software first infers the topologies for each of the bootstrap replicates and then searches for the best - scoring ml tree using every fifth bootstrap replicate tree as a starting tree . the un - rooted emrc phylogeny that consisted of those bipartitions that appear in more than half of the ml gene trees , as well as of additional compatible bipartitions that appear in less than half of the gene trees , was inferred from the consense program in the phylogeny inference package , version 3.696 ( phylip ; j. felsenstein , university of washington , seattle ; http://evolution.genetics.washington.edu/phylip.html ) . the emrc phylogeny of bipartitions with high bs was constructed using custom perl scripts . as the divergence of culicinae and anophelinae lineages is well established , all phylogenies shown in figures have been mid - point rooted at the internode that separates these two lineages for easier visualization . the coalescent species phylogeny was estimated using 100 replicates of multi - locus bootstrapping in astral ( mirarab and warnow 2015 ) ( java -xmx36000 m -jar astral.4.7.8.jar -i treecollection -o output -b bs_path -r 100 ) , and using the online version of the star software with 100 rooted bootstrap replicates trees of every gene ( liu et al . distances between trees were estimated using the normalized robinson  foulds tree distance ( robinson and foulds 1981 ) as calculated by raxml , version 8.0.20 ( stamatakis 2014 ) ( raxmlhpc - pthreads - sse3 -t 8 -f r -z treecollection -m gtrgamma -n name ) . ic , tc , and relative tc ( rtc ) ( salichos and rokas 2013 ; salichos et al . 2014 ) were calculated using raxml , version 8.0.20 ( stamatakis 2014 ) ( raxmlhpc - pthreads - sse3 -t 8 -f i -t referencetree -z treecollection -m protgammaauto -n name ) . we used these average bs and tc values to construct eight subsets of 2,007 orthogroups : four with genes that have average bs values greater than or equal to 70% ( 1,818 genes ) , 80% ( 1,379 genes ) , 90% ( 378 genes ) , or 95% ( 66 genes ) , respectively , as well as four data sets comprising the 1,818 , 1,379 , 378 , or 66 genes whose ml trees have the highest tc values , respectively . for every gene from the 2,007-gene data matrix , we also extracted all bipartitions from its ml tree that have bs values greater than or equal to 70% , 80% , 90% , and 95% , respectively . we then used each one of these four sets of highly supported bipartitions to construct emrc species phylogenies with custom perl scripts . for the 1,340-gene data matrix , in which each gene s alignment was 999 bp long , we constructed six subsets : three with genes that have average bs values that are greater than or equal to 70% ( 1,138 genes ) , 80% ( 603 genes ) , or 90% ( 45 genes ) , respectively , ( no gene had average bs values greater than or equal to 95% ) ; as well as three data sets comprising the 1,138 , 603 , or 45 genes whose ml trees have the highest tc values , respectively . for every gene from the 1,340-gene data matrix , we also extracted all bipartitions from its ml tree that have bs values greater than or equal to 70% , 80% , 90% , and 95% . we then used each one of these four sets of highly supported bipartitions to construct emrc species phylogenies with custom perl scripts . we used the complete sets of annotated orthology data of 20 culicidae mosquito genomes ( supplementary table s1 , supplementary material online ) from http://cegg.unige.ch/orthodbmoz2 ( neafsey et al . we selected 2,008 single - copy genes which that contained sequences from all 20 species as our initial data set . we aligned all genes using the mafft software , version 7.182 ( katoh and toh 2008 ) based on their amino acid sequence , using the default settings ( automatic selection of the appropriate strategy , from l - ins - i , fft - ns - i , and fft - ns-2 , according to data size ) . 2006 ) to translate amino acid sequence alignments to codon sequence alignments , and the  automated  option of trimal ( capella - gutierrez et al . 2009 ) to trim the amino acid sequence alignments . trimmed segments of the amino acid sequence alignments were deleted from their corresponding codon sequence alignments using custom perl scripts . following trimming , our data matrix consisted of 2,007 genes from 20 species and contained no missing data . to test how incongruence varied independent of gene alignment length , we also generated a data matrix that was comprised of only the first 999 bp of sequence from each gene . since there were 667 genes shorter than 999 bp , this data matrix contains 1,340 genes from 20 species , every one of which has a 999 bp long alignment and does not contain any missing data . for the codon sequence alignment of each gene , the un - rooted phylogenetic tree under the optimality criterion of ml was inferred using raxml , version 8.0.20 ( stamatakis 2014 ) , under the gtrgamma model and with the values of the nucleotide base frequencies fixed to  observed  and those of the substitution rate parameters estimated from the data ( raxmlhpc - pthreads - sse3 -t 8 -f a -x 12345 -p 12345 -n 100 -m gtrgamma -s alignment -n name ) . for the concatenation analysis , the un - rooted concatenation species phylogeny was inferred through a single ml search under the gtrgamma model in raxml , version 8.0.20 ( stamatakis 2014 ) , with the values of the nucleotide base frequencies fixed to  observed  and those of the substitution rate parameters estimated from the data ( raxmlhpc - pthreads - sse3 -t 8 -f a -x 12345 -p 12345 -n 100 -m gtrgamma -s alignment -n name ) . in all cases , robustness in inference was assessed via bootstrap resampling ( 100 replicates ) . note that the raxml software first infers the topologies for each of the bootstrap replicates and then searches for the best - scoring ml tree using every fifth bootstrap replicate tree as a starting tree . the un - rooted emrc phylogeny that consisted of those bipartitions that appear in more than half of the ml gene trees , as well as of additional compatible bipartitions that appear in less than half of the gene trees , was inferred from the consense program in the phylogeny inference package , version 3.696 ( phylip ; j. felsenstein , university of washington , seattle ; http://evolution.genetics.washington.edu/phylip.html ) . the emrc phylogeny of bipartitions with high bs was constructed using custom perl scripts . as the divergence of culicinae and anophelinae lineages is well established , all phylogenies shown in figures have been mid - point rooted at the internode that separates these two lineages for easier visualization . the coalescent species phylogeny was estimated using 100 replicates of multi - locus bootstrapping in astral ( mirarab and warnow 2015 ) ( java -xmx36000 m -jar astral.4.7.8.jar -i treecollection -o output -b bs_path -r 100 ) , and using the online version of the star software with 100 rooted bootstrap replicates trees of every gene ( liu et al . distances between trees were estimated using the normalized robinson  foulds tree distance ( robinson and foulds 1981 ) as calculated by raxml , version 8.0.20 ( stamatakis 2014 ) ( raxmlhpc - pthreads - sse3 -t 8 -f r -z treecollection -m gtrgamma -n name ) . ic , tc , and relative tc ( rtc ) ( salichos and rokas 2013 ; salichos et al . 2014 ) were calculated using raxml , version 8.0.20 ( stamatakis 2014 ) ( raxmlhpc - pthreads - sse3 -t 8 -f i -t referencetree -z treecollection -m protgammaauto -n name ) . we used these average bs and tc values to construct eight subsets of 2,007 orthogroups : four with genes that have average bs values greater than or equal to 70% ( 1,818 genes ) , 80% ( 1,379 genes ) , 90% ( 378 genes ) , or 95% ( 66 genes ) , respectively , as well as four data sets comprising the 1,818 , 1,379 , 378 , or 66 genes whose ml trees have the highest tc values , respectively . for every gene from the 2,007-gene data matrix , we also extracted all bipartitions from its ml tree that have bs values greater than or equal to 70% , 80% , 90% , and 95% , respectively . we then used each one of these four sets of highly supported bipartitions to construct emrc species phylogenies with custom perl scripts . for the 1,340-gene data matrix , in which each gene s alignment was 999 bp long , we constructed six subsets : three with genes that have average bs values that are greater than or equal to 70% ( 1,138 genes ) , 80% ( 603 genes ) , or 90% ( 45 genes ) , respectively , ( no gene had average bs values greater than or equal to 95% ) ; as well as three data sets comprising the 1,138 , 603 , or 45 genes whose ml trees have the highest tc values , respectively . for every gene from the 1,340-gene data matrix , we also extracted all bipartitions from its ml tree that have bs values greater than or equal to 70% , 80% , 90% , and 95% . we then used each one of these four sets of highly supported bipartitions to construct emrc species phylogenies with custom perl scripts . all data and analyses described in this study are deposited at figshare under the accession 10.6084/m9.figshare.1566851 . supplementary figure s1 and tables s1s3 are available at genome biology and evolution online ( http://www.gbe.oxfordjournals.org/ ) .
comparison of individual gene trees in several recent phylogenomic studies from diverse lineages has revealed a surprising amount of topological conflict or incongruence , but we still know relatively little about its distribution across the tree of life . to further our understanding of incongruence , the factors that contribute to it and how it can be ameliorated , we examined its distribution in a clade of 20 culicidae mosquito species through the reconstruction and analysis of the phylogenetic histories of 2,007 groups of orthologous genes . 
 levels of incongruence were generally low , the three exceptions being the internodes concerned with the branching of anopheles christyi , with the branching of the subgenus anopheles as well as the already reported incongruence within the anopheles gambiae species complex . 
 two of these incongruence events ( a. gambiae species complex and a. christyi ) are likely due to biological factors , whereas the third ( subgenus anopheles ) is likely due to analytical factors . similar to previous studies , 
 the use of genes or internodes with high bootstrap support or internode certainty values , both of which were positively correlated with gene alignment length , substantially reduced the observed incongruence . 
 however , the clade support values of the internodes concerned with the branching of the subgenus anopheles as well as within the a. gambiae species complex remained very low . based on these results , we infer that the prevalence of incongruence in culicidae mosquitoes is generally low , that it likely stems from both analytical and biological factors , and that it can be ameliorated through the selection of genes with strong phylogenetic signal . 
 more generally , selection of genes with strong phylogenetic signal may be a general empirical solution for reducing incongruence and increasing the robustness of inference in phylogenomic studies .