dataesr/scientific-paragraphs-categorization huggingface dataset:
We likewise observed that genes reported in D. melanogaster to have highest enrichment in the ovary largely correspond to genes that are significantly upregulated in the ovaries of Hawaiian species (Fig H in S1 Text). The 131 core ovary genes include several well-known members involved in oogenesis and germline stem cell renewal such as nanos  [34] , swallow  [35] , and oskar  [46] (Fig  2E ). We found only two genes that were identified as Hawaiian core ovary genes that are not reported in the FlyAtlas2 dataset  [29]  to be enriched in the ovary of D. melanogaster: the SET domain binding factor sbf, and Rfx, which are reported to be enriched in the heart, brain, and other non-reproductive tissues of D. melanogaster  [29]. We used the same approach to identify a core suite of 52 head-specific genes (Fig I in S1 Text).
There was no overlap between the sets of core head genes and core ovary genes. To test whether the correspondence between expression observations in Hawaiian flies and D. The x-axis shows the log 2 fold change of expression across transcripts, and points are arranged with random jitter within species on the y-axis. C, The same jitter plots as in B, now colored according to whether or not transcripts belong to a cohort of core ovary genes. These are defined as genes, grouped by BLAST similarity to D. melanogaster transcripts, for which at least one transcript is upregulated in the ovary of ten or more of the twelve species. D, A comparison of mean expression change across Hawaiian species to reported ovary-enrichment values from D. melanogaster, as reported in FlyAtlas2  [29]. Core ovary genes are marked in magenta.
These can then be divided into orthology groups by estimating gene trees and identifying maximally inclusive subtrees with no more than one sequence per taxon  [48]. We compared the representation of species across homology and orthology groups, and observed that while the representation of homology groups increases with the number of species compared, representation of orthology groups decreases  (Fig K in S1 Text). This is a known obstacle in comparative transcriptomics, attributed to many transcripts being artifactually fragmented during reference transcriptome assembly  [49]. To reduce the impact of this on our downstream analyses, we averaged TPM values across all transcripts within a homology group for each sequenced RNA library. Principle component analysis of this average expression dataset showed that the first principle component divides ovary libraries from the rest, while the second component separates samples along an axis that largely corresponds to phylogenetic distance between species (Fig L in S1 Text).
While this averaging approach reduces noise due to variable mapping affinities of fragments of the same transcript, it comes at the cost of averaging over potential variation between genuine transcripts that fall into the same homology group. Future analyses using improved assemblies for transcriptomes or genomes will likely be able to avoid this trade off and compare transcript counts directly. To test the robustness of results to this averaging, we also performed key analyses over a dataset of the identifiable strict orthologs. In both datasets (based on homology groups or based on strict orthology groups), we averaged expression values across biological replicate individuals per species. With average expression counts for homologous transcripts across species, we tested the degree to which variation across this dataset could be attributed to tissue-specific variation (here, ovary vs.
carcass), species-specific variation, or neither (residual variation). Using the linear modeling approach adapted from Breschi and colleagues (2016)  [21] , we found the proportion of variance across the dataset attributed to tissue differences decreased with phylogenetic distance, while the proportion attributed to species difference increased (Fig 3A  -3C ). In addition, we found that, when comparing ovary and carcass tissues, the Hawaiian drosophilid clade encompasses the crossover point where variation across species swamps variation across tissues (crossed lines, Fig  3A ). When comparing across the two species from the picture-wing group included in this study, an average of 45. 6% of the variation can be attributed to tissue differences. For the same comparison, 960 genes were identified as tissue-variable genes (TVGs), defined as residual variation accounting for <25% and a two-fold increase in variation attributed to tissues than to species (Fig 3B and M in S1 Text).
Visualizing the distribution of scaled changes by genes shows that most scaled changes are small and centered around zero, representing little change in gene expression bias between tissues (Fig 4B  -4C ). Using this dataset of scaled changes across genes and branches, we identified branches for which the direction of tissue bias had changed (e. g. from higher expression in the ovary than in the carcass to lower, or vice versa). Visualizing this dataset according to branches reveals that the majority of these changes in bias are located on the root and terminal branches, rather than internal branches (Fig 4D  and 4E ). This is likely because internal branches for this rapid radiation tend to be very short; even when scaling evolutionary changes to branch length, it is less probable for our analysis to identify a shift to and from ovary-biased expression on a short branch than a long branch.
Repeating the same analysis on a dataset of strict orthologs reflects the same pattern, indicating that this result is not an artefact of expression averaging across homology groups (Fig R in S1 Text). Visualizing the distribution of genes by ancestral and descendant values allows us to identify shifts in bias which represent the largest swings in expression values (Fig 4F  , points a-d ). Highlighting the top two such shifts in both directions, we identify four example genes which acquired or lost ovary-specific expression in the phylogeny of Hawaiian Drosophilidae. In the case of FMRFaR and GABA, a few Hawaiian species have gained ovary-biased expression of these genes, while most species and the ancestral state indicate non-ovary bias (Fig  4Fa-b ). In the case of vilya and the unnamed gene CG9109, each shows a pattern where one species has lost ovary bias from a biased ancestral state (Fig  4Fc-d ).
To test the degree to which the correlations observed here reflect known biological interactions between genes, we compared these measures to reported protein and genetic interactions between genes, using the database of published genetic experiments in D. melanogaster, available at http://flybase. org. We found that the mean correlation coefficient for genes that are known to physically interact as proteins was higher than for genes with no or unknown interaction (maximum p-value = <0. 001 over 100 replicates, Fig  5A ). This indicates that even with a relatively small number of observations, there is sufficient information in the matrix to detect biological signal between gene pairs. These results were calculated based on the correlation in expression bias between the ovary and carcass. However, following the same procedure using correlations in changes in head-biased expression showed no significant difference between the two groups (max.
p-value = 0. 256, Fig T in S1 Text), suggesting the strength of this signal may be dependent on the tissues being compared. We also found that genes known to interact genetically have a significantly higher mean correlation than genes with no or unknown genetic interactions (unknown vs. enhancement swings in expression within those quadrants are labeled a-d. G, The four genes with large swings from F, showing the expression bias for each transcript colored according to more expression in the ovary (red) or carcass (blue). Panels annotated with the gene symbol from the D. melanogaster sequences in the same homology group, with the exception of vilya � , which was annotated using a direct BLAST search since no D. melanogaster sequence was present in that group. https://doi. org/10. 1371/journal.
We found eight distinct homologous gene groups, comprising 12 unique D. melanogaster parent genes, that had a strong evolutionary correlation with yolk-protein genes (absolute coefficient greater than 0. 825, Fig  5D ). None of these correlated genes correspond to those listed on FlyBase  [52]  as having known interactions with yolk-protein genes in D. melanogaster (Fig  5E ). We consider these strong evolutionary correlations to be a set of new predictions about evolutionary and genetic relationships between genes which can be tested in wild and laboratory model species of Drosophila. The correlation networks for each gene family, estimated using the full dataset or using the more conservative dataset with low-confidence outliers removed, can be interrogated in the accompanying data visualization for this manuscript (https://shchurch.
However, it comes at the cost of averaging over possible biological variation in expression between genes from the same gene family. The concordance of our results with published findings from D. melanogaster suggests that the approach we have used here is robust for our dataset. However, as the quality and accessibility of genomes from diverse species continue to increase, future studies will likely be able to compare directly between orthologous genes without needing to account for fragmentation. For those future studies, a phylogenetic comparative approach like the one used here and elsewhere  [8]  can serve as an analytical framework to move expression comparisons beyond pairwise comparisons. One goal of evolutionary developmental biology is to identify changes in developmental mechanisms that underlie phenotypic differences  [12].
For DNA barcoding, DNA was extracted from one or two legs from male specimens using the Qiagen DNeasy blood and tissue extraction kit, or from the DNA of females isolated during RNA extraction (see below). We amplified and sequenced the cytochrome oxidase I (COI), II (COII) and 16S rRNA genes using the primers and protocols described in Sarikaya and colleagues (2019)  [40]. For barcode matching, we aligned sequences using MAFFT, version v7. 475  [64] , and assembled gene trees using RAxML, version 8. 2. 9  [65]. Definitive matches were considered when sequences for females formed a monophyletic clade with reference males or reference sequences from NCBI; see Table B in S1 Text. Female D. primaeva, D. macrothrix, D. sproati, and D. picticornis could be identified unambiguously using dichotomous keys.
Transcriptome assembling and expression mapping Transcriptome assembly and expression mapping was performed using the agalma pipeline, version 2. 0. 0  [48]. For the twelve reference transcriptomes, reads from separate rounds of sequencing were concatenated and inserted into the agalma catalog. Further details of transcriptome assembly and homology assessment are included in our previous manuscript  [44]. Each tissue-specific expression library was mapped to the corresponding reference transcriptome using the 'expression' pipeline in agalma, which uses the software RSEM to estimate gene and isoform count levels from RNAseq data  [66]. The agalma pipeline also includes steps to catalog the species, tissue type, and run information, which were exported as a single JavaScript object notation (JSON) file.
Phylogenetic analysis The methods for inferring a phylogenetic tree of the species here are the same as those described in our previous manuscript  [44]. In summary, we used the agalma pipeline  [48]  to infer homology, based on an all-by-all BLAST procedure. Gene orthology was inferred based on the topology of gene trees, estimated with RAxML  [65] , version 8. 2. 9. The most likely tree was inferred using IQtree  [67] , version 2. 1. 1, on a dataset partitioned by transcripts and using the default Model Finder  [68]. For this study, gene trees were additionally annotated with the software phyldog  [69].
Annotating transcripts by sequence similarity We leveraged the close relationship of these species to species of Drosophila with well-annotated genomes to annotate the transcripts considered here. For each transcript in the reference transcriptome, we performed four comparisons of sequence similarity using local BLAST:  [1]  comparing nucleotide transcript sequences to nucleotide sequences from D. melanogaster (blastn),  [2]  comparing translated nucleotide sequences to protein sequences of D. melanoagster (blastx),  [3]  comparing nucleotide sequences to a database of nucleotide sequences from D. melanogaster, D. virilis, and D. grimshawi (blastx), and  [4]  comparing translated nucleotide sequences to a database of protein sequences from the same three species (blastn). For downstream analyses, we prioritized annotations from the second comparison, but we provide all sequence similarity reports in the GitHub repository under the directory analysis/ BLAST.
Normalization and differential gene expression Transcript count tables were imported into R using the agalmar package, version 0. 0. 0. 9000. Differential gene expression analysis was performed using the package DESeq2, version 1. 34. 0. For these analyses we used only one sequencing run per library, thereby excluding duplicate sequencing runs. Analyses of differential gene expression were calculated using the default approaches in DESeq2 for estimating size factors, dispersions, and calculating log 2 fold-change and adjusted p-values using the Benjamini and Hochberg correction (Fig B, panel a, in S1 Text). Both individual and tissue were considered in the design formula. Transcripts were considered differentially expressed at a significance threshold of 0. 01. We identified a cohort of core ovary-specific genes by first identifying a parent gene for each transcript using a sequence similarity search against D.
Comparison of expression to D. melanogaster We compared our differential gene expression results to a reference database of tissue expression from D. melanogaster, known as the FlyAtlas2  [29]. We downloaded this reference in July of 2021, from http://motif. gla. ac. uk/downloads/FlyAtlas2_21. 04. 18. sql. This dataset provides data on transcript abundance and tissue enrichment, including for female ovaries. Tissue enrichment is calculated using the same methods as in the FlyAtlas2 web browser, defined as the fragments per kilobase of transcript per million mapped reads (FPKM) for a given tissue divided by that value for the reference tissue (here, female whole body), with a pseudocount of two counts added to empty values to avoid division by zero. We considered a FlyAtlas gene to be enriched in the ovary, comparable to our data, if the ovary was the maximum enrichment value across all tissues excluding the head, brain, and eye tissues, as these were separated in our RNASeq procedure (Fig B, panel a, in S1 Text).
Transforming data into comparable measurements of expression across species Transcript counts are reported in transcripts per million (TPM), but this measurement is known to not be directly comparable across species due to differences in reference transcriptome size  [7, 8]. Therefore, we normalized TPM by species using the procedure described by Munro and colleagues (2021)  [8] , where TPM values are multiplied by the number of genes in the reference, and this value is divided by 10 An additional challenge when working with reference transcriptomes is the presence of fragmented transcripts created during the assembly process  [49]. This fragmentation can result in noise in estimating the amount of transcript as reads are differentially mapped to these fragments. To reduce the impact of this noise on our analysis, we undertook a novel approach where transcripts were grouped according to inferred homology as estimated by the agalma pipeline using an all-by-all BLAST approach (Fig B, panel b, in S1 Text).
Linear modeling We performed linear modeling to calculate the relative contribution of tissue-and specieslevel differences to variation in gene expression (Fig B, panel b, in S1 Text), following the approach of Breschi and colleagues (2016)  [21]. These analyses were performed separately on datasets of ovary vs. carcass and head vs. carcass expression. Using the ANOVA script provided at https://github. com/abreschi/Rscripts/blob/master/anova. R, we built a linear model for each gene that accounts for the contribution of the organ, species, and any residual error. We then calculated the relative proportion of each factor divided by the total sum of squares for all factors. We identified groups of highly variable genes, using the same metrics defined by Breschi and colleagues (2016)  [21] , as any gene for which either tissues or species explains at least 75% of the variance.
We reconstructed the evolutionary history of tissue bias for each homology group using the species tree published in Church and Extavour, 2021  [44] , based on the same reference transcriptome data (Fig B  , panel c , in S1 Text). First, we calibrated the tree estimated using IQtree (Fig 1A of that publication) to be ultrametric using the R function chronos in the package ape, version 5. 6. 2 (using a correlated model and a lambda value of 1). We then subset this tree to only include tips for which expression data was available, and annotated this tree to be able to identify specific branches and nodes in ancestral state reconstruction analyses. Ancestral expression bias values were estimated with the R package Rphylopars, version 0. 3. 9, using the fast ancestral state reconstruction algorithm based on Ho and Ane, 2014  [70]  (Fig C, panel a, in S1 Text).
When sufficient specimens were available per species, one whole individual was used as a reference and three whole individuals were dissected into three separate tissues: the head, ovaries, and all remaining material (carcass). Reference individuals were sequenced to generate paired-end RNA reads and dissected tissues were sequenced to generate single-end RNA reads. Tissue libraries were then mapped to the assembled reference to quantify transcript expression. Teal boxes indicate data files. Dashed-line boxes indicate a repeated step.
These are defined as genes, grouped by BLAST similarity to D. melanogaster transcripts, for which at least one transcript is upregulated in the ovary of ten or more of the twelve species. D, A comparison of mean expression change across Hawaiian species to reported ovary-enrichment values from D. melanogaster, as reported in FlyAtlas2 [29]. Core ovary genes are marked in magenta. E, The boxed region shown in D, magnified and now showing only core ovary genes, annotated with the gene symbol from D. melanogaster.
Therefore, the LCO surface was not perfectly flat, and even if an optimal zone is chosen by mean of optical observation, a certain topography is expected on the analyzed surface. Scanning Micro-Raman Spectroscopy Data Processing: Raman spectra were corrected by fitting the baseline with a linear function and subtracting it to the signal (Figure  S10 , Supporting Information). Then the measured data were fitted by a sum of Laurentian functions as explained in details in Supporting Information (Figures  S11  and S12 , Supporting Information). Electrochemical Measurements: All electrochemical measurements were performed on a Biologic VMP3 Multichannel Potentiostat. The measurements were carried out in a Cascade SUMMIT200 Probe Station in dry air environment. Impedance spectroscopy measurement was performed between 1 MHz and 100 mHz.
The cell was stimulated with a 10 mV sinusoidal voltage and the current response was measured. Each measurement was averaged 5 times for each frequency. Charge protocol was a Constant Potential (CV) stimulation of the cell at 4. 3 V. The current was monitored and the charge stopped when 5% of the initial current was reached. Discharge protocol was a Constant Current (CC) stimulation. The discharge was stopped when the cell potential reaches 3 V. Galvanostatic Intermittent Titration Technique (GITT): GITT measurement was performed using a Biologic VMP3 Multichannel Potentiostat under argon atmosphere (Figure  S5 , Supporting Information). Galvanostatic Intermittent Titration Technique (GITT) Data Processing: GITT measurement data are treated as reported in Supporting Information (Figures S6 and S7, Supporting Information).
Identification of RNU4ATAC Mutations in Patients with JBTS-Like Traits. Two unrelated children suspected with having JBTS were referred to the genetic consultation. Targeted exome sequencing was performed, but no biallelic pathogenic variants or quantitative rearrangements were detected in known JBTS-related genes or other ciliary protein-encoding genes present on the Ciliome/CilDiag panels (SI Appendix, Materials and Methods). After modification of several bioinformatic parameters and the removal of a filter excluding noncoding genes, we identified in the two patients the same homozygous pathogenic variant, n. 16G>A, in RNU4ATAC (GenBank: NR_023343) (Fig. 1A  and SI Appendix, Table  S1 ). Patient P1 (F1:II-2) exhibits an MTS, a severe developmental delay with severe axial hypotonia, microcephaly, skeletal dysplasia, facial dysmorphism, retinal coloboma, nystagmus, eczema, bilateral postaxial polydactyly, genital anomalies, growth hormone deficiency, small interventricular septal defect, and immunodeficiency (SI Appendix, Table  S2 ).
Patient P2 (F2:II-3), who died at 10 d of age, showed an MTS, microcephaly, skeletal dysplasia, an atrial septal defect, hand and foot abnormalities and ocular nystagmus (and strabismus) (SI Appendix, Table  S2 ). Based on this discovery, we searched the presence of microcephaly and/or intrauterine growth retardation (IUGR) in a JBTS database of 183 patients without a molecular diagnosis after NGS analysis of a panel of 30 known JBTS genes. Six patients matching these criteria were identified and sequenced for RNU4ATAC. One of them presented with the compound heterozygous variants n. [16G>A];[51G>A] (Fig. 1A  and SI Appendix, Table  S1 ). This patient, patient P3 (F3:II-1), is a 4-y-old girl who presents with an MTS, severe developmental delay with severe hypotonia, microcephaly, antenatal growth retardation, nystagmus, dysmorphism, and skeletal dysplasia (SI Appendix, Table  S2 ).
Here, we used a different approach and performed a GO term enrichment analysis of the U12 genes vs. the U2 genes, using the entire set of genes with GO annotations (n = 618 and 19,939 for U12 and U2 genes, respectively). An enrichment is observed if more than 3% (618/20,557) of U12 genes are present in a GO term. Hence, we found 243 enriched biological processes and 96 enriched cellular components with a P value ≤0. 05 (Fig. 2A  and Datasets S1). The most enriched biological process is the membrane depolarization during action potential (P value = 5. 2e-15, with 42% of U12 genes), which partially overlaps the other GO terms neuronal action potential (P = 5. 5e-13, 40%), regulation of ion transmembrane transport (P = 1. 2e-12, 8%), calcium ion transmembrane transport (P = 4. 3e-6, 9%), and sodium ion transmembrane transport (P = 5.
2e-6, 10%). Interestingly, the nonmotile cilium assembly process (P = 6. 8e-6, 18%) was the thirteenth most significant term (red in Fig. 2A  and Dataset S1). Other significant terms explicitly related to cilia were identified, namely intraciliary transport involved in cilium assembly (P = 1. 2e-3, 15%), ciliary basal body-plasma membrane docking (P = 1. 6e-2, 7%), and protein localization to ciliary transition zone (P = 2. 8e-2, 22%). Similarly, GO term enrichment analysis of cellular components revealed an unexpected amount of U12 geneencoding proteins that localize at the basal body (P = 1. 4e-5, 11%), the ciliary base (P = 4. 3e-3, 14%), and the ciliary tip (P = 1. 0e-2, 11%) (Dataset S1). Of note, the opposite comparison of U2 vs. U12 genes did not reveal any cilium-related terms (Dataset S1).