Reusable tutorials for using cloud-based computing environments for the analysis of bacterial gene expression data from bulk RNA sequencing.

This manuscript describes the development of a resource module that is part of a learning platform named "NIGMS Sandbox for Cloud-based Learning" https://github.com/NIGMS/NIGMS-Sandbox. The overall genesis of the Sandbox is described in the editorial NIGMS Sandbox at the beginning of this Supplement. This module delivers learning materials on RNA sequencing (RNAseq) data analysis in an interactive format that uses appropriate cloud resources for data access and analyses. Biomedical research is increasingly data-driven, and dependent upon data management and analysis methods that facilitate rigorous, robust, and reproducible research. Cloud-based computing resources provide opportunities to broaden the application of bioinformatics and data science in research. Two obstacles for researchers, particularly those at small institutions, are: (i) access to bioinformatics analysis environments tailored to their research; and (ii) training in how to use Cloud-based computing resources. We developed five reusable tutorials for bulk RNAseq data analysis to address these obstacles. Using Jupyter notebooks run on the Google Cloud Platform, the tutorials guide the user through a workflow featuring an RNAseq dataset from a study of prophage altered drug resistance in Mycobacterium chelonae. The first tutorial uses a subset of the data so users can learn analysis steps rapidly, and the second uses the entire dataset. Next, a tutorial demonstrates how to analyze the read count data to generate lists of differentially expressed genes using R/DESeq2. Additional tutorials generate read counts using the Snakemake workflow manager and Nextflow with Google Batch. All tutorials are open-source and can be used as templates for other analysis.

A Uniform System for the Annotation of Vertebrate microRNA Genes and the Evolution of the Human microRNAome.

Although microRNAs (miRNAs) are among the most intensively studied molecules of the past 20 years, determining what is and what is not a miRNA has not been straightforward. Here, we present a uniform system for the annotation and nomenclature of miRNA genes. We show that less than a third of the 1,881 human miRBase entries, and only approximately 16% of the 7,095 metazoan miRBase entries, are robustly supported as miRNA genes. Furthermore, we show that the human repertoire of miRNAs has been shaped by periods of intense miRNA innovation and that mature gene products show a very different tempo and mode of sequence evolution than star products. We establish a new open access database--MirGeneDB ( http://mirgenedb.org )--to catalog this set of miRNAs, which complements the efforts of miRBase but differs from it by annotating the mature versus star products and by imposing an evolutionary hierarchy upon this curated and consistently named repertoire.

Widespread genomic breaks generated by activation-induced cytidine deaminase are prevented by homologous recombination.

Activation-induced cytidine deaminase (AID) is required for somatic hypermutation and immunoglobulin class switching in activated B cells. Because AID has no known target-site specificity, there have been efforts to identify non-immunoglobulin AID targets. We show here that AID acts promiscuously, generating widespread DNA double-strand breaks (DSBs), genomic instability and cytotoxicity in B cells with less homologous recombination ability. We demonstrate that the homologous-recombination factor XRCC2 suppressed AID-induced off-target DSBs, promoting B cell survival. Finally, we suggest that aberrations that affect human chromosome 7q36, including XRCC2, correlate with genomic instability in B cell cancers. Our findings demonstrate that AID has promiscuous genomic DSB-inducing activity, identify homologous recombination as a safeguard against off-target AID action, and have implications for genomic instability in B cell cancers.

Cetylpyridinium chloride (CPC) reduces zebrafish mortality from influenza infection: Super-resolution microscopy reveals CPC interference with multiple protein interactions with phosphatidylinositol 4,5-bisphosphate in immune function.

The COVID-19 pandemic raises significance for a potential influenza therapeutic compound, cetylpyridinium chloride (CPC), which has been extensively used in personal care products as a positively-charged quaternary ammonium antibacterial agent. CPC is currently in clinical trials to assess its effects on severe acute respiratory syndrome coronavirus 2 (SARS-CoV-2) morbidity. Two published studies have provided mouse and human data indicating that CPC may alleviate influenza infection, and here we show that CPC (0.1 µM, 1 h) reduces zebrafish mortality and viral load following influenza infection. However, CPC mechanisms of action upon viral-host cell interaction are currently unknown. We have utilized super-resolution fluorescence photoactivation localization microscopy to probe the mode of CPC action. Reduction in density of influenza viral protein hemagglutinin (HA) clusters is known to reduce influenza infectivity: here, we show that CPC (at non-cytotoxic doses, 5-10 µM) reduces HA density and number of HA molecules per cluster within the plasma membrane of NIH-3T3 mouse fibroblasts. HA is known to colocalize with the negatively-charged mammalian lipid phosphatidylinositol 4,5-bisphosphate (PIP2); here, we show that nanoscale co-localization of HA with the PIP2-binding Pleckstrin homology (PH) reporter in the plasma membrane is diminished by CPC. CPC also dramatically displaces the PIP2-binding protein myristoylated alanine-rich C-kinase substrate (MARCKS) from the plasma membrane of rat RBL-2H3 mast cells; this disruption of PIP2 is correlated with inhibition of mast cell degranulation. Together, these findings offer a PIP2-focused mechanism underlying CPC disruption of influenza and suggest potential pharmacological use of this drug as an influenza therapeutic to reduce global deaths from viral disease.

Rearrangement in the Hypervariable Region of JC Polyomavirus Genomes Isolated from Patient Samples and Impact on Transcription Factor-Binding Sites and Disease Outcomes.

JC polyomavirus (JCPyV) is the causative agent of the fatal, incurable, neurological disease, progressive multifocal leukoencephalopathy (PML). The virus is present in most of the adult population as a persistent, asymptotic infection in the kidneys. During immunosuppression, JCPyV reactivates and invades the central nervous system. A main predictor of disease outcome is determined by mutations within the hypervariable region of the viral genome. In patients with PML, JCPyV undergoes genetic rearrangements in the noncoding control region (NCCR). The outcome of these rearrangements influences transcription factor binding to the NCCR, orchestrating viral gene transcription. This study examines 989 NCCR sequences from patient isolates deposited in GenBank to determine the frequency of mutations based on patient isolation site and disease status. The transcription factor binding sites (TFBS) were also analyzed to understand how these rearrangements could influence viral transcription. It was determined that the number of TFBS was significantly higher in PML samples compared to non-PML samples. Additionally, TFBS that could promote JCPyV infection were more prevalent in samples isolated from the cerebrospinal fluid compared to other locations. Collectively, this research describes the extent of mutations in the NCCR that alter TFBS and how they correlate with disease outcome.

The Comparative Toxicogenomics Database: update 2017.

The Comparative Toxicogenomics Database (CTD; http://ctdbase.org/) provides information about interactions between chemicals and gene products, and their relationships to diseases. Core CTD content (chemical-gene, chemical-disease and gene-disease interactions manually curated from the literature) are integrated with each other as well as with select external datasets to generate expanded networks and predict novel associations. Today, core CTD includes more than 30.5 million toxicogenomic connections relating chemicals/drugs, genes/proteins, diseases, taxa, Gene Ontology (GO) annotations, pathways, and gene interaction modules. In this update, we report a 33% increase in our core data content since 2015, describe our new exposure module (that harmonizes exposure science information with core toxicogenomic data) and introduce a novel dataset of GO-disease inferences (that identify common molecular underpinnings for seemingly unrelated pathologies). These advancements centralize and contextualize real-world chemical exposures with molecular pathways to help scientists generate testable hypotheses in an effort to understand the etiology and mechanisms underlying environmentally influenced diseases.

A role for the CXCR4-CXCL12 axis in the little skate, Leucoraja erinacea.

The interaction between C-X-C chemokine receptor type 4 (CXCR4) and its cognate ligand C-X-C motif chemokine ligand 12 (CXCL12) plays a critical role in regulating hematopoietic stem cell activation and subsequent cellular mobilization. Extensive studies of these genes have been conducted in mammals, but much less is known about the expression and function of CXCR4 and CXCL12 in non-mammalian vertebrates. In the present study, we identify simultaneous expression of CXCR4 and CXCL12 orthologs in the epigonal organ (the primary hematopoietic tissue) of the little skate, Leucoraja erinacea. Genetic and phylogenetic analyses were functionally supported by significant mobilization of leukocytes following administration of Plerixafor, a CXCR4 antagonist and clinically important drug. Our results provide evidence that, as in humans, Plerixafor disrupts CXCR4/CXCL12 binding in the little skate, facilitating release of leukocytes into the bloodstream. Our study illustrates the value of the little skate as a model organism, particularly in studies of hematopoiesis and potentially for preclinical research on hematological and vascular disorders.

The Comparative Toxicogenomics Database's 10th year anniversary: update 2015.

Ten years ago, the Comparative Toxicogenomics Database (CTD; http://ctdbase.org/) was developed out of a need to formalize, harmonize and centralize the information on numerous genes and proteins responding to environmental toxic agents across diverse species. CTD's initial approach was to facilitate comparisons of nucleotide and protein sequences of toxicologically significant genes by curating these sequences and electronically annotating them with chemical terms from their associated references. Since then, however, CTD has vastly expanded its scope to robustly represent a triad of chemical-gene, chemical-disease and gene-disease interactions that are manually curated from the scientific literature by professional biocurators using controlled vocabularies, ontologies and structured notation. Today, CTD includes 24 million toxicogenomic connections relating chemicals/drugs, genes/proteins, diseases, taxa, phenotypes, Gene Ontology annotations, pathways and interaction modules. In this 10th year anniversary update, we outline the evolution of CTD, including our increased data content, new 'Pathway View' visualization tool, enhanced curation practices, pilot chemical-phenotype results and impending exposure data set. The prototype database originally described in our first report has transformed into a sophisticated resource used actively today to help scientists develop and test hypotheses about the etiologies of environmentally influenced diseases.

Natural selection canalizes expression variation of environmentally induced plasticity-enabling genes.

Many organisms survive fluctuating and extreme environmental conditions by manifesting multiple distinct phenotypes during adulthood by means of developmental processes that enable phenotypic plasticity. We report on the discovery of putative plasticity-enabling genes that are involved in transforming the gill of the euryhaline teleost fish, Fundulus heteroclitus, from its freshwater to its seawater gill-type, a process that alters both morphology and function. Gene expression that normally enables osmotic plasticity is inhibited by arsenic. Gene sets defined by antagonistic interactions between arsenic and salinity show reduced transcriptional variation among individual fish, suggesting unusually accurate and precise regulatory control of these genes, consistent with the hypothesis that they participate in a canalized developmental response. We observe that natural selection acts to preserve canalized gene expression in populations of killifish that are most tolerant to abrupt salinity change and that these populations show the least variability in their transcription of genes enabling plasticity of the gill. We found that genes participating in this highly canalized and conserved plasticity-enabling response had significantly fewer and less complex associations with transcriptional regulators than genes that respond only to arsenic or salinity. Collectively these findings, which are drawn from the relationships between environmental challenge, plasticity, and canalization among populations, suggest that the selective processes that facilitate phenotypic plasticity do so by targeting the regulatory networks that gives rise to the response. These findings also provide a generalized, conceptual framework of how genes might interact with the environment and evolve toward the development of plastic traits.

The Comparative Toxicogenomics Database: update 2013.

The Comparative Toxicogenomics Database (CTD; http://ctdbase.org/) provides information about interactions between environmental chemicals and gene products and their relationships to diseases. Chemical-gene, chemical-disease and gene-disease interactions manually curated from the literature are integrated to generate expanded networks and predict many novel associations between different data types. CTD now contains over 15 million toxicogenomic relationships. To navigate this sea of data, we added several new features, including DiseaseComps (which finds comparable diseases that share toxicogenomic profiles), statistical scoring for inferred gene-disease and pathway-chemical relationships, filtering options for several tools to refine user analysis and our new Gene Set Enricher (which provides biological annotations that are enriched for gene sets). To improve data visualization, we added a Cytoscape Web view to our ChemComps feature, included color-coded interactions and created a 'slim list' for our MEDIC disease vocabulary (allowing diseases to be grouped for meta-analysis, visualization and better data management). CTD continues to promote interoperability with external databases by providing content and cross-links to their sites. Together, this wealth of expanded chemical-gene-disease data, combined with novel ways to analyze and view content, continues to help users generate testable hypotheses about the molecular mechanisms of environmental diseases.

miRNAs: small genes with big potential in metazoan phylogenetics.

microRNAs (miRNAs) are a key component of gene regulatory networks and have been implicated in the regulation of virtually every biological process found in multicellular eukaryotes. What makes them interesting from a phylogenetic perspective is the high conservation of primary sequence between taxa, their accrual in metazoan genomes through evolutionary time, and the rarity of secondary loss in most metazoan taxa. Despite these properties, the use of miRNAs as phylogenetic markers has not yet been discussed within a clear conceptual framework. Here we highlight five properties of miRNAs that underlie their utility in phylogenetics: 1) The processes of miRNA biogenesis enable the identification of novel miRNAs without prior knowledge of sequence; 2) The continuous addition of miRNA families to metazoan genomes through evolutionary time; 3) The low level of secondary gene loss in most metazoan taxa; 4) The low substitution rate in the mature miRNA sequence; and 5) The small probability of convergent evolution of two miRNAs. Phylogenetic analyses using both Bayesian and parsimony methods on a eumetazoan miRNA data set highlight the potential of miRNAs to become an invaluable new tool, especially when used as an additional line of evidence, to resolve previously intractable nodes within the tree of life.

Toward consilience in reptile phylogeny: miRNAs support an archosaur, not lepidosaur, affinity for turtles.

Understanding the phylogenetic position of crown turtles (Testudines) among amniotes has been a source of particular contention. Recent morphological analyses suggest that turtles are sister to all other reptiles, whereas the vast majority of gene sequence analyses support turtles as being inside Diapsida, and usually as sister to crown Archosauria (birds and crocodilians). Previously, a study using microRNAs (miRNAs) placed turtles inside diapsids, but as sister to lepidosaurs (lizards and Sphenodon) rather than archosaurs. Here, we test this hypothesis with an expanded miRNA presence/absence dataset, and employ more rigorous criteria for miRNA annotation. Significantly, we find no support for a turtle + lepidosaur sister-relationship; instead, we recover strong support for turtles sharing a more recent common ancestor with archosaurs. We further test this result by analyzing a super-alignment of precursor miRNA sequences for every miRNA inferred to have been present in the most recent common ancestor of tetrapods. This analysis yields a topology that is fully congruent with our presence/absence analysis; our results are therefore in accordance with most gene sequence studies, providing strong, consilient molecular evidence from diverse independent datasets regarding the phylogenetic position of turtles.

Color-Flu Fluorescent Reporter Influenza A Viruses Allow for In Vivo Studies of Innate Immune Function in Zebrafish.

Influenza virus infection can cause severe respiratory disease and is estimated to cause millions of illnesses annually. Studies on the contribution of the innate immune response to influenza A virus (IAV) to viral pathogenesis may yield new antiviral strategies. Zebrafish larvae are useful models for studying the innate immune response to pathogens, including IAV, in vivo. Here, we demonstrate how Color-flu, four fluorescent IAV strains originally developed for mice, can be used to study the host response to infection by simultaneously monitoring infected cells, neutrophils, and macrophages in vivo. Using this model, we show how the angiotensin-converting enzyme inhibitor, ramipril, and mitophagy inhibitor, MDIVI-1, improved survival, decreased viral burden, and improved the respiratory burst response to IAV infection. The Color-flu zebrafish larvae model of IAV infection is complementary to other models where the dynamics of infection and the response of innate immune cells can be visualized in a transparent host in vivo.

Remarkably High Repeat Content in the Genomes of Sparrows: The Importance of Genome Assembly Completeness for Transposable Element Discovery.

Transposable elements (TE) play critical roles in shaping genome evolution. Highly repetitive TE sequences are also a major source of assembly gaps making it difficult to fully understand the impact of these elements on host genomes. The increased capacity of long-read sequencing technologies to span highly repetitive regions promises to provide new insights into patterns of TE activity across diverse taxa. Here we report the generation of highly contiguous reference genomes using PacBio long-read and Omni-C technologies for three species of Passerellidae sparrow. We compared these assemblies to three chromosome-level sparrow assemblies and nine other sparrow assemblies generated using a variety of short- and long-read technologies. All long-read based assemblies were longer (range: 1.12 to 1.41 Gb) than short-read assemblies (0.91 to 1.08 Gb) and assembly length was strongly correlated with the amount of repeat content. Repeat content for Bell's sparrow (31.2% of genome) was the highest level ever reported within the order Passeriformes, which comprises over half of avian diversity. The highest levels of repeat content (79.2% to 93.7%) were found on the W chromosome relative to other regions of the genome. Finally, we show that proliferation of different TE classes varied even among species with similar levels of repeat content. These patterns support a dynamic model of TE expansion and contraction even in a clade where TEs were once thought to be fairly depauperate and static. Our work highlights how the resolution of difficult-to-assemble regions of the genome with new sequencing technologies promises to transform our understanding of avian genome evolution.

MicroRNAs support the monophyly of enteropneust hemichordates.

Understanding the evolutionary history of deuterostomes requires elucidating the phylogenetic interrelationships amongst the constituent taxa. Although the monophyly and interrelationships among the three principal groups-the chordates, the echinoderms, and the hemichordates-are well established, as are the internal relationships among the echinoderm and chordate taxa, the interrelationships among the principal groups of hemichordates-the harrimaniid enteropneusts, the ptychoderid enteropneusts, and the pterobranchs-remain unresolved. Depending on the study some find enteropneusts paraphyletic with pterobranchs (e.g., Cephalodiscus) more closely related to the harrimaniid enteropneusts (e.g., Saccoglossus) than either are to the ptychoderid enteropneusts (e.g., Ptychodera), whereas other studies support a monophyletic Enteropneusta. To try and resolve between these two competing hypotheses, we turned to microRNAs, small ∼22 nt non-coding RNA genes that have been shown to shed insight into particularly difficult phylogenetic questions. Using deep sequencing we characterized the small RNA repertoires of two hemichordate species, Cephalodiscus hodgsoni and Ptychodera flava, and the crinoid echinoderm Antedon mediterranea, and combined our results with the described complements of the hemichordate Saccoglossus kowalevskii, the sea urchin Strongylocentrotus purpuratus, and the starfish Patiria miniata. Our data unambiguously support the monophyly of Enteropneusts as S. kowalevskii shares 12 miRNA sequences with P. flava that are not present in the C. hodgsoni or A. mediterranea libraries, and have never been reported from another metazoan taxon. Thus, these data resolve the phylogenetic position of pterobranchs, ultimately allowing for a better understanding of body plan evolution throughout the deuterostomes.

The Comparative Toxicogenomics Database: update 2011.

The Comparative Toxicogenomics Database (CTD) is a public resource that promotes understanding about the interaction of environmental chemicals with gene products, and their effects on human health. Biocurators at CTD manually curate a triad of chemical-gene, chemical-disease and gene-disease relationships from the literature. These core data are then integrated to construct chemical-gene-disease networks and to predict many novel relationships using different types of associated data. Since 2009, we dramatically increased the content of CTD to 1.4 million chemical-gene-disease data points and added many features, statistical analyses and analytical tools, including GeneComps and ChemComps (to find comparable genes and chemicals that share toxicogenomic profiles), enriched Gene Ontology terms associated with chemicals, statistically ranked chemical-disease inferences, Venn diagram tools to discover overlapping and unique attributes of any set of chemicals, genes or disease, and enhanced gene pathway data content, among other features. Together, this wealth of expanded chemical-gene-disease data continues to help users generate testable hypotheses about the molecular mechanisms of environmental diseases. CTD is freely available at http://ctd.mdibl.org.

Multi-spectral Fluorescent Reporter Influenza A Viruses Allow for in vivo Studies of Innate Immune Function in Zebrafish.

Influenza virus infection can cause severe respiratory disease and is estimated to cause millions of illnesses annually. Studies of the contribution of the innate immune response to influenza A virus (IAV) to viral pathogenesis may yield new antiviral strategies. Zebrafish larvae are useful models to study the innate immune response to pathogens, including IAV, in vivo. Here, we demonstrate how Color-flu, four fluorescent IAV strains originally developed for mice, can be used to study host-virus interactions by simultaneously monitoring virus particles, neutrophils, and macrophages in vivo. Using this model, we show how the angiotensin-converting enzyme inhibitor, ramipril, and mitophagy inhibitor, MDIVI-1, improved survival, decreased viral burden, and improved the respiratory burst response to IAV infection. The Color-flu zebrafish model of IAV infection is complementary to other models as it is the only model where interactions between virus particles and host cells in an intact vertebrate can be visualized in vivo.

A New Phenotype in Candida-Epithelial Cell Interaction Distinguishes Colonization- versus Vulvovaginal Candidiasis-Associated Strains.

Vulvovaginal candidiasis (VVC) affects nearly 3/4 of women during their lifetime, and its symptoms seriously reduce quality of life. Although Candida albicans is a common commensal, it is unknown if VVC results from a switch from a commensal to pathogenic state, if only some strains can cause VVC, and/or if there is displacement of commensal strains with more pathogenic strains. We studied a set of VVC and colonizing C. albicans strains to identify consistent in vitro phenotypes associated with one group or the other. We find that the strains do not differ in overall genetic profile or behavior in culture media (i.e., multilocus sequence type [MLST] profile, rate of growth, and filamentation), but they show strikingly different behaviors during their interactions with vaginal epithelial cells. Epithelial infections with VVC-derived strains yielded stronger fungal proliferation and shedding of fungi and epithelial cells. Transcriptome sequencing (RNA-seq) analysis of representative epithelial cell infections with selected pathogenic or commensal isolates identified several differentially activated epithelial signaling pathways, including the integrin, ferroptosis, and type I interferon pathways; the latter has been implicated in damage protection. Strikingly, inhibition of type I interferon signaling selectively increases fungal shedding of strains in the colonizing cohort, suggesting that increased shedding correlates with lower interferon pathway activation. These data suggest that VVC strains may intrinsically have enhanced pathogenic potential via differential elicitation of epithelial responses, including the type I interferon pathway. Therefore, it may eventually be possible to evaluate pathogenic potential in vitro to refine VVC diagnosis. IMPORTANCE Despite a high incidence of VVC, we still have a poor understanding of this female-specific disease whose negative impact on women's quality of life has become a public health issue. It is not yet possible to determine by genotype or laboratory phenotype if a given Candida albicans strain is more or less likely to cause VVC. Here, we show that Candida strains causing VVC induce more fungal shedding from epithelial cells than strains from healthy women. This effect is also accompanied by increased epithelial cell detachment and differential activation of the type I interferon pathway. These distinguishing phenotypes suggest it may be possible to evaluate the VVC pathogenic potential of fungal isolates. This would permit more targeted antifungal treatments to spare commensals and could allow for displacement of pathogenic strains with nonpathogenic colonizers. We expect these new assays to provide a more targeted tool for identifying fungal virulence factors and epithelial responses that control fungal vaginitis.

Prefoldin 5 is required for normal sensory and neuronal development in a murine model.

Molecular chaperones and co-chaperones are crucial for cellular development and maintenance as they assist in protein folding and stabilization of unfolded or misfolded proteins. Prefoldin (PFDN), a ubiquitously expressed heterohexameric co-chaperone, is necessary for proper folding of nascent proteins, in particular, tubulin and actin. Here we show that a genetic disruption in the murine Pfdn5 gene, a subunit of prefoldin, causes a syndrome characterized by photoreceptor degeneration, central nervous system abnormalities, and male infertility. Our data indicate that a missense mutation in Pfdn5, may cause these phenotypes through a reduction in formation of microtubules and microfilaments, which are necessary for the development of cilia and cytoskeletal structures, respectively. The diversity of phenotypes demonstrated by models carrying mutations in different PFDN subunits suggests that each PFDN subunit must confer a distinct substrate specificity to the prefoldin holocomplex.

MicroRNAs support a turtle + lizard clade.

Despite much interest in amniote systematics, the origin of turtles remains elusive. Traditional morphological phylogenetic analyses place turtles outside Diapsida-amniotes whose ancestor had two fenestrae in the temporal region of the skull (among the living forms the tuatara, lizards, birds and crocodilians)-and allied with some unfenestrate-skulled (anapsid) taxa. Nonetheless, some morphological analyses place turtles within Diapsida, allied with Lepidosauria (tuatara and lizards). Most molecular studies agree that turtles are diapsids, but rather than allying them with lepidosaurs, instead place turtles near or within Archosauria (crocodilians and birds). Thus, three basic phylogenetic positions for turtles with respect to extant Diapsida are currently debated: (i) sister to Diapsida, (ii) sister to Lepidosauria, or (iii) sister to, or within, Archosauria. Interestingly, although these three alternatives are consistent with a single unrooted four-taxon tree for extant reptiles, they differ with respect to the position of the root. Here, we apply a novel molecular dataset, the presence versus absence of specific microRNAs, to the problem of the phylogenetic position of turtles and the root of the reptilian tree, and find that this dataset unambiguously supports a turtle + lepidosaur group. We find that turtles and lizards share four unique miRNA gene families that are not found in any other organisms' genome or small RNA library, and no miRNAs are found in all diapsids but not turtles, or in turtles and archosaurs but not in lizards. The concordance between our result and some morphological analyses suggests that there have been numerous morphological convergences and reversals in reptile phylogeny, including the loss of temporal fenestrae.